Add matrix formuation section

project.md

@@ -267,3 +267,54 @@ sample to the left of the desired location, and for integration I just use a lef
 with these lazy techniques, the reconstruction looks quite good.
 
 ![brain](../assets/img/project_fbp_reconstruction.png)
+
+
+## Total Variation / Compressed Sensing
+
+Unfortunately this is still a work in progress. I started implementing a [modern
+approach](../assets/pdf/2011_Hansen_Jorgensen.pdf), but found I didn't have enough background to
+make it work. They gloss over some details that I tried to get around with brute force, only to find
+that the resulting problem was computationally intractable. So yesterday I started over,
+re-implementing results from one of the original compressed sensing
+[papers](../assets/pdf/2004_Candes_Romberg_Tao.pdf). However I did learn some things from my failed
+attempts...
+
+
+## Matrix Formulation
+
+Both techniques I've fully implemented so far involve only three types of operations: DFTs,
+interpolation, and multiplication (for the $$\vert \omega \vert$$ filter). All these operations are
+linear, so we can express them as matrices and reformulate each technique as a single matrix
+multiplication with a vector. In practice this fact alone is useless, since the resulting matrix
+ends up being enormous. (In fact if your image and sinogram are n by n pixels, the matrix will have
+n^4 entries.) But this formulation is used to derive the theory of many total variation versions.
+
+In the process of my first failed total variation implementation, I ended up with most of the parts
+of the Fourier and filtered back projection algorithms implemented as matrices. So let's use them
+and create a sinogram with a single matrix multiplication. As mentioned the matrices are huge, so
+here I'll work with a 32 by 32 brain image.
+
+![dft](../assets/img/project_tiny_brain.png)
+
+The DFTs were easy. Very similar to the cosine transform we used in the last problem set.
+
+![dft](../assets/img/project_matrix_dft_re.png)
+![dft](../assets/img/project_matrix_dft_im.png)
+
+To construct the polar interpolation matrix, I rewrote my interpolation routine to drop the weights
+it calculated in the appropriate entries in a matrix rather than summing things up as it goes. These
+weights depend on the type of interpolation used and the sizes of the images involved.
+
+![dft](../assets/img/project_matrix_polar_re.png)
+![dft](../assets/img/project_matrix_polar_im.png)
+
+Finally it's just some (inverse) DFTs to get the sinogram. They can all be expressed simultaneously
+in one matrix.
+
+![dft](../assets/img/project_matrix_sinogram.png)
+
+I bump up the pixel count for the polar projections so that the resampling doesn't lose as much
+information. Thus the final matrix that performs all three of these operations at once has
+$$2 \cdot 64^2$$ rows and $$32^2$$ columns, for a total of 8,388,608 entries. (I could chop this in
+half if I didn't bother computing the imaginary part of the sinogram. It should be zero but it's a
+nice sanity check.)