• Before we introduced the conversion from Euclidean coordinate to Homogeneous coordinate:
  On image plane:

  \[ \begin{bmatrix} x\\y \end{bmatrix}\Rightarrow \begin{bmatrix} x\\y\\1 \end{bmatrix} \]

  In 3D physical space:

  \[ \begin{bmatrix} x\\y\\z \end{bmatrix}\Rightarrow \begin{bmatrix} x\\y\\z\\1 \end{bmatrix} \]
• Here we introduce a new rule to convert from Homogeneous coordinate to Euclidean coordinate:
  On image plane:

  \[ \begin{bmatrix} x\\y\\w \end{bmatrix}\Rightarrow \begin{bmatrix} x/w\\y/w \end{bmatrix} \]

  In 3D physical space:

  \[ \begin{bmatrix} x\\y\\z\\w \end{bmatrix}\Rightarrow \begin{bmatrix} x/w\\y/w\\z/w \end{bmatrix} \]

• The common practice is to use another 5-parameter representation of $K$: \[ K= \begin{bmatrix} f_x & s & c_x\\ 0 & f_y & c_y\\ 0 & 0 & 1 \end{bmatrix} \]
• In practice, $K$ may be accessed by the SDK of cameras
  • For example, here is a StackOverflow post that discusses extracting the intrinsic camera matrix of the ARKit of Apple
• We will also introduce algorithms to estimate $K$ in subsequent lectures

• Recall the projection from camera frame to image plane by intrinsic camera matrix: \[ P_h'=K[I, 0]P_h \]
• We just derived that \[ P_h = \begin{bmatrix} R & t \\ 0 & 1 \end{bmatrix}P_w \]
• Composed together, we can transform from world frame to image plane: \[ P_h'=K \begin{bmatrix} I & 0 \end{bmatrix} \begin{bmatrix} R & t\\ 0 & 1 \end{bmatrix}P_w= K \begin{bmatrix} R & T \end{bmatrix} P_w \]

• $M=K \begin{bmatrix} R & t \end{bmatrix}\in\R{3\times 4}$ includes all the parameters of a pinhole camera to form images!
• We refer to $M$ as the camera matrix
• While $M$ has 12 numbers, not all matrices in $\R{3\times 4}$ are valid camera matrices
• We can use the 5 parameters in $K$ and 6 parameters in $[R,t]$ to generate all valid $M$'s
• In literature, the common expression is $M$ has 11 degree of freedoms

• For the convenience of drawing figures, we introduce the virtual image
• It is the mirror reflection of the retina image
• The projective transformation is unchanged: \( P_h'= K \begin{bmatrix} R & t \end{bmatrix} P_w \)

• This is also a least square problem
• But the function inside $\|\cdot\|^2$ is not a linear function, because $f$ is non-linear
• We call it a non-linear optimization problem

• Similarly, let $L(Q)=\|p_1-f(Q;M_1)\|^2+\|p_2-f(Q;M_2)\|^2$
• We can optimization the objective function by following $\nabla L$ using gradient descent
• In fact, there are many smarter ways to adjust the gradient direction, so that it works faster
• For example, LM algorithm (Levenberg-Marquardt) is a popular choice

1. Take photos from many views
2. Identify points of interest from images
3. Search for corresponding points in other images
4. Compute camera positions
5. Estimate 3D point positions

• Advanced approach: Extract a point feature (e.g., by neural networks)
• But let us show a very basic approach first
  1. Extract a patch around a point
  2. Compare two patches $W_1$ and $W_2$ by pixel-wise distance:
     • SSD (Sum Squared Distance) \[ \sum_{x,y}|W_1(x,y)-W_2(x,y)|^2 \]
     • NCC (Normalized Cross Correlation) \[ \frac{\sum_{x,y}(W_1(x,y)-\bar{W}_1)(W_2(x,y)-\bar{W}_2)}{\sigma_{W_1}\sigma{W_2}} \]
• Anyway, global search requires exhaustive comparison, thus would be slow.

• The intrinsic camera matrices $K_1$ and $K_2$
• and the relative pose between the cameras $(R, T)$

• If we can have some pairs of corresponding points, we would be able to estimate $F$
• $F$ can help us to find more pairs of correspondences, so that we reconstruct the 3D for more points

• Given correspondences, $W$ is a constant
• The equation is linear w.r.t $f$
• The r.h.s. of the equation is $0$. The system is called a homogeneous system
• Observe that, if $f$ is a solution, then $\lambda f\ \forall \lambda\in\R{}$ is also a solution.

• Using the language of our matrix theory, we say that $f\in \text{null}(W)$
• Recall that each row of $W$ is from a corresponding point pair: \[ W= \begin{bmatrix} x_1x'_1& x_1y'_1&x_1&y_1x'_1&y_1y'_1&y_1&x'_1&y'_1&1\\ x_2x'_2& x_2y'_2&x_2&y_2x'_2&y_2y'_2&y_2&x'_2&y'_2&1\\ x_3x'_3& x_3y'_3&x_3&y_3x'_3&y_3y'_3&y_3&x'_3&y'_3&1\\ x_4x'_4& x_4y'_4&x_4&y_4x'_4&y_4y'_4&y_4&x'_4&y'_4&1\\ x_5x'_5& x_5y'_5&x_5&y_5x'_5&y_5y'_5&y_5&x'_5&y'_5&1\\ x_6x'_6& x_6y'_6&x_6&y_6x'_6&y_6y'_6&y_6&x'_6&y'_6&1\\ x_7x'_7& x_7y'_7&x_7&y_7x'_7&y_7y'_7&y_7&x'_7&y'_7&1\\ x_8x'_8& x_8y'_8&x_8&y_8x'_8&y_8y'_8&y_8&x'_8&y'_8&1\\ \end{bmatrix}\in\R{8\times 9} \]
• So $\rank(W)=8$ (when points are at general positions), and $\dim(\null(W))=9-8=1$.

1. $W\in\R{8\times 9}$
2. $\dim(\null(W))=1$
3. $f\in \text{null}(W)$

• $W$'s null space is a line in the 9-dimensional space, and $f$ is on this line.
• In fact, any $f\in\null(W)$ might be the groundtruth $f$
• In other words, 8 points have uniquely determined $f$ up to scaling!
• That's why we have the 8-point algorithm.
• Python code to find the null space:

  
  >>> from scipy.linalg import null_space
  >>> f = null_space(W)
                              

• While 8 points are sufficient to estimate $f$ theoretically, in practice, we use many more pairs.
• This would allow us to estimate $f$ more robustly w.r.t. inaccurate correspondences.
• Practically, this is very useful and important.
• However, when there are $n>8$ point pairs, $\rank(W)=9$ and $\null(W)=0$.
• The only solution is $0$!

• We add a constraint that $\|f\|_2^2=1$ to avoid the trivial solution $f=0$
• We convert the homogeneous equation to a least square objective
• A quadratic constraint quadratic programming (QCQP) problem

• Denote $L(x,\lambda)=f(x)+\lambda g(x)$
• The optimality condition is:

  $\nabla_x L(x, \lambda)=0$

  $\Downarrow$

  \begin{align} \nabla f(x)-\lambda \nabla g(x)=0 \end{align}

• By the definition of eigenvector and eigenvalue, $f$ must be an eigenvector and $\lambda$ must be an eigenvalue
• Note that $W^TW$ is a positive definite matrix
  • positive definite matrix $A$: $x^TAx>0\ \forall x\neq 0$
• Therefore, $\lambda$ must be the smallest eigenvalue and $f$ be the corresponding eigenvector

• An eigenvector $x$ of a linear transformation $A$ is a non-zero vector that, when $A$ is applied to it, does not change direction: \[ Ax = \lambda x, x \neq 0. \]
Applying $A$ to the vector only scales the vector by the scalar value $\lambda$, called an eigenvalue.

• We can reshape $f$ to be the $3\times 3$ matrix $F$
• $F\in\R{3\times 3}$ should be a rank-2 matrix
• However, the fundamental matrix obtained by reshaping $f$ is not rank-2. We need some post-processing.
• Assume $f$ is reshaped to a $3\times 3$ matrix $\hat{F}$. The post-processing is to solve an optimization problem: \begin{align} \begin{aligned} &\underset{}{\text{minimize}}&&\|F-\hat{F}\|^2_{F}\\ &\text{subject to}&&\rank(F)=2 \end{aligned} \end{align} where $\|A\|_F=\sqrt{\sum_{ij}A^2_{ij}}$

• This problem is called low-rank approximation
• The solution is given by an algorithm called SVD decomposition: \[ F=U \begin{bmatrix} s_1 & 0 & 0 \\ 0 & s_2 & 0 \\ 0 & 0 & 0 \end{bmatrix}V^T \] where \[ \hat{F}=U \begin{bmatrix} s_1 & 0 & 0 \\ 0 & s_2 & 0 \\ 0 & 0 & s_3 \end{bmatrix}V^T \text{ such that } U^TU=I, V^TV=I, s_i>0 \text{ sorted in descending order} \]

U, S, V = np.linalg.svd(Fhat)
print('U=', U)
print('S=', S)
print('V=', V)

S[2] = 0
F = U @ np.diag(S) @ V
                    

• When $x_i$ and $y_i$ are too large or too small, $W$ will have bad numerical properties.
• We often normalize the scale of $x_i$ and $y_i$ so that they are not too or too large.
• We skip details of step.

1. Pre-processing point pair coordinates;
2. Build $W$;
3. Estimate $f$ as the eigenvector of the smallest eigenvalue of $W^TW$;
4. Build $\hat{F}$ by reshaping $f$;
5. Use SVD to obtain the rank-2 approximation of $F$ from $\hat{F}$.

• When we know correspondences with accuracy, we can estimate $F$
• When we know $F$, we can restrict the search space and confirm correspondences

for i in range(n):
    randomly choose some pairs
    repeat for m times:
        based on the inliers, estimate F
        based on F, remove pairs with big errors
                    

1. RANSAC starts from a subset of correspondences, which is randomly chosen. This random initial set affects the final result.
2. For every iteration, we determine the inlier and outlier based on the level of errors. This is usually done by a threshold $\tau$, which affects the final result.
3. To obtain good results, we often try several random initial set, and try different thresholds.

1. Take photos from many views
2. Identify points of interest from images
3. Search for corresponding points in other images
4. Compute camera positions
5. Estimate 3D point positions

• the intrinsic camera matrices $K_1$ and $K_2$
• the external camera matrices $[R, t]$

• $E=[t]_{\times} R$ is also known as Essential matrix
• If $K_1$ and $K_2$ are known, we will be able to compute $E$ from $F$: \[ E=K_1^T FK_2 \]

• Note that $[t]_{\times}$ and $R$ both have special structures
  • $[t]_{\times}$: A skew-symmetric matrix
  • An orthonormal matrix
• Mathematically, given $E$, its decomposition into the product of a skew-symmetric matrix and an orthonormal matrix is almost unique!
  • Through SVD algorithm, as well.

1. Take photos from many views
2. Identify points of interest from images
3. Search for corresponding points in other images
4. Compute camera positions
5. Estimate 3D point positions

• We might recognize the point by looking through a small window
• We want a window shift in any direction to give a large change in intensity

• We have been treating an image as a matrix
• In CV, we also often treat an image as a function:
  • $I(x,y)$ returns the pixel value at $(x,y)$
• This view allows us to compute image gradients:
  • Use finite difference to approximate partial derivative: \[ \begin{aligned} \frac{\partial f}{\partial x}(x,y) & \approx f(x+1, y)-f(x, y)\\ \frac{\partial f}{\partial y}(x,y) & \approx f(x, y+1)-f(x, y)\\ \end{aligned} \]
  • Will revisit this topic in the future

• We want to discover how $E$ behaves for small shifts
• But this is very slow to compute naively: $\mathcal{O}(\text{window_width}^2*\text{shift_range}^2*\text{image_width}^2)$
• $\mathcal{O}(11^2*11^2*600^2)=5.2 \text{ billion of these }14.6\text{ thousand per pixel in your image}$

• Ellipses: $\lambda_1>0$, $\lambda_2>0$
• Parallel lines: $\lambda_1\lambda_2=0$
• Hyperbolic: $\lambda_1\lambda_2$ is less than $0$

1. Take photos from many views
2. Identify points of interest from images
3. Search for corresponding points in other images
4. Compute camera positions
5. Estimate 3D point positions