%% Image Compression and Manipulation

%%

clear all
close all
format short
set(0, 'DefaultAxesFontSize', 16)

%%
% In Chapter 5, we discussed using MATLAB to read, write, display, and
% manipulate image data.

%%
% In the first part of this section, we illustrate a general
% technique for data compression based on a tool from linear algebra called
% ``singular value decomposition.''  A variety of names are associated with
% this technique, including ``principal component analysis,'' ``proper
% orthogonal decomposition,'' and ``empirical orthogonal functions.''
% In scientific applications, the purpose of the technique is often to
% identify significant patterns in large data sets and focus computational
% or visualization techniques on these patterns.  The compression example
% below is meant to illustrate visually how most of the perceptible
% information in a data set can sometimes be contained in a much smaller
% data set.

%%
% At the end of this section, we give some examples of commands, mostly
% from the Image Processing Toolbox, to adjust the appearance of an image.

%% Singular Value Decomposition and Compression
% The command *|svd|* computes the singular value decomposition of a matrix;
% here we will use the ``economy'' form of the output (see the online help
% for other options).  While the discussion below is for a general matrix
% |A|, we will illustrate with a specific example in MATLAB:

A = [1 2 3 4; 5 6 7 8; 9 1 2 3];
[U, S, V] = svd(A, 'econ')

%%
% To explain what the singular value decomposition is, let $n$ be the
% smaller of the number of rows in |A| and the number of columns in |A|.
% Then the three output matrices have the following properties:
%
% * |A = U*S*V'|;
%
% * |S| is an $n \times n$ diagonal matrix; its diagonal entries are, in
% decreasing order, the ``singular values'' of |A|;
%
% * |U'*U = V'*V = eye(n)| (for each matrix, its columns are an orthonormal
% set of vectors).
%
% As a consequence of the first two properties, we can decompose A as a
% sum |A = U(:,1)*S(1,1)*V(:,1)' + U(:,2)*S(2,2)*V(:,2)' + $\cdots$ +
% U(:,n)*S(1,n)*V(:,n)'|.  (Recall that in MATLAB, |U(:,j)| is column
% $j$ of |U|.)  Let's check this identity for our example:

B = zeros(size(A));
for j = 1:3
    B = B + U(:,j)*S(j,j)*V(:,j)';
end
B

%%
% If |A| is $n \times k$ where $k \geq n$, then each term in the sum is an
% $n \times k$ matrix, but it can be computed from only $n + k$ numbers:
% the $n$ entries in |U(:,j)*S(j,j)| and the $k$ entries in |V(:,j)|.
%
% One can show that the singular values of |A| are the square roots of the
% eigenvalues of |A*A'|.
% Often the singular values of |A| vary greatly in size, especially if
% |A| is large.  Since the columns of |U| and |V| are unit vectors, the
% size of |U(:,j)*S(j,j)*V(:,j)'| is determined by |S(j,j)|, and thus
% the summation above is dominated by the terms with the largest singular
% values.  If we approximate the sum by its $r$ largest terms, we can
% store the information in the sum using $(n + k)r$ numbers.  This is
% much smaller than the $n*k$ numbers in |A|, provided $r$ is much smaller
% than |n|.
%
% Since each term in the sum is a matrix with rank $1$, the sum of the
% first $r$ terms has rank (at most) $r$.  This sum is in fact the best
% rank $r$ approximation to |A| in an appropriate least-squares sense.
%
% Let's look at the rank $1$ approximation to our $3 \times 4$ example:

U(:,1)*S(1,1)*V(:,1)'

%%
% The result doesn't look much like |A|.  Let's try the rank $2$
% approximation, using a shortcut formula:

U(:,1:2)*S(1:2,1:2)*V(:,1:2)'

%%
% This is a pretty good approximation, which follows from the fact
% that the third singular value of our |A| is so much smaller than the
% other two.


%% Compression Example
% We now load a sample photograph into MATLAB with *|imread|* and check its
% size (we will display the photograph later, after converting the data to
% shades of gray): 

photoarray = imread('sample.jpg');
size(photoarray)

%%
% The values in the array are red, green, and blue intensities (integer
% values from 0 to 255) for each of 1536*2048 pixels.
%
% Next, we separate the colors and convert to floating-point values between
% 0 and 1.  We use single precision because it is sufficient for our
% calculations, and our objective is to store efficiently the essential
% data in the image.  (We are still converting each byte of the integer
% image data into four bytes of floating-point data, but the floating-point
% is necessary for the calculations we do below.)
redpix = single(photoarray(:,:,1))/255;
greenpix = single(photoarray(:,:,2))/255;
bluepix = single(photoarray(:,:,3))/255;

%%
% To illustrate the ideas below more simply, we convert to a
% black-and-white (more precisely, ``grayscale'') image.
% The following is a standard formula to convert RGB intensities to
% gray intensities:
graypix = 0.299*redpix + 0.587*greenpix + 0.114*bluepix;

%%
% Since we have converted to floating-point values, we use *|imagesc|* rather
% than *|image|* to display the photograph.  It shows part of Owachomo
% Bridge in Utah, and was taken by one of the authors.  We use |colormap|
% to get a grayscale display, and *|axis|* to turn off axis labels and
% ensure the correct aspect ratio:
imagesc(graypix)
colormap(gray), axis off image

%%
% Now we compute the singular value decomposition of the grayscale data.
[U,S,V] = svd(graypix);

%%
% To see how fast they decay, we graph the first 100 singular values.  We
% open a new figure window so that the photograph remains visible if this
% M-file is used interactively:
figure
plot(diag(S(1:100,1:100)))

%%
% The following function computes and displays the rank $n$ approximation
% to our image:
showrank = @(n) imagesc(U(:,1:n)*S(1:n,1:n)*V(:,1:n)');

%%
% Below we show the approximations of rank 1, 3, 10, and 30:
colormap(gray)
rankvals = [1 3 10 30];
set(0,'DefaultAxesFontSize', 10)
for m=1:4
    subplot(2,2,m)
    showrank(rankvals(m)), axis off image
    title(['Rank ' num2str(rankvals(m))])
end
set(0,'DefaultAxesFontSize', 16)

%%
% A rank 1 matrix can only produce a ``plaid'' image, but by the time we get
% to rank 10, the image starts to look like a photograph, and by rank 30 it
% is clear we are looking at a landscape.  Next, let's compare the rank 100
% approximation to the full-resolution image.  We clear the existing figure
% so that the new image is not shown in a subplot:
clf
showrank(100); axis off image

%%
% The rank 100 approximation is not as sharp as the original, but the loss
% of quality might not be noticed if the original were not available.
% As a rough comparison to JPEG compression, let's store the grayscale
% image as a JPEG file with default parameters and see how big it is:

imwrite(graypix, 'graysamp.jpg', 'jpg')
d = dir('graysamp.jpg'); d.bytes

%%
% The raw image data consist of $1536 \cdot 2048$ intensities that could
% be stored as 1 byte apiece, for a total of 3 megabytes.  Since our rank
% $r$ approximation requires $(1536 + 2048)r = 3584r$
% floating point numbers, and each of our numbers is stored in single
% precision (4 bytes), we can compute as follows the rank for which our
% approximation would require about the same amount of data as the JPEG
% file:

r = d.bytes/3584/4

%%
% Of course, the JPEG file uses a much more sophisticated algorithm and
% has much better quality, but we are at least
% getting a recognizable image with comparable compression, and with some
% programming effort we could reduce the number of bytes used per number.

%% Image Manipulation and the Image Processing Toolbox
% In Chapter 5, we showed how to create a mirror image and to crop an image
% by performing these operations on the array of image data.  You can also
% perform various mathematical operations on the data to enhance or
% suppress certain features of the image.  For example, for our grayscale
% image |graypix|, with values between |0| and |1|, the square-root
% function brightens the image by making all the intensities closer to |1|:
imagesc(sqrt(graypix)), axis off image

%%
% The Image Processing Toolbox includes a variety of algorithms that
% manipulate image data.  If you have this toolbox installed, it will
% appear on the menu when you open a new Help Browser tab (e.g., by typing
% *|doc|*); click to get an overview of available commands.  Some commands
% provide shortcuts for operations you could perform with basic MATLAB
% commands; for example, the command *|rgb2gray|* will produce almost the
% same result as the arithmetic commands we used above to convert a color
% image to grayscale.  Other commands perform more sophisticated
% algorithms, and we give a few examples here.

%%
% One of the
% contrast adjustment commands, *|adapthisteq|*, will both brighten dark
% areas of the image and darken light areas; it increases contrast in each
% part of the image.
%imagesc(adapthisteq(graypix)), axis off image

%%
% Like the brightened image before, visibility is improved the in shadowy
% part of the photograph, but unlike the brightened image, the clouds have
% become darker, and the rock striations have become more pronounced.

%%
% The toolbox also has a variety of filters that can ``soften'' an image;
% the simplest of these, *|imfilter|*, will convolve the image with a matrix
% that you specify.  Convolving tends to blur edges in the image, which is
% sometimes the desired effect.  On the other hand, the command *|wiener2|*
% (named for Norbert Wiener) reduces graininess in surfaces while trying to
% maintain sharp edges.
%imagesc(wiener2(graypix, [9 9])), axis off image

%%
% Here we have chosen to filter over a $9 \times 9$ grid rather than the
% default $3 \times 3$; as always, see the online help for other options
% to this and the other commands we have used.  The result de-emphasizes
% striations in the rocks, but also blurs the trees and smaller clouds.

%%
% Other toolbox commands include algorithms to deblur images; starting with
% MATLAB 2013a, the newest of these is *|imsharpen|*.