Lab assignment L3

Applications of Eigen values and Vectors

Objective

Procedure

We will interpolate between images with different lighting conditions/object positions using a function of the form I(x) = B*Y(x), where B is a linear basis for our image space and Y(x) can be acquired by an interpolation function. The steps of the lab involves: (1) computing the subspace basis B by finding the eigenvalues and vectors (Ch 4) to a sample set of images, and (2) performing data fitting over the computed set of data (Heath Ch 3,7).

Exercise 1 demonstrates the effect of applying PCA over synthetic data. Exercise 2 and 3 (to be demo in lab session) are examples that will help the understanding of the image modulation. Exercise 4 and 5 are related to the real application.

Exercise 1: Eigenvalues and Eigenvectors

To generate a random, eclipse-like 2D Guassian distribution and visualize it, the below sample code is provided:

      % Gaussian distr, Box-Mueller method
      n = 400
      U = rand(2,n);
      v = 2*pi*U(2,:);
      G = ones(2,1)*sqrt(-2.*log(U(1,:))).*[cos(v);sin(v)];
      
      plot(G(1,:),G(2,:),'.')
      axis equal
      
      % Transform
      A = [1 .5
           .5 0.5];
      t = [8; 1];
      Y = A*G+t*ones(1,n);
      plot(Y(1,:),Y(2,:),'.')
      axis equal
      hold on
      

t = 0:0.1:2*pi;
x = [cos(t); sin(t)];
plot(x(1,:),x(2,:))
axis equal

A = [1 .2
     .2 1];
y = A*x
plot(y(1,:),y(2,:))
axis equal

[X,D] = eig(A)
hold on
axis equal
plot(D(1,1)*[0 X(1,1)],D(1,1)*[0 X(2,1)],'g')
plot(D(2,2)*[0 X(1,2)],D(2,2)*[0 X(2,2)],'r')
hold off

To think about: In the above we only rotated and scaled. More generally we can also add translation y = Ax+t. You might try on your own to make a routine which can plot any ellipse given its parametric form (as in Heath Computer problem 3.5).

A more common problem in practice is that we have data that we want to somehow characterize. The data may have significant variation along some axis and less along others. A common model is that of a Gaussian distribution. If there is little variation along a particular axis, that dimension could be omitted altogether. This is the idea will be further exploited with the images below.

(Optional, TBA in Lab) A common use for PCA is to segment data into in and out sets. For instance, in images a region that looks like a solid color to the eye usually has some variation. To quantify this variation a sample of the image can be sampled. Once a gaussian model is fit, the rest of the image can be analysed as to what fits the model or not. "Bluescreening" in TV production, and fill tools in graphics/image editors can be implemented like this.
Implementation sketch:

• Select a (rectangular) region of a solid color.
• Perform PCA on the 3D RGB color space using the colors of the pixels as data.
• Use recovered PCA transform (A, t) to examine all pixels in the image. Mark pixels that are matching the Gaussian model. This step can be performed as a coordinate transform where the regular RGB color space is transformed using the A and t, s.t. the resulting coordinates are aligned and scaled with the PCA axes. Now selecting the pixels simply amounts to checking the norm of all pixels. In pixels have norm less than k (e.g. norm less than 1).

Exercise 2: A moving wave image

The effect of image motion can be obtained by modulating a sin/cos basis. Considering the following formula:

sin(ax+bt) = sin(ax)cos(bt)+cos(ax)sin(bt)

where a and b are constants, x is the spatial variable and t is time. We can interpret the above formula as the modulation of a sin and cos spatial basis with time varying coefficients - cos(bt), sin(bt). The example can be extended in 2D. Download the sincos_bs.mat data. B contains the sin/cos image basis (as columns) and Y is an example set of temporal modulation coefficients (each set of coefficients is a column, hence Y(:,1) are the coefficients for time 1, Y(:,2) for time 2 etc.). Show the modulated images on the screen in quick succession resulting from the different coefficients in Y. You can use imshow, then drawnow to force immediate rendering to the screen.
You can use the function renderim.m Ii = renderim(Y(:,i),B,imsize); imshow(Ii,[]); drawnow; pause(0.1);

Construct a similar basis that has twice the frequency of the given one and show the effect of motion through modulation.

What is the effect of changing b ?

Exercise 3: Small image motion

I(u+du,v+dv) = I(u,v) + [Ix(u,v) Iy(u,v)][du dv]' + hot

GI = [Ix Iy] represents the image gradient. Again we can interpret the above formula using image modulation. In this case the basis B is formed by [I, Ix, Iy] and the coefficients are of the form [1 du dv]. The effect of the modulation is a small image shift. In the below we compute the derivative by taking differences in spatial coordinates. It is also possible to use a time series of images, akin this example of a bus.

Repeat the procedure for another image of your choice. To calculate the image gradient (derivatives) use the Matlab function 'gradient'. Study the effect of changing the magnitudes of the coefficients in Y. What happens if you pick larger Y values (e.g. 5:10)? Why? What type of image does this procedure work best for? Explain appealing to the Taylor expansion what happens when it doesn't work.

Exercise 4: Combining the linear image subspace B with interpolation

Note: An interactive script of this exercise can be found Here (lab3_ex4.mlx). Optionally, you can choose to start with your own script.

Two precomputed data with basis B are provided, modulation coefficients Y and correxponding image index vector X.
light_pca.mat (for original images in directory: images/lighting/1)
obj_pca.mat (for images/objects/13)

Note: To put your data from Exec 3 below in the same format as the downloadable data, you need to add the meanim to B and use the transpose of the Y matrix from matlab "eigs" command, as noted under point 4, Exc 3 (It = [meanIm;B]*[1;Y(t,:)']). To understand this exercise you need to reflect on Exc 4 and 5 and notice how B holds a basis that spans image variation here. It is just of larger dimension than the two previous exercises. The columns of Y hold coefficients of the basis. You can re-generate the input images through:

load light_pca

for c = 1:size(Y,2)
  Ic = renderim(Y(:,c),B,imsize);
  imshow(Ic)
  drawnow
end

1. First of all, try to understand what each of the matrices inside the data file stands for. Additionally, use the render code provided to visualize the reconsturcted images from basis B and Y.


2. Now, densely generate new Y_new values for each of the datasets with a piecewise interpolation. The process is shown as follows:
   • Construct a vector X_new more densely sampled than X. That is, densely sample X with 4 times as many values.
   
   
   • Given X, Y, and your densely sampled X_new, perform piecewise interpolation and generate the new coefficients Y_new. Y_new will contain interpolated values of Y.
     You can use MATLAB's built-in piecewise interpolation method (interp1), and returns the interpolated values Y_new corresponding each element in the vector X_new.
   
   
   • Render the image sequence with your newly generated coefficients X_new and Y_new.
     By generating new coefficients X_new and Y_new using piecewise interpolation, you can multiply them through the basis, add the mean image (if using your own generated basis), and reshape them to the size of the original images.
   
   


3. Repeat the previous exercises with polynomial interpolation of various degrees. You can use your your polynomial interpolation function from lab 2: Y_new=polynomial_interp(X,Y,X_new,n) which fits a polynomial of degree n to Y, and returns the values of that polynomial for each element in the vector X_new. Additionally, you can also use MATLAB's built-in functions "polyfit" and "polyval".


4. Verify that your method correctly interpolates/approximates intermediate images by leaving out some components in Y. Compare how well Y_new represents Y for piecewise interpolation if you omit some components in Y.
   First try getting the first few components of Y interpolated right, (ie focus on Y(1:6,:), and ignore Y(7:m,:), (or set Y(7:m,:)=0).


5. Verify that your method correctly interpolates/approximates intermediate images by leaving out some of the original data. From the original Y, pick one (or more) t-values, and omit the sets of coefficients Y(:,t) corresponding to image t. Compute the interpolated coefficients Y_new and be sure to now include the t's corresponding to the ommited values. Does the interpolation correctly recover the Y? (ie how close is Y_new(:,t) to the dropped Y(:,t)?) Render the image at t and compare with the original. How close do they look?

Exercise 5: Finding a linear subspace for a set of images

Note: An interactive script of this exercise can be found Here (lab3_ex5.mlx). Optionally, you can choose to start with your own script.

Several sets of images are provided at images.

Note: Students can directly download the linear subspace basis for two data sets to test their results for Exercise 3.
light_pca.mat (for original images in directory: images/lighting/1)
obj_pca.mat (for images/objects/13)

They proceed in the following way:

1. Load in a sample image set (light or position) with m images, and flatten them so the pixels in each image is a column vector, for instance using ImV=Im(:);. The image as a vector will be a quite long vector, e.g. q = 10,000 for a 100 x 100 image. Then put all the m column vectors into an q x m matrix M. The images are located in images, and a suitable image loading function is provided below.
2. Perform a Normalized Eigen-vector/Principle Components Analysis (PCA) on your image data set. The procedure is shown as follows:
   
   • Find the mean of all input image vectors and make Mz a matrix where each column is an image minus the mean image (meanIm). This can be achieved by the built in Matlab function mean, but pay attention to which dimension it takes the mean over. We will synthesize/interpolate the intensity variation of a pixel over time, so the mean should be over the different image instances of the same pixel, not the mean over the pixel intensities in one image.
     Hint: Use an outer product of a mean column vector and a row vector of ones to form a matrix of the same size as Mz with the mean image repeated in each column. Optionally, visualize the acquired mean image and denote its association to the image dataset. Now you can directly subtract this mean matrix from M.
   • Calculate the covariance matrix C = (Mz' * Mz) / m;
   • Get the first K eigenvectors Y of C using eigs in matlab. (K=6 will be needed below, but for better quality try K=10 for the light, and K=20 for the objects. )
   • Calculate the basis vectors as B = Mz * Y;
   • The matrix of sample values for our coefficient functions is Y' (Mz=BY') (Pay attention to the " ' "-transposes, and note that Y in Exc 4 is transposed (because that is the way it comes out of eigs) compared to the Y in Exc 2,3,5)


3. Transform the resulting basis vectors (columns of B) into images of the same size as the input images, and display the first 5 (use absolute value). (Make sure you get the eigen vectors and images corresponding to the largest eigenvalues.) Plot a graph of the corresponding coefficients (rows of Y'); these are the functions we want to interpolate using polynomials (one poly for each row of Y').
Results should look something like these (from images of a rotating cat):

4. Then test that you can regenerate each input image k from your approximation as It = meanIm + B*Y(k,:)'. Reshape and show the images using your script from Exercise 4. (Alternatively you can add meanIm as follows: It = [meanIm;B]*[1;Y(k,:)'], where "1" is a row vector of ones. This is how the B matrix in the downloadable examples is done B = [meanIm;B])
5. Instead of using normalized PCA, perform SVD over the matrix M. You can use MATLAB's built-in SVD method to acquire the decompositions U, S, V. Again, keep only the first K eigenvectors V, and the first K eigenvalues S, by matrix slicing. As a result, the matrix of sample values for our coefficient functions Y = V(:,1:K)';', where we keep the first K eigenvectors, and the basis vectors can be computed as B = U * S(:,1:K); . Try to compare the quality of the constructed basis between both PCA and SVD and denote which method is better for approximation.