Image Compression Project

The goal of the project is to compress images using a version of the least-squares method. This notebook walks you through the process of decomposing an image into blocks, and then into vectors, and reconstituting the image from a set of vectors. In the middle, we discuss the process of finding compression matrices and how to make Mathematica convert our image vectors into shorter, compressed vectors by least-squares, and then decompressing them by multiplying by the compression matrix. We give a discussion of one “automatic” way to find the A matrix from a training set of images called Principal Component Analysis (PCA). And finally, we give specific guidance for the project itself.

You’re expected to harvest code from this notebook as desired to work on your project and to test your various candidate compression matrices. 

--Jason Cantarella, Spring 2014.

Helper Functions.

These are some functions for converting images to blocks, then to vectors, and back again. We start with a test image:

DogImage =ColorConvert[ , "Grayscale"]

ImageDimensions[DogImage]

We now use ImagePartition and ImageData, as well as some reshaping functions, to turn this image into 32 x 32 blocks, and then into 1024-vectors. Notice that we’ll cut off some of the right and bottom of the image to reduce things to the nearest multiple of 32 pixels. That’s ok.

ImagePartition[DogImage,32] // Grid

ImageToRawVectors[Img_] :=Flatten[ Map[ImageData,ImagePartition[Img,32],{2}],{{1},{2},{3,4}}];

Applying this to the dog image, we get

DogRawVectors = ImageToRawVectors[DogImage];

If we take dimensions, we can see that this is a 7x7 grid of blocks (49 vectors total).

Dimensions[DogRawVectors]

In order to make the vectors easier to work with, we’re going to zero-mean them (and store the means), and rescale each one to full contrast (and store the contrast multipliers). This will give you better results from the linear algebra. The resulting vector will have average entry zero and each coordinate contained in the interval [-1,1], where -1 represents pure black and +1 represents pure white.

ZeroMean[v_] := {Mean[v],v- Mean[v]};
ContrastAdjust[{μ_,v_}] := If[Max[Abs[v]] > 0.01,{μ,Max[Abs[v]],(1/Max[Abs[v]]) v},{μ,1,v}];
Preprocess[v_] := ContrastAdjust[ZeroMean[v]];

We’re now capable of writing the function which takes images and returns lists of vectors (together with the extra information we need to reconstruct the image: the mean and the scale factor)

ImageToData[Img_] :={#[[;;,;;,1;;2]],#[[;;,;;,3]]}&[Map[Preprocess,ImageToRawVectors[Img],{2}]];

The output of ImageToData is an array of means and scale factors (that you don’t have to worry about) and an array of 1024-vectors which you will work with. We use the functions like this:

{DogAuxData,DogVectors} = ImageToData[DogImage];

DogAuxData

The DogVectors themselves are too long to make a really meaningful display. But we can check the dimensions and make sure that this is really a 7x7 array of vectors, each of length 1024:

Dimensions[DogVectors]

NB: Why can’t we flatten this into a list of 49 vectors? Because we don’t know that the image is square. So we’ll have no way of knowing how many rows and columns the grid of blocks should have when we’re rebuilding the image unless we keep each vector nested safely in the proper place in an array.

As they say when you’re fixing your car, “installation is the reverse of removal”. So now we’re going to give code to undo this process and translate a list of vectors (together with corresponding AuxData) back into an image. 
The first step is to postprocess the vectors to return the contrast and mean to the previous values. Again, the major problem here is just reuniting the auxiliary data with the correct vectors and running functions on in without disturbing the array structure of the data.

ReMean[{μ_,v_}] := (# + μ)& /@ v;
DeContrast[{μ_,λ_,v_}] := {μ,λ v};
Postprocess[{AuxData_,Vectors_}] := Map[ReMean,Map[DeContrast,Map[{#[[1]],#[[2]],#[[3;;]]}&,Join[AuxData,Vectors,3],{2}],{2}],{2}];

Dimensions[Postprocess[{DogAuxData,DogVectors}]]

Now we need to assemble each of these vectors back into a 32 x 32 array, then to a grayscale image, then assemble these images into a final picture.

DataToImage[{AuxData_,Vectors_}] := ImageAssemble[Map[Image,Map[Partition[#,32]&,Postprocess[{AuxData,Vectors}],{2}],{2}]];

DataToImage[{DogAuxData,DogVectors}]

Just to check that the image survived all of these manipulations intact (well, except for clipping), let’s look at the original image again:

DogImage

Here’s another test image:

SecondDog =

{SecondDogAux,SecondDogVectors} = ImageToData[SecondDog];

You can see that this is a much higher resolution image (and it’s not square):

Dimensions[SecondDogVectors]

DataToImage[{SecondDogAux,SecondDogVectors}]

The Compression Process

We can now convert an image into a m x n array of 1024-vectors (Subscript[b, ij]). In the compression process, we’ll convert each (Subscript[b, ij]) vector into a corresponding (and shorter!) (Subscript[x, ij]) vector by solving the least squares problem

(min_(Subscript[x, ij])) |(Subscript[Ax, ij]) -(Subscript[b, ij])|

given an “A” matrix. Since the A matrix is tall-skinny, the (Subscript[x, ij]) vectors are shorter than the (Subscript[b, ij]) vectors (and this is where the compression comes in).

Of course, since we’re going to solve a whole bunch of least-squares problems in a row, it makes sense to compute the Moore-Penrose pseudoinverse of our A matrix first, then use it to solve all the least squares problems. 

This seems reasonable enough, but it doesn’t make much sense at all without an explicit demonstration. The idea is that the column space of the A matrix should almost contain “every” 1024 vector which might occur as a 32 x 32 block in “any” image. But where might we get such an “A” matrix? Are we really saying that there’s some deep structure in 32 x 32 matrices which come from image data which distinguish them from arbitrary 32 x 32 matrices? 

Well, yes. 

In fact, completely random noise is actually quite unlikely in a 32x32 block of an actual image. This is what it looks like:

Image[Clip[Abs[RandomVariate[NormalDistribution[],{32,32}]]],Magnification->4]

Here are some “primitive” blocks which look a little more like blocks of actual images:

LinearGradientBlock[θ_] := Module[{raw},
raw =Table[1.0i Cos[θ] +1.0 (32-j) Sin[θ],{j,1,32},{i,1,32}];
raw -= Min[raw];
raw /= Max[Abs[raw]];
Image[raw,Magnification->2]
];

Table[LinearGradientBlock[θ],{θ,0,2π,2π/10}]

We can build an A matrix with small modifications of these 10 columns. (It will be pretty terrible, but this is just an illustration.)

LinearGradientVector[θ_] := Module[{raw},
raw =Table[1.0i Cos[θ] +1.0 (32-j) Sin[θ],{j,1,32},{i,1,32}];
raw -= Mean[raw];
raw /= Max[Abs[raw]];
Flatten[raw]
];

AMatrix = Transpose[Table[LinearGradientVector[θ],{θ,0,2π-0.01,2π/10}]];

Note that we don’t want to go all the way to 2π because we’ll repeat the first vector, which will lead to a very ill-conditioned matrix.

Dimensions[AMatrix]

Now we can write the compression code:

CompressImage[A_,Vectors_] := Module[{APlus},
APlus = PseudoInverse[A];
Map[Flatten[APlus.Transpose[{#}]]&,Vectors,{2}]
];

We can see that the resulting grid of vectors contains 10-vectors now, instead of 1024-vectors.

Dimensions[CompressImage[AMatrix,DogVectors]]

CompressedDogVectors = CompressImage[AMatrix,DogVectors];

The Decompression Process

The decompression process is easy. We just compute the 1024-vectors (Subscript[Ax, ij]) by matrix multiplication.

DecompressImage[A_,CompressedVectors_] := Map[Flatten[A.Transpose[{#}]]&,CompressedVectors,{2}];

DecompressedDogVectors = DecompressImage[AMatrix,CompressedDogVectors];

Dimensions[DecompressedDogVectors]

We can view the (relatively unpleasant) results using DataToImage:

DataToImage[{DogAuxData,DecompressedDogVectors}]

You can see that this preserves some features of the original image, though not very many. In fact, the smoothness of the gradients is a real problem here-- we’ve lost all of the edges in the image. You can imagine various strategies to improve on this choice of basis blocks: one approach that worked well in the past was to carefully come up with various colums which included some hard edges, some linear gradients, and some quadratic functions.

Coming up with better A matrices: Machine Learning through PCA

Another approach to coming up with better A matrices is to use a training set of images. To do this, we’ll have to learn a little bit about principal component analysis, or PCA. There’s a lot of good stuff to read on PCA, since it’s a foundational technique in data analysis. I recommend the wikipedia page as a reasonable place to start, but we’ll take a look at some of the basic ideas here.

First, suppose we have a data set which almost lies in a subspace of higher dimensional space. Since the details of how I’m generating this data don’t matter to the discussion, this code is undocumented.

CompletelyFakeData[n_] := Module[{A},
A = {{1,2},{0,4},{2,3}};
Table[Transpose[A.Transpose[{
RandomVariate[NormalDistribution[]] {4Cos[#],Sin[#]}& [RandomVariate[UniformDistribution[{0,2π}]]]}] + 0.1*RandomVariate[NormalDistribution[],{3,1}]][[1]],{n}]
];

Rotating the picture below around should convince you that this data is pretty close to lying in a plane inside 3-space.

ListPointPlot3D[CompletelyFakeData[500]]

Suppose we build a matrix out of lots of our fake data, with a lot of samples (we’ll try 5,000). The question is how we can recover the plane from all this data.

FakeDataSet = CompletelyFakeData[5000];

Dimensions[FakeDataSet]

The first trick is that we can build a cross-correlation matrix for our data by dotting the columns of the dataset together (this is also known as building the Gram matrix for our set of 3 vectors in 5000-space).

CorrelationMatrix = Transpose[FakeDataSet].FakeDataSet;
CorrelationMatrix // MatrixForm

If we take the singular value decomposition of this matrix, we get a special basis which should reveal the planarity of our data:

{U,Σ,V} = SingularValueDecomposition[CorrelationMatrix];
U // MatrixForm
Σ // MatrixForm
V // MatrixForm

We see that the U and V matrices are the same (this isn’t surprising, because the correlation matrix is symmetric by construction), and that the singular values have really different sizes: 1 very large, 1 medium, and 1 tiny. These are telling us the importance of the columns of U in representing our data set. 

We can make an A matrix for this data directly from the U columns.

PCAMatrix = U[[1;;3,1;;2]]

To see the sense in which this “represents” the data, we’ll try writing the fake data points as linear combinations of these column vectors. First, we take the PseudoInverse of the PCAMatrix:

PCAMatrixPlus = PseudoInverse[PCAMatrix]

Then we use it to solve the various least-squares problems, keeping in mind that we have to convert the row vectors in FakeDataSet to column vectors in order to do the matrix multiplication.

LeastSquaresRep = PCAMatrixPlus.# & /@ (Transpose[{#}]&/@ FakeDataSet);

Dimensions[LeastSquaresRep]

To approximately recover the data, we multiply the least squares representation by the PCAMatrix. Of course, the results are column vectors again, so we have to reshape things in order to get lists which we can plot.

RecoveredData = Flatten[PCAMatrix.# & /@ LeastSquaresRep,{3}][[1]];

Plotting the recovered data together with the original data, we can see that the recovered data set is a remarkably good approximation to the original data, even though the recovered data is exactly 2-dimensional-- that is, we have lost almost no information by converting the data into a 2-dimensional representation, even though the original data was all 3-dimensional! This process is generally called dimension reduction and it’s a really important idea in data analysis.

ListPointPlot3D[{FakeDataSet,RecoveredData}]

We can also look at some points to see that the approximation in the recovered data is pretty good.

RecoveredData[[1;;5]]

FakeDataSet[[1;;5]]

Of course, in this example we already knew that the data was two-dimensional, so we knew how much of the U matrix to take. If you didn’t already know the dimension of your data set, you’d have to discover it by looking at the singular values and seeing how they transition from large to small. If the singular values of the correlation matrix are all comparably large, you’ve learned that your data is not approximately confined to a linear subspace of the original space. It may still be lower-dimensional (for instance, it could lie on a curved surface in 3-space) but that subject is called nonlinear dimension reduction and it’s quite a bit harder than the PCA technique we’re demonstrating here.

Many students used the PCA method to build their compression matrices. Let’s walk through this process a little bit. Suppose that you take 1000 test images, and use them to build 100,000 sample vectors, each of size 1024. The data matrix from your images is then 100,000 x 1024. For instance, you’d have:

FakeDataMatrix = RandomVariate[NormalDistribution[],{100000,1024}];

You are then taking the SVD of the correlation matrix, which is

CorrelationMatrix = Transpose[FakeDataMatrix].FakeDataMatrix;

On my (2014) laptop, this operation takes 10-20 seconds. Taking the singular value decomposition of this should be quick, because the correlation matrix is always 1024 x 1024, regardless of the size of the training set.

{U,Σ,V} = SingularValueDecomposition[CorrelationMatrix];

Of course, this data has no linear structure to find, so if we plot the singular values, we see basically a probability distribution.

ListPlot[Diagonal[Σ],AxesOrigin->{0,0}]

In fact, these singular values are all the same order of magnitude.

Max[Diagonal[Σ]]

Min[Diagonal[Σ]]

Other methods: The sky is the limit!

You’re encouraged to read, to be creative, and to try out seeming crazy ideas for improvements. Here are some starters:

1) What if you used an optimization procedure like simulated annealing to try to find the best A matrix numerically? You’d have to assemble a training set of blocks and see how well each one is recovered by a given A matrix, then combine these into an aggregate score function for the candidate A matrix. Then you’d use NMinimize (or similar) to try to vary the A matrix to find the best scoring one.

2) You could select the columns of your A matrix directly from existing images. 

3) You could use wavelets or Fourier functions to generate the columns of the A matrix.

4) If you want to start with Ryan Livingston’s (winning) set of A matrices from last class (and try to improve them), they are on the webpage. You’ll have to figure out how to load them yourself. :)

The Final Project Assignment, in detail.

In order to preserve my sanity when assembling the results of the final project, I’d like to be pretty specific about how the results should be written up and turned in. 

All projects should be turned in by email. 

The email should have the subject line “Math 4510 Final Project”, in order to make sure that I don’t miss it. 

Use Export to generate three files called 

	YourLastName50.dat : a 1024 x 512 matrix 
	YourLastName90.dat : a 1024 x 102 matrix
	YourLastName99.dat : a 1024 x 10 matrix
	
from your matrices. Attach these files to your email. 

For example, to export our example 1024x10 AMatrix, we type:

Export[NotebookDirectory[]<>"Cantarella99.dat",AMatrix]

You can load such data with the import command, such as

ImportedA = Import[NotebookDirectory[] <> "Cantarella99.dat"];

We can check that things saved and loaded correctly with:

Dimensions[ImportedA]

Dimensions[AMatrix]

ImportedA[[1;;4,1;;4]] // MatrixForm

AMatrix[[1;;4,1;;4]] // MatrixForm

Next, turn in a Mathematica notebook or PDF file called 

	YourLastNameFinalProject.pdf (or .nb)
	
containing a detailed writeup of the various experiments you performed in order to come up with your matrices. I’m particularly interested in false starts and good ideas that didn’t work out, as well as a sense of your thought process as you worked your way through the assignment, so this is more like a lab notebook than a presentation of the best method you eventually came up with. You are welcome to attach additional files (for instance, your training images or code written in another computer language).

The majority of your grade will be based on this writeup. Of course, it’s great for you to get an excellent set of compression matrices in the end, and the results of the class competition will be determined by the quality of your compression matrices. But even poor matrices can result in a good grade if you tried interesting experiments and wrote them up clearly-- in this case, you should try to determine why your method didn’t work, and form a hypothesis which could help someone else come up with better results in the future.  

Rules of the Road.

The final project is an open-book, open-internet, open-discussion research problem. You are encouraged to work in groups. However, the usual rules of scientific honesty apply.

1) Cite your sources. If you got a great idea from somebody’s webpage, or wikipedia, or reading a paper, that’s terrific! You should cite that source in your writeup. Even if you’re not directly using somebody else’s method, but were inspired by their approach, you need to write “I got this idea from the Wikipedia page on X” or “This is from a really great blog post <link>”.

2) Include your collaborators. You should acknowledge everybody who you had helpful conversations with during the process of working on your project. It’s ok to get data from each other, too, but again you need to cite it “I used the training set of images assembled by Ryan”.  

3) Write your own code, and the text of your own writeup. It is ok to ask your friends to help you look over your code to debug it. It is not ok to get them to write it for you or to copy-and-paste their code and then modify it for your purposes. Similarly, it’s ok to discuss the methods and the writeups, but you should not copy text from each other.

Please let me know if you have further questions, or trouble running this code.