ComputerVision Lab 4 -Structure from Motion Solved

This week, the lab introduces structure from motion - inferring 3D shape from 2D images. Points matches across frames are saved into a matrix and the Tomasi-Kanade factorization is used for performing SfM (Structure from Motion).

1       Finding point correspondences
For finding point correspondences between multiple images (Not necessarily from the same camera) there are two possible options that we will consider in this course.

1.    Extracting SIFT descriptors and matching these descriptors at keypoint locations followed by a chaining step. This will be discussed in Lab Assignment 6.

2.    Using Optical Flow to estimate the apparent motion vectors: vx,vy of image pixels or regions from one frame to the next. The theoretical details of Optical Flow computations will be treated in a coming lecture and are not essential for successfully completing this lab assignment.

For this lab session, we provide a script that implements the Lucas-Kanade Optical Flow-based tracker: <LKtracker.m>.

Note: You are not expected to do any work with this script for optical flow based tracking. If you are interested, you can have a look at the implementation, as extra information for the Optical Flow lecture.

For this lab we will provide:

•    The point correspondences obtained from the Lucas-Kanade Optical Flow-based tracker as mat files ( ‘Xpoints’ & ‘Ypoints’). You can use these to perform the Structure-from-motion assignment.

•    We also provide measurement matrix.txt point correspondences obtained using SIFT matches and chaining as we will do in Lab 6.

2        Structure from Motion
Using the feature points from the provided measurement matrix or the provided Xpoints and Ypoints mat file containing Optical Flow tracked points, as input, perform the affine structure from motion procedure.

The following steps will guide you through factorization for structure from motion :

1.    After you load the point correspondences, normalize them by subtracting the centroid (mean).

 

• Hint : Mean should be taken across different views of each point. So, you need the mean of each column.

2.    The matrix of measurements, D needs to be decomposed into two matrices: M and S.

Apply SVD to the 2m × n loaded points matrix to express it as D = U ∗ W ∗ V 0.

 

 

1
[... ,... ,... ] = svd(X);
The matrix D must have rank 3. Keep only the top 3 columns/rows such that:

1
% The matrix of measurements must have rank 3.
2
% Keep only the corresponding top 3 rows/columns ...

from the SVD decompositions.
3
U = ...
4
W = ...
5
V = ...
•    U is a 2m × 3 matrix,

•    W is a 3 × 3 matrix,

•    V is a n × 3 matrix.

3.    From U,W,V compute the matrices M and S. A possible decomposition is  and 

1
% Compute M and S: One possible decomposition is ...

M = U Wˆ{1/2} and S = Wˆ{1/2} V'
2
M = ...
3
S = ...
4
save('M','M')
4.    This decomposition is not unique, to resolve the affine ambiguity, we use Least Squares and the Cholesky decomposition. For this we:

•    We define a matrix L such that AiLAi = Id, where Ai are the x and y coordinates of the motion matrix M, and Id is the identity matrix.

1
% We define the starting value for L, L0 as: ...

A1 L0 A1' = Id
2
A1 = M(1:2, :);
3
L0 = pinv(A1' * A1);
4

5
% We solve L by iterating through all images ...

and finding L one which minimizes ...

Ai*L*Ai' = Id, for all i.
6
% LSQNONLIN solves non-linear least squares ...

problems. Please check the Matlab ... documentation.
7
L = lsqnonlin(@residuals, L0);
We use the Matlab lsqnonlin function which minimizes a given function (In our case the function handler: @residuals) starting from an initial solution guess, in our case L0.

We then need to define the function <residuals.m> to be optimized by the lsqnonlin: AiLAi = Id.

1
% Compute the residuals
2

3
for i = 1:size(M, 1)/2 % Loop over the ... cameras: rows are ordered as x1,y1, ... x2,y2, ..
4
% Define the x and y projections for the ...

current camera:
5
Ai = M(i*2-1 : i*2, : );
6

7
% Definite the function to be minimized: ...

Ai L Ai' - Id = 0
8
diff  i = ...
9

10
diff(i,:) = diff  i(:);
11
end
•    We decompose L using the Cholesky decomposition into L = CCT.

•    We redefine M and S as: M = MC and S = C−1S.

1
% Recover C from L by Cholesky decomposition.
2
C = chol(L,'lower');
3

4
% Update M and S with the corresponding C ...

form: M = MC and S = Cˆ{-1}S.
5
M = ...
6
S = ...
5.    With the now estimated ‘S’ matrix, you can use plot3 to visualize the points. The rows of S are the 3D coordinates (x, y and z) vectors needed as the argument to plot3.

6.    Perform these steps using the given points in the measurement matrix.txt and the points obtained from the Optical Flow tracker (Xpoints and Ypoints matrices). Compare the obtained results.