Chapter 14 Analytic Photogrammetry

14.1 Introduction Analytic photogrammetry includes the body of techniques by which, from

measurements of one or more 2D-perspective projections of a 3D object, one can make inferences about the 3D position, orientation, and lengths of the observed 3D object parts in a world reference frame.

inference problems: can be construed as nonlinear least-squares problems Photogrammetry provides a collection of methods for determining the position

and orientation of cameras and range sensors in the scene and relating camera positions and range measurements to scene coordinates.

GIS: Geographic Information System GPS: Global Positioning System

exterior orientation: determines position and orientation of camera in absolute coordinate system from projections of calibration points in scene

exterior orientation of camera: specified by all parameters of camera pose pose parameters: perspectivity center position, optical axis direction... exterior orientation specification: requires 3 rotation angles, 3 translations interior orientation: determines internal geometry of camera interior orientation of camera: parameters determining geometry of 3D rays

from measured image coordinates The parameters of interior orientation relate the geometry of ideal

perspective projection to the physics of a camera. parameters: camera constant, principal point, lens distortion .... complete specification of camera orientation: interior and exterior orientation

relative orientation: determines relative position and orientation between

2 cameras from projections of calibration points in scene relative orientation: to calibrate relation between two cameras for stereo relative orientation: relates coordinate systems of two cameras to each other,

not knowing 3D points themselves, only their projections in image relative orientation of one camera to another constitutes stereo model:

specified by 5 parameters: 3 rotation angles, 2 translations relative orientation: assumes interior orientation of each camera known

absolute orientation: determines transformation between 2 coordinate systems or position and orientation of range sensor in absolute coordinate system from coordinates of calibration points absolute orientation: to convert depth measurements in viewer-centered

coordinates to absolute coordinate system for the scene absolute orientation: orientation of stereo model in world reference frame absolute orientation: determines scale, 3 translations, 3 rotations absolute orientation: recovery of relation between two coordinate systems

( \begin{matrix}x \\ y \\ z \end{matrix} ): a point in the world reference frame

( \begin{matrix}x_{0} \\ y_{0} \\ z_{0} \end{matrix} ): camera lens position

left-hand coordinate system: clockwise angle of rotation ω: angle of rotation around x-axis of camera reference frame ω: tilt angle ϕ: angle of rotation around y-axis ϕ: pan angle κ: angle of rotation around z-axis κ: swing angle 3 × 3 rotation matrix:

R(ω) = ( \begin{matrix}1 & 0 & 0 \\ 0 & cos ω & sin ω \\ 0 & − sin ω & cos ω \end{matrix} )

R(ϕ) = ( \begin{matrix}cos ϕ & 0 & − sin ϕ \\ 0 & 1 & 0 \\ sin ϕ & 0 & cos ϕ \end{matrix} )

R(κ) = ( \begin{matrix}cos κ & sin κ & 0 \\ − sin κ & cos κ & 0 \\ 0 & 0 & 1 \end{matrix} )

R(ω, ϕ, κ) = R(κ)R(ϕ)R(ω) = ( \begin{matrix}r_{11} & r_{12} & r_{13} \\ r_{21} & r_{22} & r_{23} \\ r_{31} & r_{32} & r_{33} \end{matrix} ) =

cos ϕ cos κ sin ω sin ϕ cos κ + cos ω sin κ − cos ω sin ϕ cos κ + sin ω sin κ
− cos ϕ sin κ − sin ω sin ϕ sin κ + cos ω cos κ cos ω sin ϕ sin κ + sin ω cos κ
sin ϕ  − sin ω cos ϕ  cos ω cos ϕ

angles can be obtained directly from values r_{ij}:

sin ϕ = r_{31}

tan ω = (−r_{32})/r_{33}

tan κ = (−r_{21})/r_{11}

point (x, y, z)′ in world frame represented by (p, q, s)′ in camera frame: ′: denotes matrix or vector transpose R^{−1} = R′

( \begin{matrix}p \\ q \\ s \end{matrix} ) = R(ω, ϕ, κ) ( \begin{matrix}x − x_{0} \\ y − y_{0} \\ z − z_{0} \end{matrix} )

pinhole camera with image at distance f from camera lens, projection:

( \begin{matrix}u \\ v \end{matrix} ) = \frac{f}{s} ( \begin{matrix}p \\ q \end{matrix} )

f: camera constant, related to focal length of lens

principal point: origin of measurement image plane coordinate principal point (u_{0}, v_{0}): in image measurement plane coordinate system

( \begin{matrix}u \\ v \end{matrix} ) = ( \begin{matrix}u_{0} \\ v_{0} \end{matrix} ) + \frac{f}{s} ( \begin{matrix}p \\ q \end{matrix} )

general fundamental perspective projection equations: collinearity equation:

\frac{x − x_{0}}{z − z_{0}} = \frac{r_{11}(u − u_{0}) + r_{21}(v − v_{0}) + r_{31}f}{r_{13}(u − u_{0}) + r_{23}(v − v_{0}) + r_{33}f}

\frac{y − y_{0}}{z − z_{0}} = \frac{r_{12}(u − u_{0}) + r_{22}(v − v_{0}) + r_{32}f}{r_{13}(u − u_{0}) + r_{23}(v − v_{0}) + r_{33}f}

They show that the relationship between the measured 2D-perspective projection coordinates and the 3D coordinates is a nonlinear function of u_{0}, v_{0}, x_{0}, y_{0}, z_{0}, ω, ϕ, and κ. =====Garfield 17:27=====

14.2 Nonlinear Least-Squares Solutions noise model:

α_{k} = g_{k}(β_{1}, ..., β_{M}) + ξ_{k}, k = 1, ..., K

β_{1}, ..., β_{M}: unknown parameters governing nonlinear transformations g_{1}, ..., g_{K}: nonlinear transformations α_{1}, ..., α_{K}: observed values of g_{1}, ..., g_{k} ξ_{1}, ..., ξ_{K}: additive mean zero Gaussian random variables with covariance Σ̸

maximum likelihood solution: β_{1}, ..., β_{M} maximize P rob(α_{1}, ..., α_{K}|β_{1}, ..., β_{M}) maximum likelihood solution minimizes least-squares criterion

ϵ^{2} = (α − g)′ Σ̸^{−1}(α − g)

where

α = ( α_{1} ; α_{2} ; ⋮ ; α_{K} ) and g = ( g_{1}(β_{1}, ..., β_{M}) ; g_{2}(β_{1}, ..., β_{M}) ; ⋮ ; g_{K}(β_{1}, ..., β_{M}) )

assumption: overconstrained and unique solution, K > M

first-order Taylor series expansion of g_{k} taken around β^{t}

g_{k}(β^{t} + Δβ) = g_{k}(β^{t}) + Δg_{k}(Δβ; β^{t})

β^{t} = (β_{1}^{t}, ..., β_{M}^{t})′: tth iteration approximate solution Δβ = (Δβ_{1}, ..., Δβ_{M})′: adjustments added for a better approximate solution Δg_{k}: total derivative of g_{k}: linear function of adjustment vector Δβ

Δg_{k}(Δβ; β^{t}) = (\frac{∂g_{k}}{∂β_{1}}(β^{t}), ..., \frac{∂g_{k}}{∂β_{M}} (β^{t})) Δβ

G^{t}: Jacobian given by

G^{t} = ( \begin{matrix}\frac{∂g_{1}}{∂β_{1}} & \frac{∂g_{1}}{∂β_{2}} & ⋯ & \frac{∂g_{1}}{∂β_{M}} \\ \frac{∂g_{2}}{∂β_{1}} & \frac{∂g_{2}}{∂β_{2}} & ⋯ & \frac{∂g_{2}}{∂β_{M}} \\ ⋮ & ⋮ & ⋮ & ⋮ \\ \frac{∂g_{K}}{∂β_{1}} & \frac{∂g_{K}}{∂β_{2}} & ⋯ & \frac{∂g_{K}}{∂β_{M}} \end{matrix} )

adjustments added to the current β to form the next β:

β^{t+1} = β^{t} + Δβ = β^{t} + (G^{t′} Σ̸^{−1}G^{t})^{−1}G^{t′} Σ̸^{−1}(α − g^{t})

14.3 The Exterior Orientation Problem one-camera exterior orientation problem: known:

• (x_{n}, y_{n}, z_{n})′: set of N 3D points with known positions, n = 1, ..., N

• (u_{n}, v_{n})′: corresponding set of 2D-perspective projections, n = 1, ..., N unknown:

• rotation (ω, ϕ, κ) and translation (x_{0}, y_{0}, z_{0})′ that put the camera reference frame in the world reference frame

camera calibration problem: unknown position of camera in object frame object pose estimation problem: unknown object position in camera frame spatial resection problem in photogrammetries: 3D positions from 2D orientation (p_{n}, q_{n}, s_{n}): point (x_{n}, y_{n}, z_{n})’s coordinate in camera reference frame

( \begin{matrix}p_{n} \\ q_{n} \\ s_{n} \end{matrix} ) = R(ω, ϕ, κ) ( \begin{matrix}x_{n} − x_{0} \\ y_{n} − y_{0} \\ z_{n} − z_{0} \end{matrix} ) , n = 1, ..., N

observed perspective projection

( \begin{matrix}u_{n} \\ v_{n} \end{matrix} ) = \frac{f}{s_{n}} ( \begin{matrix}p_{n} \\ q_{n} \end{matrix} ) , n = 1, ..., N

in the notation of the previous section, β = (x_{0}, y_{0}, z_{0}, ω, ϕ, κ)

g_{2n−1}(β) = u_{n} = fp_{n}/s_{n}
g_{2n}(β) = v_{n} = fq_{n}/s_{n}

14.3.1 Standard Solution by chain rule for partial differentiation

Δu_{n} = f \frac{s_{n}Δp_{n} − p_{n}Δs_{n}}{s_{n}^{2}}

Δv_{n} = f \frac{s_{n}Δq_{n} − q_{n}Δs_{n}}{s_{n}^{2}}

in matrix form

( \begin{matrix}Δu_{n} \\ Δv_{n} \end{matrix} ) = \frac{f}{s_{n}} ( \begin{matrix}1 & 0 & −p_{n}/s_{n} \\ 0 & 1 & −q_{n}/s_{n} \end{matrix} ) ( \begin{matrix}Δp_{n} \\ Δq_{n} \\ Δs_{n} \end{matrix} )

=====p. 134===== ..... (tedious mathematics, study by yourself)

14.3.2 Auxiliary Solution Instead of iteratively adjusting the angles directly, we can reorganize the calculation such that we iteratively adjust the three auxiliary parameters of a skew symmetric matrix associated with the rotation matrix. Then we determine the adjustments of the angles in terms of the adjustments of the parameters of the skew symmetric matrix associated with the rotation matrix.

14.3.3 Quaternion Representation quaternion: one other representation for rotation matrices from any skew symmetric matrix

S = ( \begin{matrix}0 & −c & b \\ c & 0 & −a \\ −b & a & 0 \end{matrix} )

a rotation matrix R can be constructed by choosing scalar d

R = (dI + S)(dI − S)^{−1}

this definition guarantees that R′R = I:

R′R = [(dI + S)(dI − S)^{−1}]′[(dI + S)(dI − S)^{−1}] = ⋯ = I

expanding the equation for R, we have R =

( \begin{matrix}d^{2} + a^{2} − b^{2} − c^{2} & 2(ab − cd) & 2(ac + bd) \\ 2(ab + cd) & d^{2} − a^{2} + b^{2} − c^{2} & 2(bc − ad) \\ 2(ac − bd) & 2(bc + ad) & d^{2} − a^{2} − b^{2} + c^{2} \end{matrix} ) \frac{1}{a^{2} + b^{2} + c^{2} + d^{2}}

parameters a, b, c, and d can be constrained to satisfy

a^{2} + b^{2} + c^{2} + d^{2} = 1

in terms of quaternion parameters

ΔR = \frac{∂R}{∂a}Δa + \frac{∂R}{∂b}Δb + \frac{∂R}{∂c}Δc + \frac{∂R}{∂d}Δd

so that

ΔRR′ = \frac{∂R}{∂a}R′Δa + \frac{∂R}{∂b}R′Δb + \frac{∂R}{∂c}R′Δc + \frac{∂R}{∂d}R′Δd

calculating the matrices directly, we obtain

\frac{∂R}{∂a}R′ = 2 ( \begin{matrix}a & b & c \\ −b & a & −d \\ −c & d & a \end{matrix} )

..... (tedious mathematics, study by yourself) =====Oldie 33:76=====

14.4 Relative Orientation The transformation from one camera station to another can be represented by a rotation and a translation. The relation between the coordinates, r_{l} and r_{r} of a point P can be given by means of a rotation matrix and an offset vector. =====Horn, Robot Vision, Fig. 13.3===== Relative orientation is typically concerned with the determination of the position and orientation of one photograph with respect to another, given a set of corresponding image points. (x_{L}, y_{L}, z_{L}), (x_{R}, y_{R}, z_{R}): 3D positions of left and right camera lenses (ω_{L}, ϕ_{L}, κ_{L}), (ω_{R}, ϕ_{R}, κ_{R}): rotation angles specifying exterior orientation {(u_{Ln}, v_{Ln})}_{n=1}^{N}: set of N image points from left image {(u_{Rn}, v_{Rn})}_{n=1}^{N}: corresponding set of N image points from right image separation of two camera lenses (x_{R} − x_{L}): as constant controlling scale relative orientation: specified by five parameters

(y_{R} − y_{L}), (z_{R} − z_{L}), (ω_{R} − ω_{L}), (ϕ_{R} − ϕ_{L}), (κ_{R} − κ_{L}) assumption:

• camera interior orientation known

• image positions expressed to identical scale and w.r.t. principal point

14.4.1 Standard Solution rotation matrix associated with exterior orientation of left image

Q′_{L} = R(ω_{L}, ϕ_{L}, κ_{L})

rotation matrix associated with exterior orientation of right image

Q′_{R} = R(ω_{R}, ϕ_{R}, κ_{R})

f_{R}: distance between right image plane and right lens f_{L}: distance between left image plane and left lens from perspective collinearity equation

( \begin{matrix}u_{Ln} \\ v_{Ln} \\ f_{L} \end{matrix} ) = \frac{1}{λ_{Ln}}Q′_{L} ( \begin{matrix}x_{n} − x_{L} \\ y_{n} − y_{L} \\ z_{n} − z_{L} \end{matrix} ) and ( \begin{matrix}u_{Rn} \\ v_{Rn} \\ f_{R} \end{matrix} ) = \frac{1}{λ_{Rn}}Q′_{R} ( \begin{matrix}x_{n} − x_{R} \\ y_{n} − y_{R} \\ z_{n} − z_{R} \end{matrix} ) and

Hence

( \begin{matrix}x_{n} \\ y_{n} \\ z_{n} \end{matrix} ) = ( \begin{matrix}x_{L} \\ y_{L} \\ z_{L} \end{matrix} ) + λ_{Ln}Q_{L} ( \begin{matrix}u_{Ln} \\ v_{Ln} \\ f_{L} \end{matrix} ) = ( \begin{matrix}x_{R} \\ y_{R} \\ z_{R} \end{matrix} ) + λ_{Rn}Q_{R} ( \begin{matrix}u_{Rn} \\ v_{Rn} \\ f_{R} \end{matrix} )

where λ_{Ln} = \frac{s_{Ln}}{f_{L}} and λ_{Rn} = \frac{s_{Rn}}{f_{R}} ..... (tedious mathematics, study by yourself)

14.4.2 Quaternion Solution Hinsken sets up the problem in a slightly different way. Instead of determining the relative orientation of the right image with respect to the left image, he aligns a reference frame having its x-axis along the line from the left image lens to the right image lens. The relative orientation is then determined by the angles (ω_{R}, ϕ_{R}, κ_{R}), which rotate the right image into this reference frame, and the angles (ω_{L}, ϕ_{L}, κ_{L}), which rotate the left image into this reference frame.

14.5 Interior Orientation interior orientation: inner orientation: of a camera is specified by

• camera constant f: distance between image plane and camera lens

• principal point (u_{p}, v_{p}): intersection of optic axis with image plane in measurement reference frame located on image plane

• geometric distortion characteristics of the lens; assuming isotropic around the principal point

14.6 Stereo optical axes parallel to one another and perpendicular to baseline simple camera geometry for stereo photography =====Horn, Robot Vision, Fig. 13.1===== stereo pair taken by Viking Lander I on the surface of Mars =====Horn, Robot Vision, Fig. 13.2===== parallax: displacement in perspective projection by position translation (x, y, z): 3D point position (u_{L}, v_{L}): perspective projection on left image of stereo pair (u_{R}, v_{R}): perspective projection on right image of stereo pair b_{x}: baseline length in x-axis position of the left camera:

( \begin{matrix}−b_{x}/2 \\ 0 \\ 0 \end{matrix} )

position of the right camera:

( \begin{matrix}b_{x}/2 \\ 0 \\ 0 \end{matrix} )

v_{L} = v_{R}, so that the y-parallax is zero

( \begin{matrix}u_{L} \\ v_{L} \end{matrix} ) = \frac{f}{z} ( \begin{matrix}x + b_{x}/2 \\ y \end{matrix} ) and ( \begin{matrix}u_{R} \\ v_{R} \end{matrix} ) = \frac{f}{z} ( \begin{matrix}x − b_{x}/2 \\ y \end{matrix} ) and

x-parallax: the difference u_{L} − u_{R}

u_{L} − u_{R} = \frac{f}{z} [x + \frac{b_{x}}{2} − (x − \frac{b_{x}}{2})] = \frac{f}{z} (b_{x})

Hence

z = \frac{fb_{x}}{u_{L} − u_{R}}

=====Example 14.3===== relation z = fb_{x}/(u_{L} − u_{R}) close to being useless in real-world, because

• observed perspective projections are subject to measurement errors so that v_{L} ≠ v_{R} for corresponding points

• left and right camera frames may have slightly different orientations

• when two cameras used, almost always f_{R} ≠ f_{L} derive least-squares solution as in the textbook

14.7 2D-2D Absolute Orientation The relationship between two coordinate systems is easy to find if we can measure the coordinates of a number of points in both systems. It takes three measurements to tie two coordinate systems together uniquely.

(a) A single measurement leaves three degrees of freedom of motion.

(b) A second measurement removes all but one degree of freedom.

(c) third measurement rigidly attaches two coordinate systems to each other. =====Horn, Robot Vision, Fig. 13.4===== 2D-2D absolute orientation problem: 2D-2D pose detection problem: to determine from matched points more precise estimate of rotation matrix R and translation t such that y_{n} = Rx_{n} + t, n = 1, ..., N: determine R and t that minimize weighted sum of residual errors:

ϵ^{2} = N ; ∑ ; n=1 w_{n}||y_{n} − (Rx_{n} + t)||^{2}

weights w_{n}, n = 1, ..., N: satisfy w_{n} ≥ 0 and ∑_{n=1}^{N} w_{n} = 1 if no prior knowledge: can be defined to be equal, w_{n} = 1/N ..... (tedious mathematics, study by yourself)

14.8 3D-3D Absolute Orientation p_{1}, ..., p_{N}: N 3D points relative to camera reference frame q_{1}, ..., q_{N}: N 3D points relative to object model reference frame we must determine rotation matrix R and translation vector t satisfying

p_{n} = Rq_{n} + t, n = 1, ..., N

constrained least-squares problem to minimize

N ; ∑ ; n=1 w_{n}||p_{n} − (Rq_{n} + t)||^{2}

The least-squares problem can be modeled by a mechanical system in which corresponding points in the two coordinate systems are attached to each other by means of springs. The solution to the least-squares problem corresponds to the equilibrium position of the system, which minimizes the energy stored in the springs. =====Horn, Robot Vision, Fig. 13.5===== Lagrange multipliers can be introduced and a solution found ..... (tedious mathematics, study by yourself) =====joke=====

14.9 Robust M-Estimation least-squares techniques are ideal when random data perturbations

or measurement errors are Gaussian distributed this section describes some robust techniques for nonlinear regression

14.9.1 Modified Residual Method r_{i}: residual of i-th datum r_{i}: difference between i-th observation and its fitted value standard least-squares method: tries to minimize ∑_{i} r_{i}^{2} standard least-squares method: unstable if outliears present in data outlying data: strong effect to distort estmated parameters M-estimators: try to reduce outlier effect by replacing squared residuals M-estimators: min ∑_{i} ρ(r_{i}) ρ: symmetric, positive-definite function with unique minimum at zero ρ: chosen to be less increasing than square =====p. 168===== =====p. 169=====

14.9.2 Modified Weights Method

14.9.3 Experimental Results

14.10 Error Propagation x_{1}, ..., x_{N}: input parameters Δx_{1}, ..., Δx_{N}: random errors quantity y depends on input parameters through known function f:

y = f(x_{1}, ..., x_{N})

quantities x_{1} + Δx_{1}, ..., x_{N} + Δx_{N} observed and quantity y + Δy calculated:

y + Δy = f(x_{1} + Δx_{1}, ..., x_{N} + Δx_{N})

error propagation analysis: determines expected value and variance of y + Δy known information about Δx_{1}, ..., Δx_{N}: mean, variance:

E [( Δx_{1} ; ⋮ ; Δx_{N} )] = 0

E [( Δx_{1} ; ⋮ ; Δx_{N} ) (Δx_{1}, ..., Δx_{N})] = ( \begin{matrix}σ_{11} & σ_{12} & ⋯ & σ_{1N} \\ ⋮ & ⋮ & ⋮ & ⋮ \\ σ_{N1} & σ_{N2} & ⋯ & σ_{NN} \end{matrix} )

14.10.1 Implicit Form error propagation: can also be done in implicit form known function f has the form

f(x_{1}, x_{2}, ..., x_{N}, y) = 0

The quantities (x_{1} + Δx_{1}, ..., x_{N} + Δx_{N}) are observed, and the quantity y + Δy is determined to satisfy f(x_{1} + Δx_{1}, ..., x_{N} + Δx_{N}, y + Δy) = 0.

14.10.2 Implicit Form: General Case general case: y is not a scalar but a L × 1 vector β x_{1}, ..., x_{K}: are K N × 1 vectors representing true values x_{1} + Δx_{1}, ..., x_{K} + Δx_{K}: are K N × 1 vectors representing noisy observed values Δx_{1}, ..., Δx_{K}: random perturbations β: is a L × 1 vector representing unknown true parameters noiseless model:

f(x_{k}, β) = 0, k = 1, ..., K

with noisy observations, the idealized model:

f(x_{k} + Δx_{k}, β + Δβ) = 0, k = 1, ..., K

14.11 Summary we have shown how to

• take a nonlinear least-squares problem

• linearize it

• solve by iteratively solving successive linearized least-squares problems =====joke=====

Project due April 19: camera calibration: interior orientation: to determine principal point do camera motion in z-axis, find focus of expansion/contraction use lens of focal length: 16mm, 25mm, 50mm