📘

Jacobian for Bundle Adjustment

HardOptimization 📚 CV: Image Alignment and Stitching

Implement a Jacobian matrix computation for bundle adjustment optimization, a crucial step in image alignment and stitching. This process involves minimizing the reprojection error between observed and predicted image points.

Bundle adjustment is a non-linear least squares problem, where the goal is to find the optimal camera parameters and 3D point coordinates that minimize the reprojection error. The error is typically measured as the sum of squared differences between observed image points $x_{ij}$ and predicted points $\pi(C_j, P_i)$ , where $\pi$ is the projection function, $C_j$ represents the camera parameters, and $P_i$ represents the 3D point coordinates.

To compute the Jacobian matrix, we need to calculate the partial derivatives of the projection function with respect to the camera parameters and 3D point coordinates.

Compute the partial derivatives of the projection function $\pi$ with respect to the camera parameters $C_j$ .
Compute the partial derivatives of the projection function $\pi$ with respect to the 3D point coordinates $P_i$ .

E = \sum_{i,j} \|\pi(C_j, P_i) - x_{ij}\|^2

This technique is widely used in structure from motion and stereo vision applications.

Example:

Input:

camera = [500, 500, 320, 240, 0, 0, 0, 0, 0, 5]
point = [0, 0, 5]

Output:

2×13 Jacobian matrix

Reasoning:

Jacobian has derivatives of projected (u,v) w.r.t:

Intrinsics (fx, fy, cx, cy)
Rotation (rx, ry, rz)
Translation (tx, ty, tz)
3D point (X, Y, Z)

Computed via chain rule through projection.

Constraints:

camera: Camera parameters [fx, fy, cx, cy, rx, ry, rz, tx, ty, tz]
point: 3D point [X, Y, Z]
Return: 2×13 Jacobian [d(u,v)/d(camera, point)]

More from CV: Image Alignment and Stitching

Cylindrical Projection for PanoramasHard Essential Matrix EstimationHard Procrustes AnalysisMedium Gradient Descent UpdateMedium Linear Blend WeightEasy

Bundle adjustment is a nonlinear least squares problem where you optimize camera parameters and 3D point positions to minimize reprojection error across all images and points. The projection function $\pi(C_j, P_i)$ is nonlinear because it involves camera pose (rotation, translation), intrinsics, and a perspective divide, so you linearize it using a Jacobian and solve iteratively (e.g., with Gauss–Newton or Levenberg–Marquardt). The Jacobian encodes how each residual (2D reprojection error) changes with respect to each parameter (camera and point).

In practice, the Jacobian for bundle adjustment is very sparse and structured: each 2D observation only depends on one camera and one 3D point, so that residual row has nonzero entries only in the columns for that camera’s parameters and that point’s coordinates. This structure is crucial for efficiency and is what your problem is asking you to exploit: for each observation $(i,j)$ , compute the partial derivatives of the projection with respect to $C_j$ and $P_i$ , and place them into the right block locations in $J$ .

1. Background Knowledge (key concepts)

Projection model Typically:

Transform point from world to camera coordinates: $X_c = R_j P_i + t_j$
Normalize: $u = X_c^x / X_c^z,\quad v = X_c^y / X_c^z$
Apply intrinsics (for a pinhole camera):

\begin{bmatrix} f_x u + c_x \\ f_y v + c_y \end{bmatrix}

The exact form depends on your parameterization (are intrinsics fixed? any distortion?).

Residuals and Jacobian blocks For each observed 2D point $x_{ij}$ , residual is $r_{ij} = \pi(C_j, P_i) - x_{ij} \in \mathbb{R}^2$ The Jacobian row block for this residual w.r.t. all parameters is

\left[\frac{\partial \pi}{\partial C_j}\;\;\;\frac{\partial \pi}{\partial P_i}\right]$$ where $\theta$ concatenates all camera and point parameters. All other entries in that row are zero. - **Parameterizations** - Camera pose rotation is typically parameterized by **Rodrigues vector** (axis-angle, 3 params) or a minimal representation (3 DOF). - Camera translation: 3 params. - 3D point: 3 params $(X,Y,Z)$. - Intrinsics / distortion may also be parameters but are often fixed in basic tasks. --- ### 2. Algorithm / General Approach To compute the Jacobian for bundle adjustment: 1. **Fix a projection model and parameterization** for cameras and points. 2. For each observation $(i, j)$: - Compute the **intermediate quantities** (camera coordinates, normalized coordinates). - Compute **analytic partial derivatives**: - $\partial \pi / \partial X_c$ (2×3) - $\partial X_c / \partial C_j$ (3×dim($C_j$)) - $\partial X_c / \partial P_i$ (3×3) - Use **chain rule** to get: - $\partial \pi / \partial C_j$ - $\partial \pi / \partial P_i$ 3. Insert each 2×(camera-dim) block and 2×3 block into the correct positions in the global sparse Jacobian. The pattern is always: *compute projection → differentiate through steps → assemble sparse blocks.* --- ### 3. Step-by-Step Strategy (implementation-oriented) Assume: - Camera pose: $C_j = (r_j, t_j)$ with $r_j \in \mathbb{R}^3$ (Rodrigues), $t_j \in \mathbb{R}^3$. - Each 3D point: $P_i = (X_i, Y_i, Z_i)$. 1. **Indexing design** - Assign each camera $j$ a **camera parameter offset** in the parameter vector (e.g., 6 per camera if only pose). - Assign each point $i$ a **point parameter offset** (3 per point). - For each observation $(i, j)$, you know: - row indices in $J$: two consecutive rows for $r_{ij}^x$ and $r_{ij}^y$. - column ranges for camera $j$ and point $i$. 2. **Forward projection for one observation** ```python # Given camera j and point i R = rodrigues_to_rotation_matrix(r_j) # 3x3 Xc = R @ P_i + t_j # [Xc, Yc, Zc] u = Xc / Xc v = Xc / Xc x_proj = fx * u + cx y_proj = fy * v + cy # residual r = [x_proj - x_ij, y_proj - y_ij] ``` 3. **Derivatives w.r.t. camera coordinates $X_c$** - Let $X_c = (X, Y, Z)$. - Derivative of normalized coordinates:

\frac{\partial u}{\partial X_c} = \begin{bmatrix} 1/Z & 0 & -X/Z^2 \end{bmatrix},\quad \frac{\partial v}{\partial X_c} = \begin{bmatrix} 0 & 1/Z & -Y/Z^2 \end{bmatrix}

- Then:

\frac{\partial \pi}{\partial X_c} = \begin{bmatrix} f_x \frac{\partial u}{\partial X_c} \ f_y \frac{\partial v}{\partial X_c} \end{bmatrix} \in \mathbb{R}^{2\times 3}

4. **Derivatives w.r.t. 3D point $P_i$** - $X_c = R_j P_i + t_j$ so:

\frac{\partial X_c}{\partial P_i} = R_j

- Chain rule:

\frac{\partial \pi}{\partial P_i} = \frac{\partial \pi}{\partial X_c} \frac{\partial X_c}{\partial P_i} = \left(\frac{\partial \pi}{\partial X_c}\right) R_j

- This gives a 2\times3 block for point $i$. 5. **Derivatives w.r.t. camera pose $C_j$** - Translation:

\frac{\partial X_c}{\partial t_j} = I_{3\times 3} \Rightarrow \frac{\partial \pi}{\partial t_j} = \frac{\partial \pi}{\partial X_c} I = \frac{\partial \pi}{\partial X_c}

- Rotation (Rodrigues): - $X_c = R(r_j) P_i + t_j$. - Need $\partial X_c / \partial r_j$ (3×3). There is a standard closed-form expression for $\partial R / \partial r_k$; for this problem it is enough to know: - You can precompute **skew-symmetric matrix** $[P_i]_\times$. - Locally, a small rotation increment $\delta\omega$ gives $\delta X_c \approx -[X_c]_\times \,\delta\omega$. - So Jacobian w.r.t. rotation parameters often looks like:

\frac{\partial X_c}{\partial \delta \omega} \approx -[X_c]_\times

- Then:

\frac{\partial \pi}{\partial r_j} = \frac{\partial \pi}{\partial X_c} \frac{\partial X_c}{\partial r_j}

- Concatenate:

\frac{\partial \pi}{\partial C_j} = \left[ \frac{\partial \pi}{\partial r_j} ;;; \frac{\partial \pi}{\partial t_j} \right]

- This gives a 2×(pose-dim) block. 6. **Assemble the global Jacobian** - For each observation $(i, j)$: - At rows $2k, 2k+1$ (for this residual): - Insert $\partial \pi / \partial C_j$ into columns for

solution.py

Python 3.13

Test Results

0/0

Run code to see test results.