function ind = compress2D(x, y, width, height, alpha)
% COMPRESS2D: Rudimental compression algorithm for 2D matlab plots
%
%
% It is a variant of compress2Dh(x, y, width, height, alpha). compress2Dh
% is MUCH slower and A BIT more accurate than compress2D.
% Use them if you need to plot highly populated vectors when you expect
% high density of points per pixel.
%
%
% IND = COMPRESS2D(x, y, width, height, alpha) reduces the number of points
% (x, y) which is necessary to represent in a figure, eliminating redundancies.
% It returns the indexes of vectors x and y which are to be plotted.
%
% It is necessary to specify the width and the height of the desired figure
% and a resolution parameter alpha. The plane is initially fragmented in a
% number of cells which is a function of width, height and alpha:
% ceil(alpha * width) * ceil(height / width * ceil(alpha * width))
%
% alpha must be a positive number such that width * alpha < 3000: the
% smaller alpha, the higher the compression will be, the lower quality the
% final image will have. A suitable value for alpha, when the figure is at
% its maximizum size, is 0.27.
%
% ATTENTION: alpha, width and height all contribute to the quality of the
% final, compressed image. Try different combinations in order to obtain
% the desired compression at the desired quality.
%
%
% If uncertain on how to specify width and height, these lines may help:
% h = figure; % open a new figure (manually choose the desired size)
% set(h, 'Units', 'pixel'); % unit must be pixels
% pos = get(h, 'Position'); % get the Position
% close; % close figure
% pos(1:2) = []; % keep only the last two
% width = pos(1); % find the desired value of width
% height = pos(2); % find the desired value of height
% clear h pos; % clear unnecessary variables
%
% I found this function useful because I needed to save my matlab figures
% in eps format. Since my vectors were big, I obtained very heavy images
% which could not be easily imported in my LaTeX documents. I have been
% looking around for a function that could help me to compress the images
% but I could not find it. So I wrote functions compress2D and compress2Dh
% which are very rudimental and probably need important implementation, but
% they did their job well enough for my purposes. I'm happy to share them:
% suggestions on how to implement them more efficiently are all welcome.
%
% % -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
%
% % EXAMPLE1: how to use the function
% N = 1000000; % A big number
% x = rand(N, 1); % random x values
% y = x .^ 2 + rand(N, 1) / 10; % y values along the square of x, plus small random noise
% h = plot(x, y, '.'); % heavy, uncompressed image ...
% saveas(h, 'heavyuncompressedimage.eps'); % ... heavy, uncompressed file
% h1 = dir('heavyuncompressedimage.eps'); % retrieve file size and display the result
%
% disp(sprintf('The file size, for the uncompressed image, is equal to %0.2f KB', h1.bytes / 1024));
%
% h = figure; % open a new figure (manually choose the desired size)
% set(h, 'Units', 'pixel'); % unit must be pixels
% pos = get(h, 'Position'); % get the Position
% pos(1:2) = []; % keep only the last two
% width = pos(1); % find the desired value of width
% height = pos(2); % find the desired value of height
% clear pos; % clear unnecessary variables
% alpha = 0.27; % a resolution coefficient
% ind = compress2D(x, y, width, height, alpha); % find indexes for compressed image
% h = plot(x(ind), y(ind), '.'); % compressed image ...
% saveas(h, 'compressedimage.eps'); % ... small file
% h2 = dir('compressedimage.eps'); % retrieve file size and display the result
%
% disp(sprintf('The file size, for the compressed image, is equal to %0.2f KB', h2.bytes / 1024));
% disp(sprintf('Congratulations! You have saved %0.2f%% of space', 100 - h2.bytes / h1.bytes * 100));
%
% % -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
%
% % EXAMPLE2: an extremely lucky case
% N = 1000000; % A big number
% x = randsample(50, N, 1); % random x values, integers from 1 to 50
% y = x .^ 2 + rand(N, 1); % y values along the square of x, plus small random noise
% h = plot(x, y, '.'); % heavy, uncompressed image ...
% saveas(h, 'heavyuncompressedimage.eps'); % ... heavy, uncompressed file
% h1 = dir('heavyuncompressedimage.eps'); % retrieve file size and display the result
%
% disp(sprintf('The file size, for the uncompressed image, is equal to %0.2f KB', h1.bytes / 1024));
%
% h = figure; % open a new figure (manually choose the desired size)
% set(h, 'Units', 'pixel'); % unit must be pixels
% pos = get(h, 'Position'); % get the Position
% pos(1:2) = []; % keep only the last two
% width = pos(1); % find the desired value of width
% height = pos(2); % find the desired value of height
% clear pos; % clear unnecessary variables
% alpha = 3000 / width; % a resolution coefficient
% ind = compress2D(x, y, width, height, alpha); % find indexes for compressed image
% plot(x(ind), y(ind), '.'); % compressed image ...
% saveas(h, 'compressedimage.eps'); % ... small file
% h2 = dir('compressedimage.eps'); % retrieve file size and display the result
%
% disp(sprintf('The file size, for the compressed image, is equal to %0.2f KB', h2.bytes / 1024));
% disp(sprintf('Congratulations! You have saved %0.2f%% of space', 100 - h2.bytes / h1.bytes * 100));
%
% % -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
%
% % EXAMPLE3: an unlucky case
% N = 1000000; % A big number
% x = rand(N, 1); % random x values
% y = rand(N, 1); % random y values
% h = plot(x, y, '.'); % heavy, uncompressed image ...
% saveas(h, 'heavyuncompressedimage.eps'); % ... heavy, uncompressed file
% h1 = dir('heavyuncompressedimage.eps'); % retrieve file size and display the result
%
% disp(sprintf('The file size, for the uncompressed image, is equal to %0.2f KB', h1.bytes / 1024));
%
% h = figure; % open a new figure (manually choose the desired size)
% set(h, 'Units', 'pixel'); % unit must be pixels
% pos = get(h, 'Position'); % get the Position
% pos(1:2) = []; % keep only the last two
% width = pos(1); % find the desired value of width
% height = pos(2); % find the desired value of height
% clear pos; % clear unnecessary variables
% alpha = 0.4; % a resolution coefficient
% ind = compress2D(x, y, width, height, alpha); % find indexes for compressed image
% plot(x(ind), y(ind), '.'); % compressed image ...
% saveas(h, 'compressedimage.eps'); % ... small file
% h2 = dir('compressedimage.eps'); % retrieve file size and display the result
%
% disp(sprintf('The file size, for the compressed image, is equal to %0.2f KB', h2.bytes / 1024));
% disp(sprintf('Congratulations! You have saved %0.2f%% of space', 100 - h2.bytes / h1.bytes * 100));
%
% % -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
%
% % EXAMPLE4: lucky (obvious) case
% N = 1000000; % A big number
% x1 = rand(N, 1); % random x1 values
% y1 = 0.5 * ones(N, 1); % horizontal line
% x2 = 0.5 * ones(N, 1); % vertical line
% y2 = rand(N, 1); % random y2 values
% h = plot(x1, y1, '.'); % heavy, uncompressed image ...
% hold on
% plot(x2, y2, '.r');
% saveas(h, 'heavyuncompressedimage.eps'); % ... heavy, uncompressed file
% h1 = dir('heavyuncompressedimage.eps'); % retrieve file size and display the result
%
% disp(sprintf('The file size, for the uncompressed image, is equal to %0.2f KB', h1.bytes / 1024));
%
% h = figure; % open a new figure (manually choose the desired size)
% set(h, 'Units', 'pixel'); % unit must be pixels
% pos = get(h, 'Position'); % get the Position
% pos(1:2) = []; % keep only the last two
% width = pos(1); % find the desired value of width
% height = pos(2); % find the desired value of height
% clear pos; % clear unnecessary variables
% alpha = 1.50; % a resolution coefficient
% ind1 = compress2D(x1, y1, width, height, alpha); % find indexes for compressed image
% ind2 = compress2D(x2, y2, width, height, alpha); % find indexes for compressed image
% h = plot(x1(ind1), y1(ind1), '.'); % compressed image ...
% hold on
% plot(x2(ind2), y2(ind2), '.r');
% saveas(h, 'compressedimage.eps'); % ... small file
% h2 = dir('compressedimage.eps'); % retrieve file size and display the result
%
% disp(sprintf('The file size, for the compressed image, is equal to %0.2f KB', h2.bytes / 1024));
% disp(sprintf('Congratulations! You have saved %0.2f%% of space', 100 - h2.bytes / h1.bytes * 100));
%
% % -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
%
% Author: Francesco Pozzi
% E-Mail: francesco.pozzi@anu.edu.au
% Date: 08-07-2008
%
% Check inputs
ctrl = isvector(x) & isnumeric(x) & isreal(x) & isvector(y) & isnumeric(y) & isreal(y);
if ctrl
ctrl = ctrl & (length(x) == length(y));
ctrl = ctrl & all(~isinf(x)) & all(~isinf(y));
end
if ctrl
x = x(:);
y = y(:);
else
error('x and y must be vectors of real numbers and of the same size!')
end
ctrl = isscalar(width) & isscalar(height) & isscalar(alpha);
if ctrl
ctrl = (width > 0) & (height > 0) & (alpha > 0) & (alpha * width <= 3000);
end
if ~ctrl
error('Specify a valid Position and a valid value for parameter alpha!')
end
% Range of x and y:
xmin = min(x);
xmax = max(x);
ymin = min(y);
ymax = max(y);
% resol1 and resol2 are intervals between minimum and maximum for x and y
% In each area of the grid thus obtained, only one point at most will be
% represented. All others will be discarded
if (xmin == xmax) & (ymin == ymax)
ind = 1;
warning('Hey! Your entire data lie on a single 2D point! Are you sure data are correct?');
return;
elseif ymin == ymax % if horizontal line ...
resol1 = ceil(alpha * width); % resolution along x-axis
resol2 = 1; % resolution along y-axis
j = floor((x - xmin) / (xmax - xmin) * (resol1 - 1)) + 1; % discretized x
k = ones(size(j));
elseif xmin == xmax % if vertical line ...
resol1 = ceil(alpha * width); % resolution along x-axis
resol2 = ceil(height / width * resol1); % resolution along y-axis
resol1 = 1; % resolution along x-axis
k = floor((y - ymin) / (ymax - ymin) * (resol2 - 1)) + 1; % discretized y
j = ones(size(k));
else
resol1 = ceil(alpha * width); % resolution along x-axis
resol2 = ceil(height / width * resol1); % resolution along y-axis
j = floor((x - xmin) / (xmax - xmin) * (resol1 - 1)) + 1; % discretized x
k = floor((y - ymin) / (ymax - ymin) * (resol2 - 1)) + 1; % discretized y
end
ind = zeros(length(x) / 2, 1); % final indices, hopefully much less than 50% of the total
m = 0;
grids = sparse(resol1, resol2);
for i = 1:length(x)
if grids(j(i), k(i)) == 0 % for each cell, keep only one point
grids(j(i), k(i)) = i;
m = m + 1;
ind(m) = i;
end
end
ind = ind(1:m); % keep only m values found before
