Point cloud registration with error propagation

Abstract

Discussed herein are devices, systems, and methods for merging point cloud data with error propagation. A can include reducing a sum aggregate of discrepancies between respective tie points and associated 3D points in first and second 3D images, adjusting 3D error models of the first and second 3D images based on the reduced discrepancies to generate registered 3D images, and propagating an error of the first or second 3D images to the registered 3D image to generate error of the registered 3D images.

Description

TECHNICAL FIELD

Embodiments discussed herein regard devices, systems, and methods for merging three-dimensional (3D) point clouds. Error of the merged 3D point cloud can be bounded by errors of the individual 3D point clouds.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 illustrates an example of an embodiment of a system for 3D point set merging.

FIG. 2 illustrates an example diagram of an embodiment of the relationship between ground point coordinate estimates and the corresponding 3D data set observations.

FIG. 3 illustrates an example of an embodiment of a bundle adjust operation.

FIG. 4 illustrates, by way of example, a block diagram of an embodiment of a machine in the example form of a computer system within which instructions, for causing the machine to perform any one or more of the methods discussed herein, may be executed.

DETAILED DESCRIPTION

Discussed herein are methods, systems, and devices for registering a first 3D point cloud (or a portion thereof) to one or more other 3D point clouds (or a portion thereof) to generate a merged 3D point cloud. One or more the first and second 3D point clouds can include an associated error. The associated error can be propagated to the merged 3D point cloud. The error of the 3D point cloud can be used in a downstream application. Example applications include targeting and mensuration. A targeteer (one who performs targeting) can benefit from the error to better inform their targeting location choice. A mensuration of an object can benefit from the error as well.

The merged 3D point clouds can include error that is better than either of the first and second 3D point clouds individually. For example, if the first 3D point cloud includes a lower error (relative to the second 3D point cloud) in the x and y directions and the second 3D point cloud includes a lower error (relative to the first 3D point cloud) in the z direction, the merged 3D point cloud can include error bounded by the first 3D point cloud in the x and y directions and by the second 3D point cloud in the z direction. The merged point cloud can thus inherit the better of the errors between the first and second point clouds for a specified parameter.

Reference will now be made to the FIGS. to describe the methods, systems, and devices for merging 3D point clouds with error propagation.

FIG. 1 illustrates an example of a system for 3D point set merging. The system can include processing circuitry 102 that receives tie points 108, tie point error 110, a first 3D point set 104, a second 3D point set 106, tie points 108, tie point error 110, and first or second point set error 112. The first or second point set error 112 includes error for at least one of the first 3D point set 104 and the second 3D point set 106. The first or second point set error 112 can thus include error for the first 3D point set 104, the second 3D point set 106, or the first 3D point set 104 and the second 3D point set 106.

The first 3D point set 104 or the second 3D point set 106 can include a point cloud, a 3D surface, or the like. The first 3D point set 104 and the second 3D point set 106 can include (x, y, z) data for respective geographic regions. The geographic regions of the first 3D point set 104 and the second 3D point set 106 at least partially overlap. One or more of the first point set 104 and the second point set 106 can include intensity data. Intensity data can include one or more intensity values, such as red, green, blue, yellow, black, white, gray, infrared, thermal, or the like. One or more of the first point set 104 and the second the point set 106 can include error data. The error data is illustrated as being a separate input in FIG. 1, namely the first or second point set error 112. The error data can indicate an accuracy of the corresponding point of the point set.

The tie points 108 can associate respective points between the first 3D point set 104 and the second 3D point set 106. The tie points 108 can indicate a first point (x₁, y₁, z₁) in the first 3D point set 104, a second point (x₂, y₂, z₂) in the second 3D point set 106 or an error associated with the tie point 108 (shown as separate input tie point error 110). The tie point error 110 can indicate how confident one is that the first and second points correspond to the same geographic location. The tie point error 110 can include nine entries indicating a covariance (variance or cross-covariance) between three variables. The three variables can be error in the respective directions (x, y, z). A covariance matrix representation of the tie point error 110 is provided as

$\hspace{1em}\left[\begin{array}{ccc}{e}_x& e_{xy}& e_{xz}\\ e_{yx}& e_y& e_{yz}\\ e_{zx}& e_{zy}& e_z\end{array}\right]$ where the diagonal terms are respective variances in the given directions, and the off-diagonal terms are covariances between the directions.

The first or second point set error 112 can sometimes be improved, such as to be more rigorous. Sometimes, the first or second point set error 112 can be in a form that is not digestible by the bundle adjustment operation 118. The point set error 112 can be conditioned by a condition point set error operation 114 to generate an error matrix 116. The condition point set error operation 114 can include generating a covariance matrix 116 of error parameters of the first 3D point set 104 or the second 3D point set 106. The error parameters can include seven parameters. Three of the parameters can include translation in x, y, and z, respectively. Three of the parameters can be for rotation in x, y, and z (roll, pitch, and yaw), respectively. One of the parameters can be for a scale factor between the first 3D point set 104 and the second 3D point set 106. An example of the matrix 116 produced by the condition point set error operation 114 is provided as

$\hspace{1em}\left[\begin{array}{ccccccc}\stackrel{\_}{x}& \stackrel{\_}{x}\stackrel{\_}{y}& \stackrel{\_}{x}\stackrel{\_}{z}& \stackrel{\_}{x}\omega & \stackrel{\_}{x}\phi & \stackrel{\_}{x}\kappa & \stackrel{\_}{x}s\\ \stackrel{\_}{y}\stackrel{\_}{x}& \stackrel{\_}{y}& \stackrel{\_}{y}\stackrel{\_}{z}& \stackrel{\_}{y}\omega & \stackrel{\_}{y}\phi & \stackrel{\_}{y}\kappa & \stackrel{\_}{y}s\\ \stackrel{\_}{z}\stackrel{\_}{x}& \stackrel{\_}{z}\stackrel{\_}{y}& \stackrel{\_}{z}& \stackrel{\_}{z}\omega & \stackrel{\_}{z}\phi & \stackrel{\_}{z}\kappa & \stackrel{\_}{z}s\\ \omega \stackrel{\_}{x}& \omega \stackrel{\_}{y}& \omega \stackrel{\_}{z}& \omega & \mathrm{\omega \phi }& \mathrm{\omega \kappa }& \omega s\\ \phi \stackrel{\_}{x}& \phi \stackrel{\_}{y}& \phi \stackrel{\_}{z}& \mathrm{\phi \omega }& \phi & \mathrm{\phi \kappa }& \phi s\\ \kappa \stackrel{\_}{x}& \kappa \stackrel{\_}{y}& \kappa \stackrel{\_}{z}& \mathrm{\kappa \omega }& \mathrm{\kappa \phi }& \kappa & \kappa s\\ s\stackrel{\_}{x}& s\stackrel{\_}{y}& s\stackrel{\_}{z}& s\omega & s\phi & s\kappa & s\end{array}\right]$ where x is translation in x, y is translation in y, where z is translation in z, ω is roll, φ is pitch, κ is yaw, and s is scale.

The bundle adjustment operation 118 can receive the tie points 108, tie point error 110, first 3D point set 104, second 3D point set 106, and the error matrix 116 at input. The bundle adjustment operation 118 can produce a merged 3D point set 120 and a merged 3D point set error 122 as output. The bundle adjustment operation 118 can use a least squares estimator (LSE) for registration of the first 3D point set 104 and the second 3D point set 106. The operation 118 is easily extendable to merging more than two 3D data sets even though the description regards only two 3D data sets at times. The bundle adjustment operation 118 can use one or more photogrammetric techniques. The bundle adjustment operation 118 can include outlier rejection. The bundle adjustment operation 118 can determine error model parameters for the 3D data sets. Application of the error model parameters to the first 3D point set 104 and the second 3D point set 106, results in the relative alignment (registration) of the first 3D point set 104 and the second 3D point set 106.

FIG. 2 illustrates an example diagram of an embodiment of the relationship between ground point coordinate estimates {circumflex over (V)}_j and the corresponding 3D data set observations V_ij. In FIG. 2, three misregistered 3D data sets 202, 204, and 206 and a reference frame 224 are illustrated. First image observations 208, 210, 212 and a first associated ground point 214 and second image observations 216, 218, 220, and a second associated ground point 222 are illustrated. The ground point 214 can be determined using a least squares estimator. The least squares estimator can reduce (e.g., minimize) the discrepancy (across all observations and ground points) across all images. The least squares estimator can propagate an error in one or more of the 3D data sets to an error in a registered 3D data set.

This section establishes some preliminary notational conventions and symbol definitions used for developing the formulations for the bundle adjustment operation 118. The bundle adjustment operation 118 can include identifying a ground point that reduces a discrepancy between the ground point and corresponding points in respective images, and then adjusting points in the 3D data sets in a manner that reduces the discrepancy. The term “3D data set” is sometimes referred to as an “image”. For convenience, example sizes of vectors and matrices are indicated below the symbol. Thus, the symbol

$\hspace{1em}\underset{N\times M}{A}$ denotes a matrix A with N rows and M columns. Column vectors from R³ thus have the annotation 3×1. Components of a vector V are written as

$\underset{3\times 1}{V}={\left[\begin{array}{ccc}x& y& z\end{array}\right]}^T.$ If the vector includes diacritical marks or distinguishing embellishments, these are transferred to the components, as in V=[xyz]^T and V′=[x′ y′ z′]^T.

Equation modeling of the relationship between points in one 3D space to corresponding points in another 3D space is now described. A common reference space is established across all of the images. The reference space can be constructed to accommodate a simultaneous adjustment of more than two images. Correspondences can be formed between points in the reference space and the measured conjugate point locations in each image. The observation equation can be represented as Equation 1:

$\begin{array}{cc}\widetilde{\underset{3\times 1}{V}}=\left(1+s\right)\underset{3\times 3}{T}\left(\underset{3\times 1}{\hat{V}}-\stackrel{\_}{\underset{3\times 1}{V}}\right)& \mathrm{Equation}1\end{array}$

where {circumflex over (V)} is a reference-space 3D coordinate, {tilde over (V)} is the observation of {circumflex over (V)} in an image and the orientation and offset relationship between reference space and image space is taken from the orientation matrix T and offset vector V using Equation 2:

$\begin{array}{cc}T\equiv T\left(\omega ,\varphi ,\kappa \right)=\left[\begin{array}{ccc}c\varphi c\kappa & c\omega s\kappa +s\omega s\varphi c\kappa & s\omega s\kappa -c\omega s\varphi c\kappa \\ -c\varphi s\kappa & c\omega c\kappa -s\omega s\varphi s\kappa & s\omega c\kappa +c\omega s\varphi s\kappa \\ s\varphi & -s\omega c\varphi & c\omega c\varphi \end{array}\right]& \mathrm{Equation}2\end{array}$

where the symbols “c” and “s” denote trigonometric cosine and sine functions, respectively. The quantities

$\underset{3\times 1}{\theta }\equiv {\left[\begin{array}{ccc}\omega & \varphi & \kappa \end{array}\right]}^T$ refer to rotation angles (roll, pitch and yaw) about an image's x, y, and z axes respectively. The scalar s represents an isometric scale correction factor (nominally zero). The above form is conducive to modeling a simultaneous least squares adjustment of all images' offsets and orientations, provided that estimates of reference space coordinates for all conjugate image observations vectors are available. This form is more suitable and flexible than explicitly holding a single image as a reference for at least one of several reasons: (1) there are reference space ground coordinates that permit the potential use of ground control points, whose a priori covariances are relatively small (e.g., they carry high weighting in the solution); (2) the above formulation is suitable for a simultaneous adjustment for data that includes small or minimal overlap (mosaics), as well as, many images collected over the same area (stares) or any combination in between; and (3) allowing all images to adjust simultaneously provides improved geolocation accuracy of the post-adjusted and merged images. parameters.

The symbol {circumflex over (V)} will be referred to as a ground point (akin to tie point ground locations and ground control point locations in a classical photogrammetric image adjustment). The symbol {tilde over (V)} will be referred to as a ground point observation (akin to image tie point observation locations in a classical photogrammetric image adjustment).

Unlike the classical photogrammetric treatment, {circumflex over (V)} and {tilde over (V)} are both “on the ground” in the sense that they both represent ground coordinates in 3D (in the classical imagery case, the observations are in image space and are thus 2D coordinates). Further, the point may very well not be “on the ground” but could be on a building rooftop, treetop canopy, etc. However, the terminology “ground point” and “ground point observation” will be used.

If j is taken to be the index of an arbitrary ground point and i to be the index of an arbitrary image, the observation equation (Equation 1) can be written as Equation 3{tilde over (V)}_ij=(1+s_i)T_i({circumflex over (V)}_j−V_i)  Equation 3

where {circumflex over (V)}_j≡[{circumflex over (x)}_j ŷ_j {circumflex over (z)}_j]^T as the j^th ground point, V_i≡[x_i y_i z_i]T as the offset vector between image i and the reference space origin, and where

$\begin{array}{cc}{T}_i\equiv T\left({\omega }_i,{\varphi }_i,{\kappa }_i\right)=[\begin{array}{ccc}c{\varphi }_{i}c{\kappa }_i& c{\omega }_{i}s{\kappa }_i+s{\omega }_{i}s{\varphi }_{i}c{\kappa }_i& s{\omega }_{i}s{\kappa }_i-c{\omega }_{i}s{\varphi }_{i}c{\kappa }_i\\ -c{\varphi }_{i}s{\kappa }_i& c{\omega }_{i}c{\kappa }_i-s{\omega }_{i}s{\varphi }_{i}s{\kappa }_i& s{\omega }_{i}c{\kappa }_i+c{\omega }_{i}s{\varphi }_{i}s{\kappa }_i\\ s{\varphi }_i& -s{\omega }_{i}c{\varphi }_i& c{\omega }_{i}c{\varphi }_i\end{array}]& \mathrm{Equation}4\end{array}$

is the orientation matrix between image i and the reference space frame and where s_i is image i scale correction factor. Thus, {tilde over (V)}_ij is the coordinate of the i^th image's observation of ground point j.

If a particular ground point is found in two or more images, it can serve as a point which ties the images together (one of the tie points 108). These are generically referred to as tie points. A single tie point is often referred to as a collection of image observations (with coordinates) of the same point on the ground along with the corresponding ground point (with coordinates).

Since the observations over many images i are the measurements containing error, the true ground point coordinates are generally unknown. To facilitate this, an initial estimate of the ground point location can be computed. The initial estimate is provided as Equation 5 as the average of the 3D observations

$\begin{array}{cc}{\widehat{V_j}}^{\left(0\right)}\equiv \frac{1}{I_j}\sum _{i\in I_j}{\tilde{V}}_{\mathrm{ij}}^R& \mathrm{Equation}5\end{array}$

The ground points themselves are treated as derived (but unconstrained) observations and allowed to adjust in performance of the operation 118. There can be an observation of interest whose true ground coordinates are well known. These are classically called ground control points (or GCPs). Since this development can accommodate both GCPs and tie points, the more general terms of “ground point” and “ground point observation” are sometimes used (as contrasted with “tie point ground coordinate” and “tie point observation”).

The bundle adjustment operation 118 can operate on two or more images taken over a same area (with observations for tie points, sometimes called a stare scenario); two or more images taken in strips (forming a mosaic of data, with 2-way, 3-way, or m-way observations in strip overlap regions); tie points in which the corresponding ground points may appear in two or more images, incorporation of GCPs for features in imagery, providing an absolute registration; accommodation of a full covariance for tie point observations. This is conducive for tie point correlation techniques which are highly asymmetrical (e.g., as long as the asymmetry can be characterized as a measurement covariance).

The relationship between ground point coordinate estimates {circumflex over (V)}_j and the corresponding image observations {tilde over (V)}_ij can be understood as a stare scenario between three misregistered images.

For the development of the LSE formulation (and associated preprocessing) that can be performed by the bundle adjustment operation 118, more definitions are provided in Table 1.

TABLE 1
Definitions of Symbols
Symbol                           Definition
V                      R                      W                                                              3                      ×                      1 Location of the origin of the reference frame with respect to the world frame. This is thus the location of the reference frame coordinatized in the world-frame.
V                        _                                            i                                                              3                      ×                      1 Translation of ith image with respect to reference frame origin.         Vi ≡ [xi  yi  zi]T
θ                      i                                                              3                      ×                      1 Orientation angles of ith image with respect to image frame origin         θi ≡ [ωi  ϕi  κi]T
s                      i                                                              1                      ×                      1 Isometric scale factor correction for the ith image.
V                        _                                            i                                              (                        0                        )                                                                                    3                      ×                      1 Initial (zeroth-iteration) value for Vi, thus assumed to be with respect to reference frame origin.
θ                      i                                              (                        0                        )                                                                                    3                      ×                      1 Initial (zeroth-iteration) value for θi. Each element is taken to be zero.
s                      i                                              (                        0                        )                                                                                    1                      ×                      1 Initial (zeroth iteration) value for si. Nominally si(0) ≡ 0.
T                      i                                                              3                      ×                      3 Orientation matrix for ith image built from θi
V                        ^                                            j                                                              3                      ×                      1 Ground point coordinates for ground point j with respect to the reference frame origin         {circumflex over (V)}j ≡ [{circumflex over (x)}j  ŷj  {circumflex over (z)}j]T
V                        ^                                            j                                              (                        0                        )                                                                                    3                      ×                      1 Initial (zeroth-iteration) estimated value for {circumflex over (V)}j
V                        ~                                            ij                      W                                                              3                      ×                      1 Ground point observation coordinate of ground point j on image i coordinatized in the world frame (e.g., these are UTM
                                 coordinates of the ground point observation location).
V                        ~                                            ij                                                              3                      ×                      1 Ground point observation coordinate of ground point j on image i. These are implicitly assumed to be coordinatized in the local image frame for image i.
                                 {tilde over (V)}ij ≡ [{tilde over (x)}ij  {tilde over (y)}ij  {tilde over (z)}ij]T
Σ                        .                                            i                                                              7                      ×                      7 A priori covariance of ith image translation, orientation and scale correction parameter vector
W                        .                                            i                                                              7                      ×                      7 A priori parameter weight matrix for image i.           {dot over (W)}i = {dot over (Σ)}i−1
Σ                        ¨                                            j                                                              3                      ×                      3 A priori covariance of ground point j
W                        ¨                                            j                                                              3                      ×                      3 A priori weight matrix for ground point j.           {umlaut over (W)}j = {umlaut over (Σ)}j−1
Σ                        ~                                            ij                                                              3                      ×                      3 A priori covariance for observation of ground point j upon image i
W                        ~                                            ij                                                              3                      ×                      3 A priori observation weight matrix for observation of ground point j upon image i. {tilde over (W)}ij = ({tilde over (Σ)}ij)−1
General Indexing
m                                Number of images
i                                Image index. i ∈ {1, 2, . . . , m}
n                                Total number of ground points
j                                Ground point index j ∈ {1, 2, . . . , n}
q                                Total number of ground point observations
b                                Ground point observation index b ∈ {1, 2, . . . , q}
(p)                              Non-linear least squares iteration index
Gi                               The index set of all ground points appearing in image i.
                                 Thus Gi ⊆ {1, 2, . . . , n}
Oj                               The index set of observations of ground point j over all images.
                                 Thus Oj ⊆ {1, 2, . . . , q}
Ij                               The index set of images upon which ground point j is an observation.
                                 Thus Ij ⊆ {1, 2, . . . , m}
MbG                              Mapping of observation index to Ground point index. MbG gives the
                                 ground point index (∈ {1, 2, . . . , n) for a specified observation index
                                 b ∈ {1, 2, . . . , q}.
|S|                              Cardinality of set S (e.g., the number of index elements in set S).

Ground point observations can be indexed by ground point j and image i (as in {tilde over (V)}_ij) or by linear indexing, b (as in {tilde over (V)}_b). Use of the subscripting depends upon the context. In the former, it is of interest to characterize the fact that a particular ground point j appears on a particular image i. In the latter, it is of interest to enumerate all observations independent of to which image or to which ground point they refer.

Since some 3D point set data is presented in a “world” space coordinate system (e.g., Universal Transverse Mercator (UTM) map projection) and since the observation Equation 3 is image dependent, some coordinate frame definitions and transformations can aid understanding.

If it is assumed that ground point observation locations are specified in world coordinates, it is of interest to transform the ground point observation locations to be “image” relative. Further, it can be of interest to obtain the ground locations and image offsets themselves to be relative to a “local” reference coordinate frame.

A motivation for a local reference coordinate frame can be to remove large values from the coordinates. For example, UTM coordinates can typically be in the hundreds of thousands of meters. This makes interpretation of the coordinates more difficult, for example, when examining a report of updated coordinate locations. A motivation for an image-relative coordinate frame can be so that the interpretation of the orientation angles comprising the T_i matrices can be relative to the center of the data set. This is contrasted with the origin of rotation being far removed from the data set (e.g., coincident with the local reference frame origin in the mosaic scenario).

In both cases, the transformations between coordinate frames simply involve a 3D translation. The mnemonics W, R and I are used to denote the “world”, “reference” and “image” coordinate frames, respectively. To facilitate the transformations, the following convention is established. A superscript on a vector denotes the coordinate frame to which it is referred. Thus {tilde over (V)}_ij^W corresponds to the world space coordinates of a particular tie point observation, while {tilde over (V)}_ij^R and {tilde over (V)}_ij^I represent the same tie point observation but referred to the reference frame and image frame, respectively.

Following the above convention, the symbol can V_A^B represent “the location of the origin of frame A coordinatized in frame B”. Thus, V_R^V can represent the location of the reference frame in the world coordinate system (e.g., UTM coordinates of the origin of the reference frame). The relationship between an arbitrary vector V^R coordinatized in the reference frame and the same vector V^W coordinatized in the world frame can be represented by Equation 6V^W=V^R+V_R^W  Equation 6

The reference frame can be established as an average of all of the world-space coordinates of tie points. This offset (denoted V_R^W) can be determined using Equation 7

$\begin{array}{cc}{V}_R^W=\frac{1}{q}\sum _{b=1}^q{\tilde{V}}_b^W& \mathrm{Equation}7\end{array}$

For simplicity, it can be assumed that the reference frame origin, referred to by the world frame, can be computed by a process external to the bundle adjustment operation 118 (e.g., by the process that assembles the tie points 108 for use in the bundle adjustment operation 118).

The image frame (e.g., a frame defined on a per-image basis) can be the world coordinates of the center of an image. Under the assumption that there are bounding coordinates in the image data (specifying the min and max extents of the data in world-frame X, Y and Z), the center of the data can thus be taken to be the respective averages of the min and max extents. Since this image frame refers to world space, the computed offset is denoted V_I^W. If bounding coordinates are not available, value for V_I^W is taken as the average of the tie point locations over the specific image i, as described in Equation 8

$\begin{array}{cc}{V}_I^W=\frac{1}{G_i}\sum _{j\in G_i}{\tilde{V}}_{ij}& \mathrm{Equation}8\end{array}$

The image frame offset in reference space coordinates is taken to be the initial value for V⁽⁰⁾ on a per image basis. Thus, for each image i, an external process can compute reference frame coordinates according to Equation 9)V_i^(0)R=V_I^W−V_R^W  Equation 9

Since the tie point observation values can be input in world coordinates and since the observation equation domain assumes reference frame coordinates, some preprocessing of the input data can help make it consistent with that assumed by the observation equation (Equations 1 or 3). The tie point observation coordinates can be converted from world space to reference space. This can be performed for each observation per Equation 10.{tilde over (V)}_ij^R={tilde over (V)}_ij^W−V_R^W  Equation 10

Next, since the true ground point coordinates used in Equation 3 can be unknown, they can be estimated. The ground point coordinates can be assumed to be coordinatized in the reference frame. The initial estimated values for the ground coordinates of each tie point can be computed as an average of the ground point observations over all images in which it appears as described by Equation 11

$\begin{array}{cc}{\widehat{V_j}}^{\left(0\right)}\equiv \frac{1}{I_j}\sum _{i\in I_j}{\tilde{V}}_{\mathrm{ij}}^R& \mathrm{Equation}11\end{array}$

Since the true locations of the tie point ground coordinates can be treated as unknown, the a priori covariance can reflect this by treating the errors in the ground coordinates as numerically unconstrained (in units of meters squared) as described by Equation 12{umlaut over (Σ)}_j≡diag([10¹²10¹²10¹²])  Equation 12

The tie point observation coordinates for use in the observation equation can be converted to image-relative coordinates using Equation 13.{tilde over (V)}_uj={tilde over (V)}_ij^R−V_i⁽⁰⁾  Equation 13

Next, a least squares formulation and solution are discussed. Since the observation equation, Equation 1 or 3, is non-linear in the orientation angles that form T_i, the least squares problem becomes a non-linear least squares problem. Equation 3 can be linearized. Solving the linearized equation can be a multidimensional root finding problem (in which the root is the vector of solution parameters).

For simplification in notation of the linearization, consider a fixed image and a fixed ground point. Let the unknown error model parameters (offset, orientation, and scale correction) be represented by Equation 14:

$\begin{array}{cc}\underset{7\times 1}{X}={\left[\begin{array}{cccccc}\stackrel{\_}{x}& \stackrel{\_}{y}& \stackrel{\_}{z}& \omega & \phi & {\kappa }_S\end{array}\right]}^T& \mathrm{Equation}14\end{array}$

The observation equation for ground point observation V can be written as Equation 15{tilde over (V)}=(1+s)T({circumflex over (V)}−V)  Equation 15

where T is the true image orientation matrix, V is the true image translation vector, {circumflex over (V)} is the true ground point coordinate and {tilde over (V)} is the corresponding ground point observation coordinate.

If one wishes to include the ground point coordinates {circumflex over (V)} as additional observations, the solution for X and {circumflex over (V)} can be cast as a root solving problem based on Equation 16F(X;{circumflex over (V)})=0  Equation 16whereF(X;{circumflex over (V)})={tilde over (V)}−(1+s)T({circumflex over (V)}−V)  Equation 17

In vector form, the function, F, can be represented by Equation 18

$\begin{array}{cc}\left[\begin{array}{c}{f}_1\\ f_2\\ f_3\end{array}\right]=\left[\begin{array}{c}\tilde{x}\\ \tilde{y}\\ \tilde{z}\end{array}\right]-\left(1+s\right)T(\left[\begin{array}{c}\hat{x}\\ \hat{y}\\ \hat{z}\end{array}\right]-\left[\begin{array}{c}\stackrel{\_}{x}\\ \stackrel{\_}{y}\\ \stackrel{\_}{z}\end{array}\right])& \mathrm{Equation}18\end{array}$

The function F can be approximated using a first-order Taylor series expansion of F about initial estimates X⁽⁰⁾ and {circumflex over (V)}⁽⁰⁾ as in Equation 19

$\begin{array}{cc}0=F\left(X;\hat{V}\right)\approx F\left(X^{\left(0\right)};{\hat{V}}^{\left(0\right)}\right)+\frac{\partial F^{\left(0\right)}}{\partial X}\stackrel{.}{\Delta }+\frac{\partial F^{\left(0\right)}}{\partial \hat{V}}\ddot{\Delta }& \mathrm{Equation}19\end{array}$

where X⁽⁰⁾ is an initial approximation of X, {circumflex over (V)}⁽⁰⁾ is an initial approximation of {circumflex over (V)}, the Jacobians

$\frac{\partial F^{\left(0\right)}}{\partial X}\mathrm{and}\frac{\partial F^{\left(0\right)}}{\partial \hat{V}}$ are the partial derivatives of F evaluated at X⁽⁰⁾ and {circumflex over (V)}⁽⁰⁾ respectively, {dot over (Δ)}≡[Δx Δy Δz Δω Δφ Δκ ΔS]^T is a vector of corrections to X, {umlaut over (Δ)}≡[Δ{circumflex over (x)} Δŷ Δ{circumflex over (z)}] is a vector of corrections to {circumflex over (V)}. The values for X⁽⁰⁾ and {circumflex over (V)}⁽⁰⁾ are discussed in Table 1

The Jacobians can be written as Equations 20 and 21

$\begin{array}{cc}\underset{3\times 7}{\stackrel{.}{B}}\equiv \frac{\partial F}{\partial X}=\hspace{1em}\hspace{1em}\left[\begin{array}{c}\frac{\partial f_1}{\partial X}\\ \frac{\partial f_2}{\partial X}\\ \frac{\partial f_3}{\partial X}\end{array}\right]=[\begin{array}{ccccccc}\frac{\partial f_1}{\partial \stackrel{\_}{x}}& \frac{\partial f_1}{\partial \stackrel{\_}{y}}& \frac{\partial f_1}{\partial \stackrel{\_}{z}}& \frac{\partial f_1}{\partial \omega }& \frac{\partial f_1}{\partial \phi }& \frac{\partial f_1}{\partial \kappa }& \frac{\partial f_1}{\partial s}\\ \frac{\partial f_2}{\partial \stackrel{\_}{x}}& \frac{\partial f_2}{\partial \stackrel{\_}{y}}& \frac{\partial f_2}{\partial \stackrel{\_}{z}}& \frac{\partial f_2}{\partial \omega }& \frac{\partial f_2}{\partial \phi }& \frac{\partial f_2}{\partial \kappa }& \frac{\partial f_2}{\partial s}\\ \frac{\partial f_3}{\partial \stackrel{\_}{x}}& \frac{\partial f_3}{\partial \stackrel{\_}{y}}& \frac{\partial f_3}{\partial \stackrel{\_}{z}}& \frac{\partial f_3}{\partial \omega }& \frac{\partial f_3}{\partial \phi }& \frac{\partial f_3}{\partial \kappa }& \frac{\partial f_3}{\partial s}\end{array}]& \mathrm{Equation}20\\ \underset{3\times 3}{\ddot{B}}\equiv \frac{\partial F}{\partial \hat{V}}=\left[\begin{array}{c}\frac{\partial f_1}{\partial \hat{V}}\\ \frac{\partial f_2}{\partial \hat{V}}\\ \frac{\partial f_3}{\partial \hat{V}}\end{array}\right]=\left[\begin{array}{ccc}\frac{\partial f_1}{\partial \hat{x}}& \frac{\partial f_1}{\partial \hat{y}}& \frac{\partial f_1}{\partial \hat{z}}\\ \frac{\partial f_2}{\partial \hat{x}}& \frac{\partial f_2}{\partial \hat{y}}& \frac{\partial f_2}{\partial \hat{z}}\\ \frac{\partial f_3}{\partial \hat{x}}& \frac{\partial f_3}{\partial \hat{y}}& \frac{\partial f_3}{\partial \hat{z}}\end{array}\right]& \mathrm{Equation}21\end{array}$

Note that the dot symbols are merely notations, following the classical photogrammetric equivalent, and do not intrinsically indicate “rates,” as is sometimes denoted in other classical physics contexts.

In matrix notation, Equation 19 can be written as

$\begin{array}{cc}\left[\begin{array}{cc}\frac{\partial F^{\left(0\right)}}{\partial X}& \frac{\partial F^{\left(0\right)}}{\partial \hat{V}}\end{array}\right]\left[\begin{array}{c}\stackrel{.}{\Delta }\\ \ddot{\Delta }\end{array}\right]=-F\left(X^{\left(0\right)};{\hat{V}}^{\left(0\right)}\right)\mathrm{or}& \mathrm{Equation}22\\ \left[\begin{array}{cc}\underset{3\times 7}{\stackrel{.}{B}}& \underset{3\times 3}{\ddot{B}}\end{array}\right]\left[\begin{array}{c}\underset{7\times 1}{\stackrel{.}{\Delta }}\\ \underset{3\times 1}{\ddot{\Delta }}\end{array}\right]=\underset{3\times 1}{\varepsilon }& \mathrm{Equation}23\end{array}$

Since the problem is nonlinear, the estimation of the parameter vector can be iterated (via a multi-dimensional extension of the Newton-Raphson method for root finding, or other technique). The solution can include relinearization at each iteration. The relinearization can be performed about the most recent estimate of the parameter vector. The linearized form of Equation 22 at iteration (p) can be represented as in Equation 24.

$\begin{array}{cc}F\left(X;V\right)\approx F\left(X^{\left(p\right)};{\hat{V}}^{\left(p\right)}\right)+\frac{\partial F^{\left(p\right)}}{\partial X}\stackrel{.}{\Delta }+\frac{\partial F^{\left(p\right)}}{\partial \hat{V}}\ddot{\Delta }& \mathrm{Equation}24\end{array}$

where X^(p) is the p^th iteration estimate of the parameter vector X, {circumflex over (V)}^(p) is the p^th iteration estimate of {circumflex over (V)},

$\frac{\partial F^{\left(p\right)}}{\partial X}$ is the Jacobian of F with respect to X evaluated at X^(p),

$\frac{\partial F^{\left(p\right)}}{\partial \hat{V}}$ is the Jacobian of F with respect to {circumflex over (V)} evaluated at {circumflex over (V)}^(p)), {dot over (Δ)} is a vector of corrections to X for the p^th iteration, and {umlaut over (Δ)} is a vector of corrections to {circumflex over (V)} for the p^th iteration.

With each iteration, the parameter and ground point vectors can be updated with the most recent correction as in Equations 25 and 26.X^(p)=X^(p-1)+{dot over (Δ)}  Equation 25{circumflex over (V)}^(p)={circumflex over (V)}^(p-1)+{umlaut over (Δ)}  Equation 26

For the initial iteration, initial values for X⁽⁰⁾ and {circumflex over (V)}⁽⁰⁾ can be estimated as discussed previously. The system represented by Equation 24 is now linear in {dot over (Δ)} and {umlaut over (Δ)}. A linear solver can be used to solve for the parameters.

For a particular image i and a particular ground point j, Equation 23 can be written as Equation 27

$\begin{array}{cc}\left[\begin{array}{cc}\underset{3\times 7}{{\stackrel{.}{B}}_{\mathrm{ij}}}& \underset{3\times 3}{{\ddot{B}}_{\mathrm{ij}}}\end{array}\right]\left[\begin{array}{c}\underset{7\times 1}{{\stackrel{.}{\Delta }}_i}\\ \underset{3\times 1}{{\ddot{\Delta }}_j}\end{array}\right]=\underset{3\times 1}{{\varepsilon }_{\mathrm{ij}}}& \mathrm{Equation}27\end{array}$

The discrepancy vector for the p^th iteration is thus be represented as in Equation 28ε_ij^(p)=−F(X^(p);{circumflex over (V)}^(p))  Equation 28and thusε_ij^(p)=−[{tilde over (V)}_ij−(1+s_i)T_i({circumflex over (V)}_j−V_i)]  Equation 29

To accommodate a simultaneous solution of all images and ground points, Equation 27 can be extended as

$\begin{array}{cc}\left[\begin{array}{cccccccc}{\stackrel{.}{B}}_{11}& 0& \dots & 0& {\ddot{B}}_{11}& 0& \dots & 0\\ {\stackrel{.}{B}}_{12}& 0& \dots & 0& 0& {\ddot{B}}_{12}& \dots & 0\\ ⋮& ⋮& ⋮& 0& ⋮& ⋮& \ddots & 0\\ 0& {\stackrel{.}{B}}_{21}& \cdots & 0& {\ddot{B}}_{21}& 0& ⋮& 0\\ ⋮& ⋮& ⋮& 0& ⋮& ⋮& ⋮& ⋮\\ 0& 0& \cdots & B_{\mathrm{mn}}& 0& 0& \cdots & {\ddot{B}}_{\mathrm{mn}}\end{array}\right]\left[\begin{array}{c}{\stackrel{.}{\Delta }}_1\\ {\stackrel{.}{\Delta }}_2\\ ⋮\\ {\stackrel{.}{\Delta }}_m\\ {\ddot{\Delta }}_1\\ {\ddot{\Delta }}_2\\ ⋮\\ {\ddot{\Delta }}_n\end{array}\right]=\left[\begin{array}{c}{\varepsilon }_{11}\\ {\varepsilon }_{12}\\ ⋮\\ {\varepsilon }_{21}\\ ⋮\\ {\varepsilon }_{\mathrm{mn}}\end{array}\right]& \mathrm{Equation}30\end{array}$

Equation 30 can be re-written as Equation 31BΔ=E  Equation 31

then the normal equation matrix can be represented as Equation 32 or Equation 33(B^TB)Δ=B^TE  Equation 32ZΔ=H  Equation 33

It can be less efficient to form B as in Equation 30, for one or more of the following reasons: (1) B is very sparse; (2) the quantities {dot over (B)}_ij and {umlaut over (B)}_ij are nonzero if and only if ground point j is observed on image i. For this reason, the classical development of the normal matrix B^TB and right-hand side vector B^TE uses summations over the appropriate indexing. These summations are provided in the normal matrix partitioning below.

The foregoing equations form a foundation for the present problem that is sufficient for development of the normal equations, examination of the normal matrix structure and formulation of the normal equation solution.

The normal equation can be written as in Equation 34ZΔ=H  Equation 34

The matrices can be partitioned as in Equations 35-37

$\begin{array}{cc}Z=\left[\begin{array}{cc}\stackrel{.}{N}+\stackrel{.}{W}& \stackrel{\_}{N}\\ {\stackrel{\_}{N}}^T& \ddot{N}+\ddot{W}\end{array}\right]& \mathrm{Equation}35\\ \Delta =\left[\begin{array}{c}\stackrel{.}{\Delta }\\ \ddot{\Delta }\end{array}\right]& \mathrm{Equation}36\\ H=\left[\begin{array}{c}\stackrel{.}{K}\\ \ddot{K}\end{array}\right]-\left[\begin{array}{c}\stackrel{.}{W}\stackrel{.}{C}\\ \ddot{W}\ddot{C}\end{array}\right]& \mathrm{Equation}37\end{array}$

The quantities {dot over (K)}, {umlaut over (K)}, Ċ and {umlaut over (C)} are described in more details elsewhere herein.

Combining Equations 35, 36 and 37 yields Equation 38

$\begin{array}{cc}Z=\left[\begin{array}{cc}\stackrel{.}{N}+\stackrel{.}{W}& \stackrel{\_}{N}\\ {\stackrel{\_}{N}}^T& \ddot{N}+\ddot{W}\end{array}\right]\left[\begin{array}{c}\stackrel{.}{\Delta }\\ \ddot{\Delta }\end{array}\right]=\left[\begin{array}{c}\stackrel{.}{K}\\ \ddot{K}\end{array}\right]-\left[\begin{array}{c}\stackrel{.}{W}\stackrel{.}{C}\\ \ddot{W}\ddot{C}\end{array}\right]& \mathrm{Equation}38\end{array}$

The matrix Z can thus be represented as Equation 39

$\begin{array}{cc}Z=N+W=\left[\begin{array}{cc}\underset{7m\times 7m}{\stackrel{.}{N}}& \underset{7m\times 3n}{\stackrel{\_}{N}}\\ \underset{3n\times 6m}{{\stackrel{\_}{N}}^T}& \underset{3n\times 3n}{\ddot{N}}\end{array}\right]+\left[\begin{array}{cc}\stackrel{.}{\underset{7m\times 7m}{W}}& 0\\ 0& \ddot{\underset{3n\times 3n}{W}}\end{array}\right]& \mathrm{Equation}39\end{array}$

The matrix {dot over (N)} can be written as Equation 40

$\begin{array}{cc}\stackrel{.}{\underset{7m\times 7m}{N}}=\left[\begin{array}{cccc}\underset{7\times 7}{{\stackrel{.}{N}}_1}& 0& \dots & 0\\ 0& \underset{7\times 7}{{\stackrel{.}{N}}_2}& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & \underset{7\times 7}{{\stackrel{.}{N}}_m}\end{array}\right]& \mathrm{Equation}40\end{array}$

and analogously {dot over (W)} can be written as Equation 41

$\begin{array}{cc}\stackrel{.}{\underset{7m\times 7m}{W}}=\left[\begin{array}{cccc}\underset{7\times 7}{{\stackrel{.}{W}}_1}& 0& \dots & 0\\ 0& \underset{7\times 7}{{\stackrel{.}{W}}_2}& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & \underset{7\times 7}{{\stackrel{.}{W}}_m}\end{array}\right]& \mathrm{Equation}41\end{array}$

The block entries of {dot over (N)}_i can be defined as in Equation 42

$\begin{array}{cc}\underset{7\times 7}{{\stackrel{.}{N}}_i}=\sum _{j\in G_i}{\stackrel{.}{B}}_{\mathrm{ij}}^T{\tilde{W}}_{\mathrm{ij}}{\stackrel{.}{B}}_{\mathrm{ij}}& \mathrm{Equation}42\end{array}$

The subscripts ij on the {dot over (B)}_ij matrices indicate that they are a function of image i and ground point j.

The matrix N can be expanded as in Equation 43

$\begin{array}{cc}\ddot{\underset{3n\times 3n}{N}}=\left[\begin{array}{cccc}\underset{3\times 3}{{\ddot{N}}_1}& 0& \dots & 0\\ 0& \underset{3\times 3}{{\ddot{N}}_2}& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & \underset{3\times 3}{{\ddot{N}}_n}\end{array}\right]& \mathrm{Equation}43\end{array}$

{dot over (W)} can be expanded as in Equation 44:

$\begin{array}{cc}\ddot{\underset{3n\times 3n}{W}}=\left[\begin{array}{cccc}\underset{3\times 3}{{\ddot{W}}_1}& 0& \dots & 0\\ 0& \underset{3\times 3}{{\ddot{W}}_2}& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & \underset{3\times 3}{{\ddot{W}}_n}\end{array}\right]& \mathrm{Equation}44\end{array}$

The block entries of Equation 43 can be defined as in Equation 45

$\begin{array}{cc}\underset{3\times 3}{{\ddot{N}}_i}=\sum _{i\in O_j}{\ddot{B}}_{\mathrm{ij}}^T{\tilde{W}}_{\mathrm{ij}}{\ddot{B}}_{\mathrm{ij}}& \mathrm{Equation}45\end{array}$

The matrix N from Equation 39 can be expanded as in Equation 46

$\begin{array}{cc}\stackrel{\_}{\underset{7m\times 3n}{N}}=\left[\begin{array}{cccc}\underset{7\times 3}{{\stackrel{\_}{N}}_{11}}& \underset{7\times 3}{{\stackrel{\_}{N}}_{12}}& \dots & \underset{7\times 3}{{\stackrel{\_}{N}}_{1n}}\\ \underset{7\times 3}{{\stackrel{\_}{N}}_{21}}& \underset{7\times 3}{{\stackrel{\_}{N}}_{22}}& \dots & \underset{7\times 3}{{\stackrel{\_}{N}}_{2n}}\\ ⋮& ⋮& \ddots & ⋮\\ \underset{7\times 3}{{\stackrel{\_}{N}}_{m1}}& \underset{7\times 3}{{\stackrel{\_}{N}}_{m2}}& \dots & \underset{7\times 3}{{\stackrel{\_}{N}}_{\mathrm{mn}}}\end{array}\right]& \mathrm{Equation}46\end{array}$

The block entries of N from Equation 45 can be defined as in Equation 47

$\begin{array}{cc}\underset{7\times 3}{{\stackrel{\_}{N}}_{\mathrm{ij}}}={\stackrel{.}{B}}_{\mathrm{ij}}^T{\tilde{W}}_{\mathrm{ij}}{\ddot{B}}_{\mathrm{ij}}i\in \left\{1,\dots ,m\right\},j\in \left\{1,\dots ,n\right\}& \mathrm{Equation}47\end{array}$

In a similar development the right hand side matrix H from Equation 34 can be expanded as in Equation 48

$\begin{array}{cc}H=\left[\begin{array}{c}{\stackrel{.}{K}}_1\\ {\stackrel{.}{K}}_2\\ ⋮\\ {\stackrel{.}{K}}_m\\ {\ddot{K}}_1\\ {\ddot{K}}_2\\ ⋮\\ {\ddot{K}}_n\end{array}\right]-\left[\begin{array}{cc}\stackrel{.}{W}& 0\\ 0& \ddot{W}\end{array}\right]\left[\begin{array}{c}\stackrel{.}{C}\\ \ddot{C}\end{array}\right]& \mathrm{Equation}48\end{array}$

The subblocks of H can be defined as in Equations 49 and 50

$\begin{array}{cc}\underset{7\times 1}{{\stackrel{.}{K}}_i}=\sum _{j\in G_i}{\stackrel{.}{B}}_{\mathrm{ij}}^T{\tilde{W}}_{\mathrm{ij}}{\varepsilon }_{\mathrm{ij}}& \mathrm{Equation}49\\ \underset{3\times 1}{{\ddot{K}}_j}=\sum _{i\in O_j}{\ddot{B}}_{\mathrm{ij}}^T{\tilde{W}}_{\mathrm{ij}}{\varepsilon }_{\mathrm{ij}}& \mathrm{Equation}50\end{array}$

with the discrepancy vector ε_ij defined as in Equation 29 and with{umlaut over (C)}^(p)=Ċ^(p-1)−Ċ⁽⁰⁾  Equation 51{umlaut over (C)}^(p)={umlaut over (C)}^(p-1)−{umlaut over (C)}⁽⁰⁾  Equation 52

The values for Ċ⁽⁰⁾ and {umlaut over (C)}⁽⁰⁾ are the initial parameter values. The initial values for the translation parameters portion of Ċ⁽⁰⁾ can be taken to be the V_i^(0)R as computed in Equation 9. The initial values for the rotation parameters portion of Ċ⁽⁰⁾ can be taken to be zero.

The initial values of {umlaut over (C)}⁽⁰⁾ can be taken to be the values of the ground point coordinates {circumflex over (V)}_j⁽⁰⁾ as computed in accord with Equation 11.

The solution to the normal equation matrix on iteration (p) can be determined as in Equation 53Δ^(p)=Z⁻¹H  Equation 53

At each iteration, the parameters can be updated via Equations 51 and Equation 52 and the normal matrix can be formed and solved again. The process can continue until the solution converges. Examples of the convergence criterion can be discussed in the following section.

Since the solution is iterated, a convergence criterion can be established. An example of a convergence criterion is to compute the root-mean-square (RMS) of the residuals as in Equation 54

$\begin{array}{cc}{R}^{\left(p\right)}=\sqrt{\frac{{\varepsilon }^T\varepsilon }{q-7m}}& \mathrm{Equation}54\end{array}$

The value in the denominator of Equation 54 represents the number of degrees of freedom (e.g., the number of observation equations minus the number of estimated parameters).

Since typically q>>7m Equation 54 can be estimated as in Equation 55

$\begin{array}{cc}{R}^{\left(p\right)}=\sqrt{\frac{{\varepsilon }^T\varepsilon }{q}}& \mathrm{Equation}55\end{array}$

The condition q>>7m can be guaranteed with sufficient redundancy of ground point observations as compared with the number of images (e.g., enough tie points are measured between the images so that the aforementioned condition is satisfied).

Convergence happens when the residuals settle to the same values on consecutive iterations. The convergence criterion can be|R^(p)−R^(p-1)|<δ  Equation 56

where δ is a prescribed tolerance.

A rigorous formulation for the standard error of unit weight (to be used in error propagation discussed elsewhere) is provided in Equation 57

$\begin{array}{cc}{\left[{\sigma }^{\left(p\right)}\right]}^2=\frac{\sum _{j=1}^n\sum _{i\in I_j}{\varepsilon }_{ij}^T{\tilde{W}}_{ij}{\varepsilon }_{ij}+\sum _{i=1}^m{\stackrel{.}{C}}_i^T{\stackrel{.}{W}}_i{\stackrel{.}{C}}_i+\sum _{j=1}^n{\ddot{C}}_j^T{\ddot{W}}_j{\ddot{C}}_j}{\mathrm{ndof}}& \mathrm{Equation}57\end{array}$

where ndof is the number of degrees of freedom—the number of observation equations minus the number of error model solution parameters:ndof=q−7m  Equation 58

Since blundered points can be effectively removed from the solution via deweighting, the number of observations remaining effectively doesn't include the blunders. To be strictly correct, the value for q in Equation 58 can be the number of non-blundered observations.

The full form of the matrix Equation 34 can be reduced under the assumption that the errors in the ground point locations are uncorrelated. Under this assumption, the error covariance matrix of the ground point locations {umlaut over (Σ)} becomes a block-diagonal matrix of 3×3 matrix blocks. Since it is a sparse matrix, its inverse is easily computed by inverting the 3×3 diagonal blocks. The development in this section reformulates the normal equations taking advantage of this. The result is a reduced normal equation matrix in which the size of the normal matrix is 6m×6m instead of (6m+3n)×(6m+3n). This gives the obvious advantage that the size of the normal matrix is much smaller and remains invariant with the number of ground points.

The reduced system formation is sometimes referred to as a “ground point folding,” since the ground point portion of the reduced normal matrix is incorporated into the image portion. The development of the reduced normal equation begins with the original normal equation from Equation 34 and repeated as Equation 59ZΔ=H  Equation 59

To facilitate ground point folding into a reduced normal equation matrix, Equation 59 can be re-written as Equation 60

$\begin{array}{cc}\left[\begin{array}{cc}\stackrel{.}{\underset{7m\times 7m}{Z}}& \underset{7m\times 3n}{\stackrel{\_}{Z}}\\ \underset{3n\times 6m}{{\stackrel{\_}{Z}}^T}& \ddot{\underset{3n\times 3n}{Z}}\end{array}\right]\left[\begin{array}{c}\stackrel{.}{\underset{7m\times 1}{\Delta }}\\ \ddot{\underset{3n\times 1}{\Delta }}\end{array}\right]=\left[\begin{array}{c}\stackrel{.}{\underset{7m\times 1}{H}}\\ \ddot{\underset{3n\times 1}{H}}\end{array}\right]& \mathrm{Equation}60\end{array}$

whereŻ={dot over (N)}+{dot over (W)}  Equation 61Z=N  Equation 62{umlaut over (Z)}={umlaut over (N)}+{umlaut over (W)}  Equation 63{dot over (H)}={dot over (K)}−{dot over (W)}Ċ  Equation 64{umlaut over (H)}={umlaut over (K)}−{umlaut over (W)}{umlaut over (C)}  Equation 65

Suppose a matrix system ZΔ=H is partitioned into blocks of the appropriate sizes as

$\begin{array}{cc}\left[\begin{array}{cc}A& B\\ C& D\end{array}\right]\left[\begin{array}{c}a\\ b\end{array}\right]=\left[\begin{array}{c}c\\ d\end{array}\right]& \mathrm{Equation}66\end{array}$

where the matrices A and D are both square.

Further, assume that matrix D is non-singular and can be represented as a sparse block diagonal matrix. Then[A−BD⁻¹C][a]=[c−BD⁻¹d]  Equation 67

Applying Equation 67 to Equation 59 provides the reduced normal matrix equation[Ż−Z{umlaut over (Z)}⁻¹Z^T][{dot over (Δ)}]=[{dot over (H)}−Z{umlaut over (Z)}^−1]  Equation 68

The reduced normal equation matrix can be written as in Equation 69M{dot over (Δ)}=C  Equation 69

where M≡[Ż−Z{umlaut over (Z)}⁻¹Z^T] and C≡[{dot over (H)}−Z{umlaut over (Z)}⁻¹{umlaut over (H)}].

Next it is of interest to examine the form of the components of the reduced system for an efficient implementation. Let {circumflex over (Z)}≡Z{umlaut over (Z)}⁻¹Z^T. Then

$\begin{array}{cc}\widehat{\underset{7m\times 7m}{Z}}=[\begin{array}{ccc}\underset{\left(7\times 3\right)}{{\stackrel{\_}{Z}}_{11}}& \dots & \underset{\left(7\times 3\right)}{{\stackrel{\_}{Z}}_{1n}}\\ ⋮& \ddots & ⋮\\ \underset{\left(7\times 3\right)}{{\stackrel{\_}{Z}}_{m1}}& \dots & \underset{\left(7\times 3\right)}{{\stackrel{\_}{Z}}_{\mathrm{mn}}}\end{array}][\begin{array}{cccc}\underset{\left(3\times 3\right)}{{\ddot{Z}}_1^{-1}}& 0& \dots & 0\\ 0& \underset{\left(3\times 3\right)}{{\ddot{Z}}_2^{-1}}& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & \underset{\left(3\times 3\right)}{{\ddot{Z}}_n^{-1}}\end{array}][\begin{array}{ccc}\underset{\left(3\times 7\right)}{{\stackrel{\_}{Z}}_{11}^T}& \dots & \underset{\left(3\times 7\right)}{{\stackrel{\_}{Z}}_{m1}^T}\\ ⋮& \ddots & ⋮\\ \underset{\left(3\times 7\right)}{{\stackrel{\_}{Z}}_{1n}^T}& \dots & \underset{\left(3\times 7\right)}{{\stackrel{\_}{Z}}_{\mathrm{mn}}^T}\end{array}]& \mathrm{Equation}70\end{array}$

By extension

$\begin{array}{cc}\widehat{\underset{7m\times 7m}{Z}}=[\begin{array}{cccc}\sum _{j=1}^n{\stackrel{\_}{Z}}_{1j}{\ddot{Z}}_j^{-1}{\stackrel{\_}{Z}}_{1j}^T& \sum _{j=1}^n{\stackrel{\_}{Z}}_{1j}{\ddot{Z}}_j^{-1}{\stackrel{\_}{Z}}_{2j}^T& \dots & \sum _{j=1}^n{\stackrel{\_}{Z}}_{1j}{\ddot{Z}}_j^{-1}{\stackrel{\_}{Z}}_{\mathrm{mj}}^T\\ \sum _{j=1}^n{\stackrel{\_}{Z}}_{2j}{\ddot{Z}}_j^{-1}{\stackrel{\_}{Z}}_{1j}^T& \sum _{j=1}^n{\stackrel{\_}{Z}}_{2j}{\ddot{Z}}_j^{-1}{\stackrel{\_}{Z}}_{2j}^T& \dots & \sum _{j=1}^n{\stackrel{\_}{Z}}_{2j}{\ddot{Z}}_j^{-1}{\stackrel{\_}{Z}}_{\mathrm{mj}}^T\\ ⋮& ⋮& \ddots & ⋮\\ \sum _{j=1}^n{\stackrel{\_}{Z}}_{\mathrm{mj}}{\ddot{Z}}_j^{-1}{\stackrel{\_}{Z}}_{1j}^T& \sum _{j=1}^n{\stackrel{\_}{Z}}_{\mathrm{mj}}{\ddot{Z}}_j^{-1}{\stackrel{\_}{Z}}_{2j}^T& \dots & \sum _{j=1}^n{\stackrel{\_}{Z}}_{\mathrm{mj}}{\ddot{Z}}_j^{-1}{\stackrel{\_}{Z}}_{\mathrm{mj}}^T\end{array}]& \mathrm{Equation}71\end{array}$

The blocks of {circumflex over (Z)} in Equation 71 can be the equivalent N_ij as defined in equation 47.

The assumption that errors in the a priori ground points are uncorrelated yields Equation 72

$\begin{array}{cc}\ddot{\underset{3n\times 3n}{W}}=\left[\begin{array}{cccc}{\ddot{\sum _{3\times 3}}}_1^{-1}& 0& \dots & 0\\ 0& {\ddot{\sum _{3\times 3}}}_2^{-1}& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & {\ddot{\sum _{3\times 3}}}_n^{-1}\end{array}\right]& \mathrm{Equation}72\end{array}$

where {umlaut over (Σ)}_j⁻¹ is the inverse of the a priori covariance matrix for ground point j. Thus

$\begin{array}{cc}{\ddot{\underset{3\times 3}{Z}}}_j^{-1}={\left[{\ddot{N}}_j+{\ddot{W}}_j\right]}^{-1}& \mathrm{Equation}73\end{array}$

The general row and column term for {circumflex over (Z)} can then be given by

$\begin{array}{cc}\underset{7\times 7}{{\hat{Z}}_{r,c}}=\sum _{j=1}^n{\stackrel{\_}{Z}}_{\mathrm{rj}}{\ddot{Z}}_j^{-1}{\stackrel{\_}{Z}}_{\mathrm{cj}}^T& \mathrm{Equation}74\end{array}$

and, by the definition of Z_ij, {circumflex over (Z)}_r,c is zero if and only if images r and c have no ground points in common. Also note that {circumflex over (Z)} is block symmetric. Thus, in its formation, only the upper block triangle need be formed, followed by reflection of the upper right triangle to the lower left triangle for completion of the matrix formation.

The matrix M can thus be written as in Equation 75

$\begin{array}{cc}\underset{\left(7m\times 7m\right)}{M}=\left[\begin{array}{cccc}\underset{\left(7\times 7\right)}{{\stackrel{.}{Z}}_1}& 0& \dots & 0\\ 0& \underset{\left(7\times 7\right)}{{\stackrel{.}{Z}}_2}& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & \underset{\left(7\times 7\right)}{{\stackrel{.}{Z}}_p}\end{array}\right]-[\begin{array}{cccc}\sum _{j=1}^n{\stackrel{\_}{Z}}_{1j}{\ddot{Z}}_j^{-1}{\stackrel{\_}{Z}}_{1j}^T& \sum _{j=1}^n{\stackrel{\_}{Z}}_{1j}{\ddot{Z}}_j^{-1}{\stackrel{\_}{Z}}_{2j}^T& \dots & \sum _{j=1}^n{\stackrel{\_}{Z}}_{1j}{\ddot{Z}}_j^{-1}{\stackrel{\_}{Z}}_{\mathrm{mj}}^T\\ \sum _{j=1}^n{\stackrel{\_}{Z}}_{2j}{\ddot{Z}}_j^{-1}{\stackrel{\_}{Z}}_{1j}^T& \sum _{j=1}^n{\stackrel{\_}{Z}}_{2j}{\ddot{Z}}_j^{-1}{\stackrel{\_}{Z}}_{2j}^T& \dots & \sum _{j=1}^n{\stackrel{\_}{Z}}_{2j}{\ddot{Z}}_j^{-1}{\stackrel{\_}{Z}}_{\mathrm{mj}}^T\\ ⋮& ⋮& \ddots & ⋮\\ \sum _{j=1}^n{\stackrel{\_}{Z}}_{\mathrm{mj}}{\ddot{Z}}_j^{-1}{\stackrel{\_}{Z}}_{1j}^T& \sum _{j=1}^n{\stackrel{\_}{Z}}_{\mathrm{mj}}{\ddot{Z}}_j^{-1}{\stackrel{\_}{Z}}_{2j}^T& \dots & \sum _{j=1}^n{\stackrel{\_}{Z}}_{\mathrm{mj}}{\ddot{Z}}_j^{-1}{\stackrel{\_}{Z}}_{\mathrm{mj}}^T\end{array}]& \mathrm{Equation}75\end{array}$

The reduced matrix M can be formed by first storing the diagonal entries of Ż and then subtracting the summed entries of the subtrahend in Equation 75 (namely the {umlaut over (Z)}_r,c defined in Equation 74).

Since the subblocks of the subtrahend are merely summations over the ground point indexes, j, the matrix, M, can be built by iterating over the ground points (assuming the minuend of Equation 75 on-diagonals were formed in advance) and subtracting out the contributions for a particular ground point in the appropriate place within M.

The constant column vector C can be formed similarly with some of the same matrices:

$\begin{array}{cc}C=\left[\begin{array}{c}{\stackrel{.}{H}}_1\\ ⋮\\ {\stackrel{.}{H}}_m\end{array}\right]-\left[\begin{array}{c}\sum _{j=1}^n{\stackrel{\_}{Z}}_{1j}{\ddot{Z}}_j^{-1}{\ddot{H}}_j\\ ⋮\\ \sum _{j=1}^n{\stackrel{\_}{Z}}_{\mathrm{mj}}{\ddot{Z}}_j^{-1}{\ddot{H}}_j\end{array}\right]& \mathrm{Equation}76\end{array}$

After the matrices M and C are built, the solution vector for the adjustable parameters from the reduced system can be computed as{dot over (Δ)}=M⁻¹C  Equation 77

The solution vector can be decomposed into per-image-adjustable vectors {dot over (Δ)}_i for each image i as in Equation 78:

$\begin{array}{cc}\stackrel{.}{\Delta }=\left[\begin{array}{c}{\stackrel{.}{\Delta }}_1\\ {\stackrel{.}{\Delta }}_2\\ ⋮\\ {\stackrel{.}{\Delta }}_m\end{array}\right]& \mathrm{Equation}78\end{array}$

After the solution vector {dot over (Δ)} for the image-adjustable parameters is obtained, the solution vector {umlaut over (Δ)} for corrections to the ground point positions can be extracted (or “unfolded”) from the reduced system. To formulate the extraction, Equation 67 can be used to obtain Equation 79

$\begin{array}{cc}\ddot{\underset{3n\times 1}{\Delta }}={\ddot{Z}}^{-1}\left[\ddot{H}-{\stackrel{\_}{Z}}^T\stackrel{.}{\Delta }\right]\mathrm{If}& \mathrm{Equation}79\\ \ddot{\underset{3n\times 1}{\Delta }}=\left[\begin{array}{c}\underset{3\times 1}{{\ddot{\Delta }}_1}\\ \underset{3\times 1}{{\ddot{\Delta }}_2}\\ ⋮\\ \underset{3\times 1}{{\ddot{\Delta }}_n}\end{array}\right]& \mathrm{Equation}80\end{array}$

represents the correction vector for the ground points then

$\begin{array}{cc}\left[\begin{array}{c}\underset{3\times 1}{{\ddot{\Delta }}_1}\\ \underset{3\times 1}{{\ddot{\Delta }}_2}\\ ⋮\\ \underset{3\times 1}{{\ddot{\Delta }}_n}\end{array}\right]=[\begin{array}{cccc}\underset{\left(3\times 3\right)}{{\ddot{Z}}_1^{-1}}& 0& \dots & 0\\ 0& \underset{\left(3\times 3\right)}{{\ddot{Z}}_2^{-1}}& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & \underset{\left(3\times 3\right)}{{\ddot{Z}}_n^{-1}}\end{array}][[\begin{array}{c}\underset{3\times 1}{{\ddot{H}}_1}\\ \underset{3\times 1}{{\ddot{H}}_2}\\ ⋮\\ \underset{3\times 1}{{\ddot{H}}_n}\end{array}]-\hspace{1em}[\begin{array}{ccc}\underset{\left(3\times 7\right)}{{\stackrel{\_}{Z}}_{11}^T}& \dots & \underset{\left(3\times 7\right)}{{\stackrel{\_}{Z}}_{m1}^T}\\ ⋮& \ddots & ⋮\\ \underset{\left(3\times 7\right)}{{\stackrel{\_}{Z}}_{1n}^T}& \dots & \underset{\left(3\times 7\right)}{{\stackrel{\_}{Z}}_{\mathrm{mn}}^T}\end{array}]\left[\begin{array}{c}\underset{7\times 1}{{\stackrel{.}{\Delta }}_1}\\ \underset{7\times 1}{{\stackrel{.}{\Delta }}_2}\\ ⋮\\ \underset{7\times 1}{{\stackrel{.}{\Delta }}_n}\end{array}\right]]& \mathrm{Equation}81\end{array}$

where {dot over (Δ)}_i is the adjustable parameter correction vector for image i. Thus

$\begin{array}{cc}{\ddot{\underset{3\times 1}{\Delta }}}_j={\ddot{Z}}_j^{-1}\left[{\ddot{H}}_j-\sum _{i\in I_j}{\stackrel{\_}{Z}}_{ij}^T{\stackrel{.}{\Delta }}_i\right]& \mathrm{Equation}82\end{array}$

where I_j is as defined as the index set of images upon which ground point j is an observation.

This section provides formulations for extraction of a posteriori error covariances for ground points. If a priori sensor model error estimates are available (and reliable), the errors may be propagated to the space of the registration error models. In this case, the error propagation is a rigorous predicted error for the accuracy of the a posteriori ground point locations.

The a posteriori error covariances of the image parameters are the appropriate subblocks of the inverse of the reduced normal matrix M⁻¹ from Equation 69 (after application of the variance of unit weight, as described at the end of this section). For the full normal matrix solution, the a posteriori error covariance can be the inverse of the normal matrix, Z⁻¹, times the variance of unit weight. For the reduced system, however, the a posteriori error covariances of the ground points can be extracted from M⁻¹ by unfolding. To facilitate this, the full normal matrix can be written as

$\begin{array}{cc}Z=\left[\begin{array}{cc}\stackrel{.}{\underset{\left(7m\times 7m\right)}{Z}}& \underset{\left(7m\times 3n\right)}{\stackrel{\_}{Z}}\\ \underset{\left(3n\times 7m\right)}{{\stackrel{\_}{Z}}^T}& \underset{\left(3n\times 3n\right)}{\ddot{Z}}\end{array}\right]& \mathrm{Equation}83\end{array}$

Denote the inverse matrix blocks as

$\begin{array}{cc}{Z}^{-1}=\left[\begin{array}{cc}\sum ^{.}& \sum ^{\_}\\ {\sum ^{\_}}^T& \ddot{\sum }\end{array}\right]& \mathrm{Equation}84\end{array}$

Note that, {dot over (Σ)} and {umlaut over (Σ)} as defined are distinctly different from those defined in previous sections. (The symbols in the present section are a posteriori covariances and those in previous sections are a priori covariances). However, this subtle distinction is not problematic if the appropriate context is adhered.

$\begin{array}{cc}\underset{3n\times 3n}{\ddot{\sum }}={\ddot{Z}}^{-1}+{\ddot{Z}}^{-1}{\stackrel{\_}{Z}}^T\stackrel{.}{\Sigma }\stackrel{\_}{Z}{\ddot{Z}}^{-1}& \mathrm{Equation}85\end{array}$

The a posteriori covariance between ground points r and c can be represented as block element

$\underset{3\times 3}{\underset{r,c}{\ddot{\sum }}}\mathrm{of}\ddot{\sum }.$ With n as the number or ground points and m as the number of images,

$\begin{array}{cc}\ddot{\Sigma }=\left[\begin{array}{ccc}\underset{11}{\ddot{\sum }}& \dots & \underset{1n}{\ddot{\sum }}\\ ⋮& \ddots & ⋮\\ \underset{n1}{\ddot{\sum }}& \dots & \underset{\mathrm{nn}}{\ddot{\sum }}\end{array}\right]=\left[\begin{array}{cccc}{\ddot{Z}}_1^{-1}& 0& \dots & 0\\ 0& {\ddot{Z}}_2^{-1}& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & {\ddot{Z}}_n^{-1}\end{array}\right]+\hspace{1em}\hspace{1em}[\begin{array}{cccc}{\ddot{Z}}_1^{-1}& 0& \dots & 0\\ 0& {\ddot{Z}}_2^{-1}& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & {\ddot{Z}}_n^{-1}\end{array}]\hspace{1em}{\stackrel{\_}{Z}}^T\stackrel{.}{\Sigma }\stackrel{\_}{Z}\hspace{1em}[\begin{array}{cccc}{\ddot{Z}}_1^{-1}& 0& \dots & 0\\ 0& {\ddot{Z}}_2^{-1}& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & {\ddot{Z}}_n^{-1}\end{array}]& \mathrm{Equation}86\end{array}$

The r^th row of {umlaut over (Σ)} involves only {umlaut over (Z)}_r⁻¹ of the first {umlaut over (Z)}⁻¹ matrix in term two of Equation 86. Similarly, the c^th column of {umlaut over (Σ)} involves only {umlaut over (Z)}_c⁻¹ of the second {umlaut over (Z)}⁻¹ matrix in term two. Thus{umlaut over (Σ)}_r,c=δ(r,c){umlaut over (Z)}_c⁻¹+{umlaut over (Z)}_r⁻¹Z^T{dot over (Σ)}Z{umlaut over (Z)}_c⁻¹  Equation 87

where the delta function can be

$\begin{array}{cc}\delta \left(r,c\right)\equiv \{\begin{array}{c}1\mathrm{if}r=c\\ 0\mathrm{otherwise}\end{array}& \mathrm{Equation}88\end{array}$

Now the form of the (r, c) block of Z^T {dot over (Σ)}Z is derived. Let

$\begin{array}{cc}\underset{3m\times 3m}{G}\equiv {\stackrel{\_}{Z}}^T\stackrel{.}{\Sigma }\stackrel{\_}{Z}=[\begin{array}{ccc}\underset{\left(3\times 7\right)}{{\stackrel{\_}{Z}}_{11}^T}& \dots & \underset{\left(3\times 7\right)}{{\stackrel{\_}{Z}}_{m1}^T}\\ ⋮& \ddots & ⋮\\ \underset{\left(3\times 7\right)}{{\stackrel{\_}{Z}}_{1n}^T}& \dots & \underset{\left(3\times 7\right)}{{\stackrel{\_}{Z}}_{\mathrm{mn}}^T}\end{array}]\hspace{1em}[\begin{array}{cccc}\underset{\left(7\times 7\right)}{\stackrel{.}{\sum _{11}}}& \underset{\left(7\times 7\right)}{\stackrel{.}{\sum _{12}}}& \dots & \underset{\left(7\times 7\right)}{\stackrel{.}{\sum _{1m}}}\\ \underset{\left(7\times 7\right)}{\stackrel{.}{\sum _{21}}}& \underset{\left(7\times 7\right)}{\stackrel{.}{\sum _{22}}}& \dots & \underset{\left(7\times 7\right)}{\stackrel{.}{\sum _{2m}}}\\ ⋮& ⋮& \ddots & ⋮\\ \underset{\left(7\times 7\right)}{\stackrel{.}{\sum _{m1}}}& \underset{\left(7\times 7\right)}{\stackrel{.}{\sum _{m2}}}& \dots & \underset{\left(7\times 7\right)}{\stackrel{.}{\sum _{mm}}}\end{array}][\begin{array}{ccc}\underset{\left(7\times 3\right)}{{\stackrel{\_}{Z}}_{11}^{}}& \dots & \underset{\left(7\times 3\right)}{{\stackrel{\_}{Z}}_{1n}^{}}\\ ⋮& \ddots & ⋮\\ \underset{\left(7\times 3\right)}{{\stackrel{\_}{Z}}_{m1}^{}}& \dots & \underset{\left(7\times 3\right)}{{\stackrel{\_}{Z}}_{\mathrm{mn}}^{}}\end{array}]& \mathrm{Equation}89\end{array}$

The r^th row of G involves only the r^th row of Z^T and the c^th column of G involves only the c^th column of Z. Thus

$\begin{array}{cc}\underset{3\times 3}{G_{r,c}}=\left[\begin{array}{cc}{\overline{Z}}_{1r}^T\stackrel{.}{\underset{11}{\Sigma }}+\dots +{\overline{Z}}_{mr}\underset{m1}{\stackrel{.}{\Sigma }}& \begin{array}{cc}{\overline{Z}}_{1r}^T\stackrel{.}{\underset{12}{\Sigma }}+\dots +{\overline{Z}}_{mr}\underset{m2}{\stackrel{.}{\Sigma }}& \dots \end{array}\end{array}\right]\stackrel{\_}{Z}=\left[\begin{array}{ccc}\underset{s=1}{\sum ^m}{\overline{Z}}_{\mathrm{sr}}^T\stackrel{.}{\sum _{s1}}& \underset{s=1}{\sum ^m}{\overline{Z}}_{\mathrm{sr}}^T\stackrel{.}{\sum _{s2}}& \dots \end{array}\right]\left[\begin{array}{c}{\overline{Z}}_{1c}\\ {\overline{Z}}_{2c}\\ ⋮\\ {\overline{Z}}_{mc}\end{array}\right]=\sum _{t=1}^m\left[\sum _{s=1}^m{\overline{Z}}_{\mathrm{sr}}^T\stackrel{.}{\sum _{st}}\right]{\overline{Z}}_{tc}& \mathrm{Equation}90\end{array}$

Now Z_ij=0 if ground point j is not an observation on image i.

ThusT_r,c=Σ_t∈Ic[Σ_s∈IrZ_sr^T{dot over (Σ)}_st]Z_tc  Equation 91

where I_j is the index set of images upon which ground point j is an observation. Substituting Equation 91 into Equation 87 yields{umlaut over (Σ)}_r,c=δ(r,c){umlaut over (Z)}_c⁻¹+{umlaut over (Z)}_r⁻¹[Σ_t∈Ic[Σ_s∈IrZ_sr^T{dot over (Σ)}_st]Z_tc]{umlaut over (Z)}x₂⁻¹  Equation 92

The a posteriori covariance is usually defined by scaling the inverse of the normal matrix by an estimate of the variance of unit weight. An estimate of the variance of unit weight is denoted as [σ^(p)]² and is provided in Equation 57. Thus, the above formulation can be used, but instead defining

$\begin{array}{cc}\left[\begin{array}{cc}\sum ^{.}& \sum ^{\_}\\ {\sum ^{\_}}^T& \ddot{\sum }\end{array}\right]={\left[{\sigma }^{\left(p\right)}\right]}^2\left[Z^{-1}\right]& \mathrm{Equation}93\end{array}$

For a full normal matrix solution, Z⁻¹ is readily available, thus the a posteriori covariance of the error model parameters and ground points can be the right hand side of Equation 93.

The right hand summand of Equation 92 includes the factor [σ^(p))]² since it includes {dot over (Σ)}_st. However, the left hand summand does not include the factor. This can be compensated for by a modified form of Equation 92{umlaut over (Σ)}_r,c=[σ^(p)]²δ(r,c){umlaut over (Z)}_c⁻¹+{umlaut over (Z)}_r⁻¹[Σ_t∈Ic[Σ_s∈IrZ_sr^T{dot over (Σ)}_st]Z_tc]{umlaut over (Z)}_c⁻¹  Equation 94

If the standard error of unit weight σ^(p) is deemed to be unreliable (e.g., is much greater than unity) this may be an indicator of improper (or incorrect) a priori error covariance in the process. One can still, however, be able to provide a reliable error estimate from the least squares process by simply forcing the standard error to one (e.g., by setting σ^(p)←1.0 in Equations 93 and 94.

FIG. 3 illustrates an example of an embodiment of the operation 118. The operation 118 can include 3D data set registration with error propagation. The operation 118, as illustrated, includes initializing solution and corrections, at operation 302; computing discrepancies, at operation 304; computing a normal equation solution, at operation 306; updating parameters based on the computed normal equation solution, at operation 308; computing residuals, at operation 310; extracting propagated error, at operation 312; and compensating misregistration of the first 3D point set 104 and the second 3D point set 106, at operation 314.

The operation 302 can include setting the solution vector X and the correction vector ΔX to the zero vector¹: ¹ The solution vector X is set to the fixed-point location for the linearization. If an a priori estimate is available, it is used here in place of the zero vector.

$X=\underset{7\times 1}{0}$ $\Delta X=\underset{7\times 1}{0}$

The operation 304 can include computing the discrepancy vector for each observation as provided in Equation 29. The operation 306 can include building the normal equations matrices and solving for the correction vector as provided in Equation 53. The operation 308 can include updating the parameter vector for the current iteration as provided in Equations 51 and 52. Details of the operation 308 for unfolding of the ground points for the folded normal equation solution is provided via pseudocode below.

The operation 310 can include computing the residuals (final discrepancies) as provided in Equation 29. The operation 312 can include computing a standard error of unit weight via Equation 57. Note that the standard error can the square root of the left hand side of the Equation 57 (e.g., σ^(p)=√{square root over ([σ^(p)]²))}.

If the delta between the current and previous standard error of unit weight is less than the convergence criterion in absolute value, the solution nominally converged. To accommodate blunder rejection, the convergence criterion check can be augmented with a check to see if the blunder weights should be used in continuation of the solution (“useBW”, indicating to use “blunder-checking weighting”). If convergence occurs and useBW is true, this is an indicator to perform blunder checking, and this time using a normalized residual computation in order to check for blunders on the next iteration.

If useBW is true, blunders can be computed. If there are blunders remaining, the blunder “cycle” number is incremented and the process is repeated with the correction vector reset to a priori values (e.g., go to operation 302). If there are no blunders remaining, a check can be performed to see if the number of post convergence blunder cycles can be set to zero. This check can be performed to effectively force one more solution after all blunders have been eliminated.

If useBW is false and it is currently a first iteration of a blundering solution, useBW can be set to true. This has the effect of forcing the normalized residual blunder iteration for determining the blunders on subsequent iterations. In this case, a solution has converged but normalized blunder residuals have not been computed. Setting useBW to true can forces this to happen on the next solution iteration. The solution can be iterated by going to the operation 306. If there are no more blunders and the number of blunders is not zero, this indicates the “non-blunder iteration” solution has converged.

The operation 118 can include providing a report that includes an iteration number, current correction vector ΔX, current iteration estimates of parameters and ground points (e.g., as computed in equations 51 and 52), standard error of unit weight (e.g., as provided in Equation 55). The operation 118 can include a check for non-convergence by examining the current iteration number with a maximum number of iterations, M. If the number of iterations exceeds the maximum, stop the iteration process. The solution did not converge. An exception can be raised and the operation 118 can be complete.

The following is a pseudocode outline for the operation 118 for computing the full normal equations solution. This first pseudocode does not include ground point folding.

The names in {braces} allude to method (e.g., function) names in a software implementation. Also, within the development below, ground point indexes, ground point observation indexes and image indexes are assumed to be zero-relative. For efficiency in the implementation, the following elements can be cached in a per-image workspace object, which is updated with each iteration:

• • 1. Values for trig functions sω, cω, sϕ, cϕ, sκ, cκ as given in Equation 4.
  • 2. Three 3×3 partial derivative matrices of T with respect to the three angles as given in the following Equations 95-97:

$\begin{array}{cc}{T}^{\omega }\equiv \frac{\partial T}{\partial \omega }=\left[\begin{array}{ccc}0& -s\omega s\kappa +c\omega s\varphi c\kappa & c\omega s\kappa +s\omega s\varphi c\kappa \\ 0& -s\omega c\kappa -c\omega s\varphi s\kappa & c\omega c\kappa -s\omega s\varphi s\kappa \\ 0& -c\omega c\varphi & -s\omega c\varphi \end{array}\right]& \mathrm{Equation}95\\ T^{\varphi }\equiv \frac{\partial T}{\partial \varphi }=\left[\begin{array}{ccc}-s\varphi c\kappa & s\omega c\varphi c\kappa & -c\omega c\varphi c\kappa \\ s\varphi s\kappa & -s\omega c\varphi s\kappa & c\omega c\varphi s\kappa \\ c\varphi & s\omega s\varphi & -c\omega s\varphi \end{array}\right]& \mathrm{Equation}96\\ T^{\kappa }\equiv \frac{\partial T}{\partial \kappa }=\left[\begin{array}{ccc}-c\mathrm{\varphi s\kappa }& c\omega c\kappa -s\omega s\varphi s\kappa & s\omega c\kappa +c\omega s\varphi s\kappa \\ -c\varphi c\kappa & -c\omega s\kappa -s\omega s\mathrm{\varphi c\kappa }& -s\omega s\kappa +c\omega s\varphi c\kappa \\ 0& 0& 0\end{array}\right]& \mathrm{Equation}97\end{array}$

• • 3. Rotation matrix T_i as given in Equation 4
  • 4. 3×1 vector V_i

The pseudocode begins by setting the non-linear least squares iteration index (p) to zero.

1.  {initializeData} Store nominal initial values for the solved
    parameters for each image i ∈ {0, . . . , m − 1}. The initial values for Vi(0) and
    θi(0).
Equation 98
C                        .                                            i                                              (                        0                        )                                                              ≡                                                                  [                                                                                                                                                                              V                                  _                                                                i                                                                  (                                  0                                  )                                                                                                                                                                                                                                        θ                                i                                                                  (                                  0                                  )                                                                                                                                                                                                                                        s                                i                                                                  (                                  0                                  )                                                                                                                                                                    ]                                                                    7                        ×                        1
2.  {initializeData} Compute
    a. Initial ground point coordinates via Equation 11 for each ground
       point j ∈ {0, . . . , n − 1}. These form {umlaut over (C)}j(0).
    b. Initial image cache data as described previously.
3.  {outputIterationReport} Output initial iteration report (input
    parameters, initial ground point and ground point observation coordinates.)
4.  {initializeMatrices} Block partition the {dot over (N)} portion of Z into m
    subblocks, each of size 6 × 6.
5.  {initializeMatrices} Block partition the {umlaut over (N)} portion of Z into n
    subblocks, each of size 3 × 3.
6.  {initializeMatrices} Block partition H similarly
7.  {computeAndStoreDiscrepancies} For each ground point index j ∈
    {0, . . . , n − 1} and each observation of that ground point i ∈ Oj, compute the
    discrepancy vector εij given in Equation 29 as follows.
    for j ∈ {0, . . . , n − 1}
    {
       Fetch the most recent ground point position {circumflex over (V)}j.
       for i ∈ Ij and observation index b ∈ Oj
       {
       a) Retrieve image-cached values for Ti and Vi
       b) Retrieve ground point observation {tilde over (V)}ij = {tilde over (V)}b
       c) Apply the observation equation to obtain the projected value for {circumflex over (V)}j.
Equation 99
Vij = (1 + si)Ti({circumflex over (V)}j − Vi)
       d) Compute and store the discrepancy vector for observation b as in
          Equation 29
Equation 100
εb ≡ εij = −[{tilde over (V)}ij − Vij]
       }// end for i
    }//end for j
8.  Compute the standard error of unit weight as in Equation 57.
9.  {buildNormalMatrices + initialize Weights} Zero Z and H
    and initialize Z (and likewise H) with the block weight matrices on the
    diagonal. This involves setting the blocks of Z to the subblocks of {dot over (W)} and {umlaut over (W)}.
    and setting the subblocks (subrows) of H to −{dot over (W)}Ċ and −{umlaut over (W)}{umlaut over (C)}.
10. {sumInPartials} Loop over ground points and images containing the
    ground points and sum in the contributions of the {dot over (B)} and {umlaut over (B)} matrices into Z
    and H.
    for j ∈ {0, . . . , n − 1}
    {
       for i ∈ Ij
       {
          Retrieve εb ≡ εij as computed in Equation 100
          Compute {dot over (B)}ij and {umlaut over (B)}ij as in Equations 101 and 102.
Equation 101
B                          .                                                ij                                                                    3                        ×                        7                                                              =                                          [                                                                                                                                  (                                                              1                                +                                                                  s                                  i                                                                                            )                                                        ⁢                                                          T                              i                                                                                                            3                            ×                            3                                                                          ⁢                                                                                                  ⁢                                                                                                            (                                                              1                                +                                                                  s                                  i                                                                                            )                                                        ⁢                                                          A                              ij                                                                                                            3                            ×                            3                                                                          ⁢                                                                                                  ⁢                                                                                                            -                                                              T                                i                                                                                      ⁢                                                          Y                              ij                                                                                                            3                            ×                            1                                                                                              ]
Equation 102
B                          ¨                                                ij                                                                    3                        ×                        3                                                              =                                                                  -                                                  (                                                      1                            +                                                          s                              i                                                                                )                                                                    ⁢                                                                                          ⁢                                                                        [                                                      T                            i                                                    ]                                                                          3                          ×                          3
          Retrieve observation weight matrix {tilde over (W)}ij
          {dot over (N)}i: Sum {dot over (B)}ijT{tilde over (W)}ij{dot over (B)}ij into Z.block(i, i)
          {umlaut over (N)}j: Sum {umlaut over (B)}ijT{tilde over (W)}ij{umlaut over (B)}ij into Z.block(m + j, m + j)
          N              ij: Sum {dot over (B)}ijT{tilde over (W)}ij{umlaut over (B)}ij into Z.block(i, m + j)
          Ċi: Sum {dot over (B)}ijT{tilde over (W)}ijεij into H.block(i, 0)
          {umlaut over (C)}j: Sum {umlaut over (B)}ijT{tilde over (W)}ijεij into H.block(m + j, 0)
       } //end i
    } // end j
11. {solveNormalEquation} Form the lower transpose of the Z matrix and
    solve the system Δ = Z−1H. Note that the normal equation system is a
    symmetric system (e.g., the normal matrix Z is symmetric). Thus, a
    symmetric system solver can be used. In the case of a symmetric system
    solver, it may not be necessary to form the lower triangle.
12. {updateParameters} Update all the parameters by extracting the
    corrections from the Δ matrix as in Equations 51 and 52.
13. If (p) ≠ 0 compare with the previous RMS of residuals and check for
    convergence. The convergence condition can be
Equation 103
|R(p) − R(p−1)| < δ
14. {computePostConvergenceResiduals +
    checkForBlunders} If convergence has been reached, perform
    automatic blunder editing. If convergence has not been reached, increment
    the iteration index
Equation 104
(p) ← (p + 1)
    and go to step 7.

What follows is pseudocode for the operation 118 for building the reduced normal equations system, computing corrections to ground point positions and performing error propagation via extraction of data from the reduced normal matrix. This portion of the pseudocode includes ground point coordinate folding.

As in the full normal solution provided in the previous pseudocode, the same per-image elements are cached in a workspace object and updated with each iteration. The algorithm for the reduced solution can be broken into two major portions: priming and folding. Priming involves storing of weights and the contributions along the diagonal of the full normal equation matrix (and corresponding data for the right hand column vector H). This corresponds to the Ż portion of Z. Thus, priming involved formation of the minuends of Equation 75 and Equation 76. Folding can include incorporation of the subtrahends of the aforementioned Equations.

To provide an efficient implementation, a ground point workspace can be created. The workspace can include the following elements: {umlaut over (Z)}_j, {umlaut over (H)}_j, {umlaut over (Z)}_j⁻¹. These things are indexed by ground point for the ground point workspace. The technique can begin by setting the non-linear least squares iteration index (p) to zero.

1.  {initializeData} Store nominal initial values for the solved
    parameters for each image i ∈ {0, ... , m − 1}. The initial values for Vi(0) and
    θi(0) can also be set along with Ċi(0) as in Equation 98.
2.  {initializeData} Compute
    a. Initial ground point coordinates via Equation 11 for each ground
    point j ∈ {0, ... , n − 1}. These form {umlaut over (C)}j(0).
    b. Initial image cache data as described above.
3.  {outputIterationReport} Output initial iteration report (input
    parameters, initial ground point, and ground point observation coordinates).
4.  {initializeMatrices} Block partition the reduced normal matrix
    Minto m subblocks, each of size 6 × 6. Block partition the reduced column
    vector C similarly
5.  {computeAndStoreDiscrepancies} For each ground point index j ∈
    {0, ... , n − 1} and each observation of that ground point i ∈ Oj, compute the
    discrepancy vector εij given in Equation 29 as:
    for j ∈ {0, ... , n − 1}
    {
    Fetch the most recent ground point position {circumflex over (V)}j.
    for i ∈ Ij and observation index b ∈ Oj
    {
    a) Retrieve image-cached values for Ti and Vi
    b) Retrieve ground point observation {tilde over (V)}ij = {tilde over (V)}b
    c) Apply the observation equation to obtain the projected value for {circumflex over (V)}j
    from Equation 99.
    d) Compute and store the discrepancy vector for observation b as in
    Equation 100
    }// end for i
    }//end for j
6.  Compute the standard error of unit weight.as in Equation 57.
7.  {buildNormalMatrices}
    a. Zero the reduced matrices M and C
    b. {initialize Weights} Initialize M (and likewise C) with the
    block weight matrices on the diagonal. This involves setting the
    blocks of M to the subblocks of {dot over (W)} and setting the subblocks
    (subrows) of C to −{dot over (W)} Ċ.
    c. {sumInPartialsAndFoldGPs} Form the main diagonal and
    ground point matrices {umlaut over (Z)}j by iterating over ground points. Perform
    folding for ground point j
    for j ∈ {0, ... , n − 1}
    {
    PRIMING:
    Store {umlaut over (W)}j into {umlaut over (Z)}j of GPWS
    Store −{umlaut over (W)}j{umlaut over (C)}j into {umlaut over (H)}j of GPWS
    for i ∈ Ij (where Ij is set of image indexes upon which GP j is an
    observation)
    {
    Form partial derivatives:
    Build {dot over (B)}ij as in Equation 101
    Build {umlaut over (B)}ij as in Equation 102
    Retrieve discrepancy vector εij as computed in Equation 100.
    Retrieve observation weight matrix {tilde over (W)}ij
    Sum in contribution of GP j's obs in image i within M:
    Sum {dot over (B)}ijT{tilde over (W)}ij{dot over (B)}ij into M.block(i,i)
    Sum {dot over (B)}ijT{tilde over (W)}ijεij into C.block(i,0)
    Sum in i's contribution to {umlaut over (Z)}j and {umlaut over (H)}j:
    Sum {umlaut over (B)}ijT{tilde over (W)}ij{umlaut over (B)}ij into {umlaut over (Z)}j
    Sum {umlaut over (B)}ijT{tilde over (W)} εij into {umlaut over (H)}j
    }// end i
    Invert {umlaut over (Z)}j and store into GPWS as {umlaut over (Z)}j−1
    FOLDING INTO M (note : iteration loop over j is still valid)
    for r ∈ Ij
    {
    Form Zrj = {dot over (B)}ijT{tilde over (W)}ij{umlaut over (B)}ij as in Equations 69 and 47
    for c ∈ Ij|c ≥ r
    {
    Form ZcjT
    Sum in −Zrj{umlaut over (Z)}j−1ZcjT into M.block(r, c).
    }//end c
    Sum in −Zrj{umlaut over (Z)}j−1{umlaut over (H)}j into C.block(r, 0).
    }//end r
    }// end j
8.  Complete the lower diagonal entries of M and solve {dot over (Δ)} = M−1C. As in the
    full normal equation solution, note that M is symmetric and thus a symmetric
    system solver is in order.
9.  First use pseudocode provided below to compute corrections to ground
    points. Then update all the parameters from the {dot over (Δ)} vector.
10. If (p) ≠ 0 compare with the previous RMS of residuals and check for
    convergence. The convergence condition is
    |R(p) − R(p − 1)| < ε Equation 105
11. {computePostConvergenceResiduals + checkForBlunders} If convergence
    has been reached, perform automatic blunder editing as detailed elsewhere.
    After there are no more blunders, proceed with error propagation
    {propagateErrors}. If convergence has not been reached, increment the
    iteration index as in Equation 104 and go to step 5.

After the solution vector {dot over (Δ)} is obtained, unfolding the ground point corrections is a matter of employing Equation 80, replicated here for reference:

Equation 80
Δ                          ¨                                                j                                                                    3                        ×                        1                                                              =                                                                                            Z                          ¨                                                j                                                  -                          1                                                                    ⁡                                              [                                                                                                            H                              ¨                                                        j                                                    -                                                                                    Σ                                                              i                                ∈                                                                  I                                  j                                                                                                                      ⁢                                                                                          Z                                _                                                            ij                              T                                                        ⁢                                                                                          Δ                                .                                                            i                                                                                                      ]
{unfoldGroundPoints}
for j ∈ {0, . . . , n − 1}
{
Retrieve {umlaut over (Z)}j−1 and {umlaut over (H)}j from GPWS
Store {umlaut over (H)}j into a new subtrahend matrix S (i.e. initialize S to {umlaut over (H)}j)
for r ∈ Ij
{
Form Zrj
Sum −ZrjT{dot over (Δ)}r into S.
}// end r
Compute {umlaut over (Δ)}j = {umlaut over (Z)}j−1S
}//end j
The general cross error covariance between ground point indexes r
and c can obtained by evaluation of Equation 94.
{relativeGroundPointCov}
Retrieve {umlaut over (Z)}r−1 and {umlaut over (Z)}c−1 from ground point workspace.
Obtain indexing sets Ir and Ic (image indexes of ground points r and c)
Allocate                    ⁢                                                                                  ⁢                    matrix                    ⁢                                                                                  ⁢                                                                                            Σ                          ¨                                                                          r                          ,                          c                                                                                            3                        ×                        3                                                              ⁢                                                                                  ⁢                    and                    ⁢                                                                                  ⁢                    initialize                    ⁢                                                                                  ⁢                    to                    ⁢                                                                                  ⁢                                          zero                      .                                                                                          ⁢                      This                                        ⁢                                                                                  ⁢                    is                    ⁢                                                                                  ⁢                    the                    ⁢                                                                                  ⁢                    output                    ⁢                                                                                  ⁢                    of                    ⁢                                                                                  ⁢                    this                    ⁢                                                                                  ⁢                                          function                      .
Allocate                    ⁢                                                                                  ⁢                    matrix                    ⁢                                                                                  ⁢                                          P                                              3                        ×                        7                                                              ⁢                                                                                  ⁢                    and                    ⁢                                                                                  ⁢                    initialize                    ⁢                                                                                  ⁢                    to                    ⁢                                                                                  ⁢                    zero
for t ∈ Ic
{
Allocate⁢⁢matrix⁢⁢Q3×7⁢⁢and⁢⁢initialize⁢⁢to⁢⁢zero
for s ∈ Ir
{
Form Zsr and Ztc from Equation 47. (Note that for the reduced system, N = Z)
Extract {dot over (Σ)}st = {dot over (Σ)}.getBlock(s, t), where {dot over (Σ)} is as defined in Equation 93
Set Q ← Q + ZstT{dot over (Σ)}st
}//end s
Set P ← P + QZrc
}//end t
Compute {umlaut over (Σ)}r,c = {umlaut over (Z)}r−1P{umlaut over (Z)}c−1
if(r = c)
{
Set {umlaut over (Σ)}r,c ← [σ(p)]2{umlaut over (Z)}c−1 + {umlaut over (Σ)}r,c
}
return {umlaut over (Σ)}r,c
//end {relativeGroundPointCov}

The full 3n×3n ground error covariance matrix

$\underset{3n\times 3n}{\ddot{\sum }}$ may be obtained by invoking the method for r∈{1, 2, . . . , n} and for c∈{r, r+1, . . . , n}. Note that the indexing for c starts with r since the full ground covariance matrix is symmetric (i.e., build the upper triangle of

$\underset{3n\times 3n}{\ddot{\sum }}$ and “reflect about the diagonal” to obtain the lower symmetric portion).

What follows regards how to perform operation 314. The operation 314 proceeds given the outputs of the LSE techniques discussed. The compensation applies the inverse of the observation equation, accommodating the various relative frame offsets to arrive at compensated world space coordinates from misregistered world space coordinatesV_reg^W=Compensate_i(V_misreg^W)  Equation 106

The motivation for providing the inputs and outputs in world space coordinates can be that is the native space of the inputs and desired space of the outputs for each element of each image's point cloud.

For an arbitrary misregistered vector V_misreg^W on image i, the compensation formula can be performed as in Equation 107

$\begin{array}{cc}{V}_{reg}^W=\frac{1}{\left(1+s_1^{\left(p\right)}\right)}{T}_i^T\left(V_{\mathrm{misreg}}^W-{\overline{V_i}}^{{\left(0\right)}^R}-V_R^W\right)+{\overline{V_i}}^{\left(p\right)}+V_R^W& \mathrm{Equation}107\end{array}$

where T_i is constructed from the solution vector θ_i^(p) and the other symbols in Equation 107 are defined elsewhere. Note that the values for −V^(0)R−V_R^W and +V_i^(p)+V_R^W and T_i^T can be precomputed on a per image basis when applying Equation 107 for a time-efficient implementation.

FIG. 4 illustrates, by way of example, a block diagram of an embodiment of a machine in the example form of a computer system 400 within which instructions, for causing the machine to perform any one or more of the methods discussed herein, may be executed. In a networked deployment, the machine may operate in the capacity of a server or a client machine in server-client network environment, or as a peer machine in a peer-to-peer (or distributed) network environment. The machine may be a personal computer (PC), a tablet PC, a set-top box (STB), a Personal Digital Assistant (PDA), a cellular telephone, a web appliance, a network router, switch or bridge, or any machine capable of executing instructions (sequential or otherwise) that specify actions to be taken by that machine. Further, while only a single machine is illustrated, the term “machine” shall also be taken to include any collection of machines that individually or jointly execute a set (or multiple sets) of instructions to perform any one or more of the methodologies discussed herein.

The example computer system 400 includes a processor 402 (e.g., a central processing unit (CPU), a graphics processing unit (GPU) or both), a main memory 404 and a static memory 406, which communicate with each other via a bus 408. The computer system 400 may further include a video display unit 410 (e.g., a liquid crystal display (LCD) or a cathode ray tube (CRT)). The computer system 400 also includes an alphanumeric input device 412 (e.g., a keyboard), a user interface (UI) navigation device 414 (e.g., a mouse), a mass storage unit 416, a signal generation device 418 (e.g., a speaker), a network interface device 420, and a radio 430 such as Bluetooth, WWAN, WLAN, and NFC, permitting the application of security controls on such protocols.

The mass storage unit 416 includes a machine-readable medium 422 on which is stored one or more sets of instructions and data structures (e.g., software) 424 embodying or utilized by any one or more of the methodologies or functions described herein. The instructions 424 may also reside, completely or at least partially, within the main memory 404 and/or within the processor 402 during execution thereof by the computer system 400, the main memory 404 and the processor 402 also constituting machine-readable media.

While the machine-readable medium 422 is shown in an example embodiment to be a single medium, the term “machine-readable medium” may include a single medium or multiple media (e.g., a centralized or distributed database, and/or associated caches and servers) that store the one or more instructions or data structures. The term “machine-readable medium” shall also be taken to include any tangible medium that is capable of storing, encoding or carrying instructions for execution by the machine and that cause the machine to perform any one or more of the methodologies of the present invention, or that is capable of storing, encoding or carrying data structures utilized by or associated with such instructions. The term “machine-readable medium” shall accordingly be taken to include, but not be limited to, solid-state memories, and optical and magnetic media. Specific examples of machine-readable media include non-volatile memory, including by way of example semiconductor memory devices, e.g., Erasable Programmable Read-Only Memory (EPROM), Electrically Erasable Programmable Read-Only Memory (EEPROM), and flash memory devices; magnetic disks such as internal hard disks and removable disks; magneto-optical disks; and CD-ROM and DVD-ROM disks.

The instructions 424 may further be transmitted or received over a communications network 426 using a transmission medium. The instructions 424 may be transmitted using the network interface device 420 and any one of a number of well-known transfer protocols (e.g., HTTP). Examples of communication networks include a local area network (“LAN”), a wide area network (“WAN”), the Internet, mobile telephone networks, Plain Old Telephone (POTS) networks, and wireless data networks (e.g., WiFi and WiMax networks). The term “transmission medium” shall be taken to include any intangible medium that is capable of storing, encoding or carrying instructions for execution by the machine, and includes digital or analog communications signals or other intangible media to facilitate communication of such software.

ADDITIONAL NOTES AND EXAMPLES

The subject matter can be understood by way of some examples.

Example 1 can include a method for registering a first three-dimensional (3D) image to a second 3D image with error propagation, the method comprising reducing a sum aggregate of discrepancies between respective tie points and associated 3D points in the first and the second 3D images; adjusting 3D error models of the first and second 3D images based on the reduced discrepancies to generate registered 3D images, and propagating an error of the first or second 3D images to the registered 3D image to generate error of the registered 3D images.

In Example 2, Example 1 can further include conditioning the error of the first and second 3D images before propagating the error.

In Example 3, Example 2 can further include, wherein the conditioned error includes errors in translation in x, errors in translation in y, errors in translation in z, errors in roll, errors in yaw, errors in pitch, and errors in scale between the first and second images.

In Example 4, at least one of Examples 1-3 can further include, wherein reducing the sum aggregate of discrepancies includes using a least squares estimator between the tie points and the associated 3D points in the first and second 3D images.

In Example 5, at least one of Examples 1-4 can further include, wherein the tie points include respective tie point errors and reducing the sum aggregate of discrepancies is further determined based on the tie point errors.

In Example 6, at least one of Examples 1-5 can further include, wherein the error of the registered 3D images is an improvement over the initial errors of the first 3D image and the second 3D image.

In Example 7, at least one of Examples 1-6 can further include, wherein reducing the discrepancies include folding ground points associated with the tie points into a representation of the first and second 3D images.

In Example 8, Example 7 can further include unfolding corrections to the ground points from the representation.

In Example 9, at least one of Examples 1-8 can further include, wherein the first 3D image and the second 3D image are views of first and second geographical regions, respectively, that overlap.

Example 10 can include a non-transitory machine-readable medium including instructions that, when executed by a machine, cause a machine to perform operations for registering a first three-dimensional (3D) image to a second 3D image with error propagation, the operations comprising reducing a sum aggregate of discrepancies between respective tie points and associated 3D points in the first and the second 3D images, adjusting 3D error models of the first and second 3D images based on the reduced discrepancies to generate registered 3D images, and propagating an error of the first or second 3D images to the registered 3D image to generate error of the registered 3D images.

In Example 11, Example 10 can further include, wherein the operations further comprise conditioning the error of the first and second 3D images before propagating the error.

In Example 12, Example 11 can further include, wherein the conditioned error includes errors in translation in x, errors in translation in y, errors in translation in z, errors in roll, errors in yaw, errors in pitch, and errors in scale between the first and second images.

In Example 13, at least one of Examples 10-12 can further include, wherein reducing the sum aggregate of discrepancies includes using a least squares estimator between the tie points and the associated 3D points in the first and second 3D images.

In Example 14, at least one of Examples 10-13 can further include, wherein the tie points include respective tie point errors and reducing the sum aggregate of discrepancies is further determined based on the tie point errors.

In Example 15, at least one of Examples 10-14 can further include, wherein the error of the registered 3D images is an improvement over the initial errors of the first 3D image and the second 3D image.

Example 16 can include a memory including first and second three-dimensional (3D) images of first and second geographical regions stored thereon, and processing circuitry coupled to the memory, the processing circuitry configured to reduce a sum aggregate of discrepancies between respective tie points and associated 3D points in the first and the second 3D images, adjust 3D error models of the first and second 3D images based on the reduced discrepancies to generate registered 3D images, and propagate an error of the first or second 3D images to the registered 3D image to generate error of the registered 3D images.

In Example 17, Example 16 can further include, wherein the processing circuitry is further configured to reduce the discrepancies by folding ground points associated with the tie points into a representation of the first and second 3D images.

In Example 18, Example 17 can further include, wherein the processing circuitry is further configured to unfold corrections to the ground points from the representation.

In Example 19, at least one of Examples 16-18 can further include, wherein the first 3D image and the second 3D image are views of first and second geographical regions, respectively, that overlap.

In Example 20, at least one of Examples 16-19 can further include, wherein the processing circuitry is further configured to condition the error of the first and second 3D images before propagating the error and wherein the conditioned error includes errors in translation in x, errors in translation in y, errors in translation in z, errors in roll, errors in yaw, errors in pitch, and errors in scale between the first and second images.

Although an embodiment has been described with reference to specific example embodiments, it will be evident that various modifications and changes may be made to these embodiments without departing from the broader spirit and scope of the invention. Accordingly, the specification and drawings are to be regarded in an illustrative rather than a restrictive sense. The accompanying drawings that form a part hereof, show by way of illustration, and not of limitation, specific embodiments in which the subject matter may be practiced. The embodiments illustrated are described in sufficient detail to enable those skilled in the art to practice the teachings disclosed herein. Other embodiments may be utilized and derived therefrom, such that structural and logical substitutions and changes may be made without departing from the scope of this disclosure. This Detailed Description, therefore, is not to be taken in a limiting sense, and the scope of various embodiments is defined only by the appended claims, along with the full range of equivalents to which such claims are entitled.

Claims

1. A method for registering a first three-dimensional (3D) image to a second 3D image with error propagation, the method comprising: iteratively reducing, based on respective iterative adjustments to tie points, a sum aggregate of discrepancies between respective tie points and associated 3D points in the first and the second 3D images, each tie point of the tie points indicating a 3D location on a ground surface, a corresponding 3D point in the first 3D image and a corresponding 3D point in the second 3D image;adjusting 3D error models of the first and second 3D images based on the reduced sum aggregate to generate registered 3D images; andpropagating the adjusted error models of the first or second 3D images to the registered 3D image to generate error of the registered 3D images.

2. The method of claim 1, further comprising conditioning the error of the first and second 3D images before propagating the error.

3. The method of claim 2, wherein the conditioned error includes errors in translation in x, errors in translation in y, errors in translation in z, errors in roll, errors in yaw, errors in pitch, and errors in scale between the first and second images.

4. The method of claim 1, wherein reducing the sum aggregate of discrepancies includes using a least squares estimator between the tie points and the associated 3D points in the first and second 3D images.

5. The method of claim 1, wherein the tie points include respective tie point errors and reducing the sum aggregate of discrepancies is further determined based on the tie point errors.

6. The method of claim 1, wherein the error of the registered 3D images is an improvement over initial errors of the first 3D image and the second 3D image.

7. The method of claim 1, wherein reducing the sum aggregate includes folding ground points associated with the tie points into a representation of the first and second 3D images.

8. The method of claim 7, further comprising unfolding corrections to the ground points from the representation.

9. The method of claim 1, wherein the first 3D image and the second 3D image are views of first and second geographical regions, respectively, that overlap.

10. A non-transitory machine-readable medium including instructions that, when executed by a machine, cause a machine to perform operations for registering a first three-dimensional (3D) image to a second 3D image with error propagation, the operations comprising: iteratively reducing, based on respective iterative adjustments to tie points, a sum aggregate of discrepancies between respective tie points and associated 3D points in the first and the second 3D images, each tie point of the tie points indicating a 3D location on a ground surface, a corresponding 3D point in the first 3D image and a corresponding 3D point in the second 3D image;adjust 3D error models of the first and second 3D images based on the reduced sum aggregate to generate registered 3D images; andpropagating the adjusted error models of the first or second 3D images to the registered 3D image to generate error of the registered 3D images.

11. The non-transitory machine-readable medium of claim 10, wherein the operations further comprise conditioning the error of the first and second 3D images before propagating the error.

12. The non-transitory machine-readable medium of claim 11, wherein the conditioned error includes errors in translation in x, errors in translation in y, errors in translation in z, errors in roll, errors in yaw, errors in pitch, and errors in scale between the first and second images.

13. The non-transitory machine-readable medium of claim 10, wherein reducing the sum aggregate of discrepancies includes using a least squares estimator between the tie points and the associated 3D points in the first and second 3D images.

14. The non-transitory machine-readable medium of claim 10, wherein the tie points include respective tie point errors and reducing the sum aggregate of discrepancies is further determined based on the tie point errors.

15. The non-transitory machine-readable medium of claim 10, wherein the error of the registered 3D images is an improvement over initial errors of the first 3D image and the second 3D image.

16. A system comprising: a memory including first and second three-dimensional (3D) images of first and second geographical regions stored thereon;processing circuitry coupled to the memory, the processing circuitry configured to:iteratively reduce, based on respective iterative adjustments to tie points, a sum aggregate of discrepancies between respective tie points and associated 3D points in the first and the second 3D images, each tie point of the tie points indicating a 3D location on a ground surface, a corresponding 3D point in the first 3D image and a corresponding 3D point in the second 3D image;adjust 3D error models of the first and second 3D images based on the reduced sum aggregate to generate registered 3D images; andpropagate the adjusted error models of the first or second 3D images to the registered 3D image to generate error of the registered 3D images.

17. The system of claim 16, wherein the processing circuitry is further configured to reduce the sum aggregate by folding ground points associated with the tie points into a representation of the first and second 3D images.

18. The system of claim 17, wherein the processing circuitry is further configured to unfold corrections to the ground points from the representation.

19. The system of claim 16, wherein the first 3D image and the second 3D image are views of first and second geographical regions, respectively, that overlap.

20. The system of claim 16, wherein the processing circuitry is further configured to condition the error of the first and second 3D images before propagating the error and wherein the conditioned error includes errors in translation in x, errors in translation in y, errors in translation in z, errors in roll, errors in yaw, errors in pitch, and errors in scale between the first and second images.