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 π(Cj,Pi) 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 Cj and Pi, 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:
Xc=RjPi+tj
- Normalize:
u=Xcx/Xcz,v=Xcy/Xcz
- Apply intrinsics (for a pinhole camera):
[fxu+cxfyv+cy]
The exact form depends on your parameterization (are intrinsics fixed? any distortion?).
- Residuals and Jacobian blocks
For each observed 2D point xij, residual is
rij=π(Cj,Pi)−xij∈R2
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
−Chainrule:
\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×3 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