Calibration for multi-camera and multisensory systems

Abstract

A method and apparatus for calibrating an image capture device are provided. The method includes capturing one or more of a single or Multiview image set by the image capture device, detecting one or more calibration features in each set by a processor, initializing each of the one or more calibration parameters a corresponding default value, extracting one or more relevant calibration parameters, computing an individual cost term for each of the identified relevant calibration parameters, and scaling each of the relevant cost terms. The method continues with combining all the cost terms once each of the calculated relevant cost terms have been scaled, determining if the combination of the cost terms has been minimized, adjusting the calibration parameters if it is determined that that the combination of the cost terms has not been minimized, and returning to the step of extracting one or more of the relevant calibration parameters.

Claims

1. A method for calibrating an image capture device, comprising: capturing one or more of a single or Multiview image set by the image capture device; detecting one or more calibration features in each set by a processor; initializing each of a one or more calibration parameters a corresponding default value; extracting one or more relevant calibration parameters; computing an individual cost term for each of the identified relevant calibration parameters; scaling each of the relevant cost terms; combining all of the cost terms once each of the calculated relevant cost terms have been scaled; determining if the combination of the cost terms has been minimized; adjusting the calibration parameters if it is determined that that the combination of the cost terms has not been minimized; and returning to the step of extracting one or more of the relevant calibration parameters.

2. The method of claim 1, wherein one of the individual cost functions comprises a Euclidian distance between two or more 3D features.

3. The method of claim 2, wherein the step of calculating the Euclidian distance further comprises the steps of: converting one or more 2D features to 3D; computing a Euclidian distance between each of two or more pairs of 3D features; computing a variance across all observations of the same dimension; and determining whether the computed variance has been minimized.

4. The method of claim 1, wherein one of the individual cost functions comprises a reprojection error.

5. The method of claim 4, wherein the step of calculating the reprojection error further comprises the steps of: converting one or more 2D features to 3D; reprojecting one or more 3D features to 2D; computing the reprojection error for one or more of the reprojected features; and determining whether the reprojection error has been minimized.

6. The method of claim 1, wherein one of the individual cost functions comprises a a difference observed in a disparity of reprojection points.

7. The method of claim 6, wherein the step of calculating the reprojection error further comprises the steps of: converting one or more 2D features to 3D; reprojecting one or more 3D features to 2D; computing the disparity of one or more reprojected points; determining between the observed and reprojected disparities; and determining whether the difference has been minimized.

8. The method of claim 1, wherein one of the individual cost functions comprises a distance between two or more epipolar feature lines corresponding to two or more matched 2D features.

9. The method of claim 8, wherein the step of calculating the distance between the two or more epipolar feature lines further comprises the steps of: extracting a fundamental matrix from the one or more calibration parameters; computing the epipolar lines corresponding to matched 2D features; computing the distance between each computer epipolar line; and determining whether the distance has been minimized.

10. The method of claim 1, wherein one of the individual cost functions comprises a sampson error.

11. The method of claim 10, wherein the step of calculating the sampson error further comprises the steps of: extracting a fundamental matrix from the one or more calibration parameters; computing the sampson error for one or more observed features from the fundamental matrix; and determining whether the sampson error has been minimized.

12. The method of claim 1, wherein one of the individual cost functions comprises a difference error between two or more parameters from a calculated essential matrix.

13. The method of claim 12, wherein the step of calculating the difference error further comprises the steps of: computing a fundamental matrix from the one or more matched 2D features; computing an essential matrix from the fundamental matrix; decomposing the essential matrix into rotation and translation elements; determining a difference between two or more of the parameters of the essential matrix; and determining whether the differences have been minimized.

14. The method of claim 1, wherein one of the individual cost functions comprises an aggregation of a dimensioning error and a reprojection error.

15. The method of claim 14, wherein the step of aggregating the dimensioning error and the reprojection error further comprises the steps of: converting one or more 2D features to 3D; computing a Euclidian distance between each of two or more pairs of 3D features; computing a difference between each of the 3D measurements and the ground truth; scaling a dimensioning error to one or more pixels using a resolution scaler based on pixel resolution; computing a reprojection error; aggregating the dimensioning error and the reprojection error determining whether the errors have been minimized.

16. A system for calibrating an image capture device, comprising: a calibration target having a known size, and further including one or more known detectable features included thereon; a calibration jig to which the calibration target, the calibration jig allowing movement of the calibration target one or more known distances relative to the image capture device; and a processor for analyzing one or more of a single or Multiview image set captured by the image capture device; detecting one or more calibration features in each set; initializing each of a one or more calibration parameters a corresponding default value; extracting one or more relevant calibration parameters; computing an individual cost term for each of the identified relevant calibration parameters; scaling each of the relevant cost terms; combining all of the cost terms once each of the calculated relevant cost terms have been scaled, determining if the combination of the cost terms has been minimized; adjusting the calibration parameters if it is determined that that the combination of the cost terms has not been minimized; and returning to the step of extracting one or more of the relevant calibration parameters.

17. The system of claim 16, further comprising: one or more extendible elements to adjust a size of the calibration target.

18. The system of claim 16, wherein the calibration target may be adjusted in accordance with its pan, tilt or rotation relative to the calibration jig.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) FIG. 1 is an example of a calibration target and jig constructed in accordance with an embodiment of the invention;

(2) FIG. 2 is a perspective view depicting preferred motions of the calibration target and jig of FIG. 1;

(3) FIG. 3 is a plan view depicting alternative flower calibration elements constructed in accordance with alternative embodiments of the present invention;

(4) FIG. 4 is a perspective view depicting implementation of a calibration process with an object of unknown size in accordance with an alternative embodiment of the invention;

(5) FIG. 5 is a flowchart depicting a generic calibration process for estimating camera parameters in accordance with an embodiment of the invention;

(6) FIG. 6 is a flowchart depicting a cost function based upon dimensional consistency in accordance with an embodiment of the invention;

(7) FIG. 7 is a flowchart depicting a cost function based upon reprojection error in accordance with an embodiment of the invention;

(8) FIG. 8 is a flowchart depicting a cost function based upon disparity error in accordance with an embodiment of the invention;

(9) FIG. 9 is a flowchart depicting a cost function based upon epipolar error in accordance with an embodiment of the invention;

(10) FIG. 10 is a flowchart depicting a cost function based upon Sampson distance in accordance with an embodiment of the invention;

(11) FIG. 11 is a flowchart depicting a cost function based upon fundamental matrix in accordance with an embodiment of the invention;

(12) FIG. 12 is a flowchart depicting a cost function based upon dimensions and reprojection in accordance with an embodiment of the invention; and

(13) FIG. 13 is a plan view depicting measurement error during calibration with a calibration target positioned at increasing distances from a camera.

DETAILED DESCRIPTION

(14) Auto Calibration within the Context of an Observed Constancy

(15) We revisit the problem of auto calibration with an advantage: if a given distance, defined as a dimension, is constant, then it can be used to improve calibration. Incorporating such information into a new cost function becomes very critical to address the issue of better and more accurate calibration. In order to define such a cost function, a series of observations have to first be made:

(16) Given two points in three-dimensional space, the distance between these two points will remain constant across multiple views taken by a camera rig. If one looks at a population of data from a calibrated stereo rig, then that population will be “well behaved” in so far as adhering to the overall disparity graph of the stereo rig. This means that the relationship between disparity and depth should follow and adhere to the overall theoretical curve. Given a hypothesis of a theoretical curve, if the data associated with real world observations does not match that curve, then the curve is fundamentally incorrect. Outliers between the theoretical curve and the data represent a good measurement of how far off the observed data is from the curve. A cost function can be defined that encapsulates information about outliers relative to the theoretical curve.

(17) Calibration Targets

(18) The first step in any calibration algorithm is to acquire images of the calibration target from one or more image acquisition angles. Two options are presented for a calibration target: The first consists of an object with a known size, preferably containing one or more patterns with detectable features; and the second consists of an object with unknown, albeit fixed size. Both options may produce reliable camera calibration.

(19) Calibration with a Known Pattern Size

(20) An example of a calibration target 110 with a known size is shown in FIG. 1. This calibration target 110 consists of one or more patterns of detectable features 120, such as a checkerboard or pattern of circular dots. Including multiple patterns on the same target is a method of obtaining additional high accuracy calibration features from a single image. Preferably, this target is mounted to a calibration jig 130 that is able to move the calibration target to multiple fixed distances and/or orientations (pan, tilt, rotation). Referring next to FIG. 2, a calibration element 240 is shown including a support element 242 removably fixable to jig 130, a central calibration element 244, and one or more extendable calibration elements 246 for adjusting a size of calibration element 240. As is shown in FIG. 2, extendible. As is further shown in FIG. 2, calibration target 240 is preferably moveable along jig 130 in a direction 248 perpendicular to a plane of a camera to be calibrated. Additionally, calibration target 240 may be adjusted (pan, tilt, rotation) according to the movement thereof in any of the corresponding directions 250, 260, 270.

(21) Alternative flower type calibration targets are shown in FIG. 3, which include not only different shapes, but also different sizes, thus further representing calibration target 240 of FIG. 2. As is further shown in FIG. 3, initial calibration target 240 may be provided with an additional outer layer, as is shown in calibration target 310 (including offset calibration patterns 315) or calibration target 320 (including aligned calibration patterns 325). In an alternative embodiment of the invention, a non-regular pattern may also be provided in order to fill a field of view of a camera to be calibrated, as is shown with calibration element 330, and extending portion 335. Variations include multiple rings of petals, petals at different orientations, and petals that are only partially expanded.

(22) Calibration with an Unknown Pattern Size

(23) Additionally, an object of unknown size may be used as the calibration target. This target again must contain detectable features, but these features need not be a part of a repeating pattern. The only constraint is that the features remain in a fixed location on calibration target in all image captures. For example, as shown in FIG. 4, a cuboid object, such a box of unknown dimensions 410 may be used as the calibration target. In this case, the corners 420 and/or edges 430 of the box are of a fixed and detectable feature. Consistent locations of these features thus allow for calibration in accordance with any of the methods or systems described in this application. Additionally, an object of unknown size may be used as the calibration target. This target, similarly to the “box” example note above, may be utilized as long as the target includes one or more easily-distinguishable features, and as long as those distinguishable features are maintained at consistent positions.

(24) Cost Functions

(25) A general framework for estimating calibration parameters is depicted in FIG. 5. As previously mentioned, the first step 510 is to acquire multiple single or multi-view image sets viewing a calibration target. The features of the calibration target, such as the centroids of a circular pattern or the corners of a checkerboard or cuboidal object, are detected in each image at step 515. The calibration parameters are then initialized to default values, which may be random, derived of hardware tolerances, or estimated from features of the detected features in closed form, at step 520. Next, the cost terms for each relevant cost function are computed at step 525 using the current estimate of the camera parameters. As is further shown in FIG. 5, step 525 inquires whether all cost terms have been analyzed. If the inquiry at step 525 is answered in the negative, and it is therefore determined that all cost terms have not yet been analyzed, processing continues at step 530 where one or more relevant calibration parameters are extracted relevant to the cost terms being currently considered. Processing continues at 535 where an individual cost term is computed, and finally to step 540 where the computer cost term is scaled or normalized. At this step, each term is either normalized, weighted, or scaled to a common unit system allowing the individual terms to be averaged, aggregated, or otherwise combined to obtain a final cost term.

(26) If, on the other hand, the inquiry at step 525 is answered in the affirmative, and it is therefore determined that all cost terms have been considered, processing instead passes to step 545 all cost terms are preferably combined. Of course, any desirable or relevant subset of all cost terms may be employed at this point. An inquiry is then initiated at step 550 as to whether the cost has been minimized. If this inquiry at step 550 is answered in the affirmative, and it is determined that cost has been minimized, then processing ends at step 560.

(27) If, on the other hand, the inquiry at step 550 is answered in the negative, and it is determined that the cost has not been minimized, processing passes to step 555 where a second inquiry is initiated to determine whether one or more particular exit criterium has been reached. Once again, if this inquiry is answered in the affirmative, and it is determined that an exit criterium has been met, processing again ends at step 560.

(28) If, on the other hand, the inquiry at step 555 is answered in the negative, and it is therefore determined that an exit criterium has not been reached, processing then passes to step 565 where one or more calibration parameters is adjusted. The camera parameters are adjusted using an optimization algorithm until the cost has been minimized or an exit criterion like maximum number of iterations has been reached.

(29) After adjustment, processing then returns to step 525 for additional computation of individual cost terms.

(30) Mathematically, the general form for online and offline calibration recovers the calibration parameters which minimize the cost, ε, over a set of N datasets:

(31) $\begin{array}{cc}{\mathrm{.Math.}}^2=\frac{1}{N}\underset{n=1}{\overset{N}{.Math.}}{\mathrm{.Math.}}_n^2& \mathrm{Equation}1\end{array}$
whereby, each dataset contains observations of matched features in the scene, and the cost of each dataset is determined by a weighted (α) combination of K terms which target individual calibration parameters:

(32) $\begin{array}{cc}{\mathrm{.Math.}}_n^2=\frac{1}{K}\underset{k=1}{\overset{K}{.Math.}}{\alpha }_k{\mathrm{.Math.}}_k^2& \mathrm{Equation}2\end{array}$

(33) Dimensional Consistency—This type of cost function utilizes the fact that the dimensions of a physical or virtual object do not vary with viewpoint or distance from the camera to estimate the camera parameters. As shown in FIG. 6, the cost function computes the Euclidean distance between a set or subset of 3D calibration features, then minimizes the variation in all observations of the same measurement in the image capture set using the current camera parameters. The parameters are then are adjusted in an iterative manner using an optimization algorithm. This approach is suitable for both single and multi-view camera setups.

(34) Referring next to FIG. 6, a process for determining a cost function based upon dimensional consistency (or constancy) is shown. The first step 610 in FIG. 6 is to acquire multiple single or multi-view image sets viewing a calibration target. The features of the calibration target, such as the centroids of a circular pattern or the corners of a checkerboard or cuboidal object, are detected in each image at step 615. The calibration parameters are then initialized to default values, which may be random, or they may be derived from camera hardware tolerances, or they may be estimated from features of the detected features in closed form, at step 620.

(35) Processing then passes to step 625 where one or more 2D features is converted into 3D features. A Euclidian distance is then computed between the 3D features in step 630, followed by step 635, where a variance across all observations of the same dimension is computed. It is then queried at step 640 whether the variance computed at step 635 has been minimized. If this inquiry is answered in the affirmative, and it is therefore determined that the variance has been minimized, processing ends at step 655. If on the other hand it is determined at step 640 that the variance has not been minimized, then processing passes to step 645 where it is further queried whether an exit criterium has been reached. If the query at step 645 is answered in the affirmative, and it is therefore determined that an exit criterium has been reached, processing ends at step 655. If, on the other hand, it is determined at step 645 that an exit criterium has not yet been reached, processing passes to step 650 where the calibration parameters are adjusted. Calibration then returns to step 625 for further processing.

(36) Mathematically, this approach can be described in the following lines:

(37) Let m.sub.L,i,j=[x.sub.L,i,j y.sub.L,i,j 1].sup.T and m.sub.R,i,j=[x.sub.R,i,j y.sub.R,i,j 1].sup.T represent the image coordinates of matched feature j in left (L) and right (R) frames at time i. The 3D coordinates of the feature M.sub.i,j=[X.sub.i,j Y.sub.i,j Z.sub.i,j 1].sup.T are computed from the singular value decomposition of:

(38) $\begin{array}{cc}{m}_{i,j,L}=\left[K_L.Math.0]M_{i,j}=H_L{M}_{i,j}{m}_{i,j,R}=K_R[R.Math.t\right]M_{i,j}=H_R{M}_{i,j}& \mathrm{Equation}3\end{array}$
where K.sub.R and K.sub.L describe the intrinsic camera parameters (focal length and principle point) and R and t describe the extrinsic parameters (rotation and translation). Given J tracked features, there exists Q dimensions between every feature pair:

(39) $\begin{array}{cc}Q=\frac{J!}{2\left(J-2\right)!}=\frac{J\left(J-1\right)\left(J-2\right)!}{2\left(J-2\right)!}=\frac{J\left(J-1\right)}{2}& \mathrm{Equation}4\end{array}$
where the dimension between any two features (j and j+l) is given as the Euclidean distance between the 3D coordinates:
L.sub.i,j,j+l=∥M.sub.i,j−M.sub.i,j+l∥  Equation 5

(40) When the same features are tracked over all I frames, the dimension (distance between features) should remain constant. Accordingly, this cost term recovers the extrinsic parameters which minimize the variance, σ.sup.2, across all dimensions:

(41) $\begin{array}{cc}\mathrm{.Math.}=\frac{2}{J\left(J-1\right)}\underset{j=1}{\overset{J}{.Math.}}\underset{l=1}{\overset{J-l}{.Math.}}{\sigma }_{j,j+l}^2& \mathrm{Equation}6\end{array}$
where the dimension variance between features j and j+l is given as a function of the mean μ.sub.j,j+l:

(42) $\begin{array}{cc}{\sigma }_{j,j+l}^2=\frac{1}{l}\underset{i=0}{\overset{I}{.Math.}}{\left(L_{i,j,j+l}-{\mu }_{j,j+l}\right)}^2& \mathrm{Equation}7\\ {\mu }_{j,j+l}=\frac{1}{I}\underset{i=0}{\overset{I}{.Math.}}{L}_{i,j,j+l}& \mathrm{Equation}8\end{array}$
This term specifically recovers the yaw element of the rotation matrix (vergence angle), but cannot recover pitch, roll, or translation very accurately. In such a case, it is pertinent to insure tight tolerances for pitch, roll or translation during a manufacturing process.

(43) Reprojection Error—This approach is one of the most widely used techniques in industry and academia for camera calibration and can be used for both single and multi-view image capture setups. The cost function in this case is based on the difference between the 2D locations of the detected calibration features and reprojection of the same features from 3D to 2D using the current camera parameter. 2D to 3D conversion (and vice versa) can leverage stereo triangulation and/or projection of known feature coordinates in 3D. The main drawback of reprojection error is that it is fundamentally front-heavy in the analysis. The data that is associated with a given target upfront have more pixel resolution and more associated details than the data in the back, since pixel resolution is higher in the front. This implies that the reprojection error will be schewed to the front of the data. However, most of the error comes from the same, fixed size patterns that are moved to the back of the field of view. However, the reprojection error in the back is significantly lower than the front, contributing to a lower per-pixel value. One work-around would be to normalize via individual camera pixel resolution, so that the reprojection error is no longer in pixels, but rather in real-world coordinate systems, such as inches or centimeters.

(44) Referring next to FIG. 7, a process for determining a cost function based upon reprojection error is shown. The first step 710 in FIG. 7 is to acquire multiple single or multi-view image sets viewing a calibration target. The features of the calibration target, such as the centroids of a circular pattern or the corners of a checkerboard or cuboidal object, are detected in each image at step 715. The calibration parameters are then initialized to default values, which may be random, derived of hardware tolerances, or estimated from features of the detected features in closed form, at step 720.

(45) Processing then passes to step 725 where one or more 2D features are converted into 3D features. 3D features are then reprojected to 2D features in step 730, followed by step 735, where a reprojection error is computed. It is then queried at step 740 whether the error computed at step 735 has been minimized. If this inquiry is answered in the affirmative, and it is therefore determined that the variance has been minimized, processing ends at step 755. If on the other hand it is determined at step 740 that the variance has not been minimized, then processing passes to step 745 where it is further queried whether an exit criterium has been reached. If the query at step 745 is answered in the affirmative, and it is therefore determined that an exit criterium has been reached, processing ends at step 755. If on the other hand it is determined at step 745 that an exit criterium has not yet been reached, processing passes to step 750 where the calibration parameters are adjusted. Processing then returns to step 725 for further processing.

(46) Mathematically, this approach can be described as:

(47) Given the matched 2D feature, m.sub.L,i,j and m.sub.R,i,j, the triangulated 3D coordinates, M.sub.i,j, are reprojected onto the image plane using (p.sub.L,i,j and p.sub.R,i,j) which expands to:

(48) $\begin{array}{cc}{p}_{L,i,j}=\left[\begin{array}{c}{x}_{L,i,j}\\ y_{L,i,j}\\ 1\end{array}\right]=\left[\begin{array}{ccc}{f}_{L,x}& 0& c_{L,x}\\ 0& f_{L,y}& c_{L,x}\\ 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{c}{X}_{i,j}\\ Y_{i,j}\\ Z_{i,j}\\ 1\end{array}\right]=\hspace{1em}\hspace{1em}\left[\begin{array}{cccc}{f}_{L,x}& 0& c_{L,x}& 0\\ 0& f_{L,y}& c_{L,x}& 0\\ 0& 0& 1& 0\end{array}\right]\left[\begin{array}{c}{X}_{i,j}\\ Y_{i,j}\\ Z_{i,j}\\ 1\end{array}\right]p_{R,i,j}=\hspace{1em}\left[\begin{array}{c}{x}_{R,i,j}\\ y_{R,i,j}\\ 1\end{array}\right]=\left[\begin{array}{ccc}{f}_{R,x}& 0& c_{R,x}\\ 0& f_{R,y}& c_{R,y}\\ 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}{r}_{11}& r_{12}& r_{13}& t_x\\ r_{21}& r_{22}& r_{23}& t_y\\ r_{31}& r_{32}& r_{33}& t_z\end{array}\right]\left[\begin{array}{c}{X}_{i,j}\\ Y_{i,j}\\ Z_{i,j}\\ 1\end{array}\right]& \mathrm{Equation}9\end{array}$
Where the f describes the focal length and c the optical center. The reprojection error cost term recovers pitch and roll by minimizing:

(49) $\begin{array}{cc}\mathrm{.Math.}=\underset{j=1}{\overset{J}{.Math.}}\underset{i=1}{\overset{I}{.Math.}}{.Math.p_{L,i,j}-m_{L,i,j}.Math.}^2+{.Math.p_{R,i,j}-m_{R,i,j}.Math.}^2& \mathrm{Equation}10\end{array}$

(50) Disparity Error—The disparity error approach depicted in FIG. 8 is closely related to the reprojection error based approach. Only, the differences in the reprojected and detected disparities are used instead of the differences in the reprojected and detected 2D features as the error term for the cost function. Disparity is defined as the difference in the 2D locations of a matched feature across multiple-views. As a result, this approach can be only used in a multi-view image capture setup.

(51) Referring next to FIG. 8, a process for determining a cost function based upon disparity error is shown. The first step 810 in FIG. 8 is to acquire multiple multi-view image sets viewing a calibration target. The features of the calibration target, such as the centroids of a circular pattern or the corners of a checkerboard or cuboidal object, are detected in each image at step 815. The calibration parameters are then initialized to default values, which may be random, derived of hardware tolerances, or estimated from features of the detected features in closed form, at step 820.

(52) Processing then passes to step 825 where one or more 2D features are converted into 3D features. 3D features are then reprojected to 2D features in step 830, followed by step 835, where a disparity of the reprojected points is computed. At step 838 a difference between the observed and reprojected disparities is determined. It is then queried at step 840 whether the difference determined at step 838 has been minimized. If this inquiry is answered in the affirmative, and it is therefore determined that the difference has been minimized, processing ends at step 855. If on the other hand it is determined at step 840 that the difference has not been minimized, then processing passes to step 845 where it is further queried whether an exit criterium has been reached. If the query at step 845 is answered in the affirmative, and it is therefore determined that an exit criterium has been reached, processing ends at step 855. If on the other hand it is determined at step 845 that an exit criterium has not yet been reached, processing passes to step 850 where the calibration parameters are adjusted. Processing then returns to step 825 for further processing.

(53) Mathematically, this approach can be described as follows:

(54) The disparity error is an extension of the reprojection error, which compares the observed disparity of a matched feature (m.sub.L,i,j and m.sub.R,i,j):
d.sub.obs,i,j=m.sub.L,i,j−m.sub.R,i,j  Equation 11
to an idealized disparity computed from the reprojection of the 3D point, M.sub.i,j, onto the image plane (p.sub.L,i,j and p.sub.R,i,j):
d.sub.ideal,i,j=[K.sub.L|0]M.sub.i,j−K.sub.R[R|t]M.sub.i,j=p.sub.L,i,j−p.sub.R,i,j  Equation 12
where the cost term recovers the pitch and roll, minimizing to:

(55) $\begin{array}{cc}\mathrm{.Math.}=\underset{j=1}{\overset{J}{.Math.}}\underset{i=1}{\overset{I}{.Math.}}{\left(d_{\mathrm{obs},i,j}-d_{\mathrm{ideal},i,j}\right)}^2& \mathrm{Equation}13\end{array}$

(56) Epipolar Error—The epipolar error is another multi-view based approach that defines cost as the distance between the 2D location of a matched calibration feature and the corresponding epipolar line determined from the current parameter estimate. The iterative process then adjusts the parameters to minimize this difference. The details of the approach are described in FIG. 9.

(57) Referring next to FIG. 9, a process for determining a cost function based upon epipolar error is shown. The first step 910 in FIG. 9 is to acquire multiple multi-view image sets viewing a calibration target. The features of the calibration target, such as the centroids of a circular pattern or the corners of a checkerboard or cuboidal object, are detected in each image at step 915. The calibration parameters are then initialized to default values, which may be random, derived of hardware tolerances, or estimated from features of the detected features in closed form, at step 920.

(58) Processing then passes to step 925 a fundamental matrix is preferably extracted from the calibration parameters. Epipolar lines corresponding to matched 2D features are then computed in step 930, followed by step 935, where a distance between each feature and line is computed. It is then queried at step 940 whether the distance computed at step 935 has been minimized. If this inquiry is answered in the affirmative, and it is therefore determined that the distance has been minimized, processing ends at step 955. If on the other hand it is determined at step 940 that the distance has not been minimized, then processing passes to step 945 where it is further queried whether an exit criterium has been reached. If the query at step 945 is answered in the affirmative, and it is therefore determined that an exit criterium has been reached, processing ends at step 955. If on the other hand it is determined at step 945 that an exit criterium has not yet been reached, processing passes to step 950 where the calibration parameters are adjusted. Processing then returns to step 925 for further processing.

(59) Mathematically, the 2D coordinates of a matched feature in left and right images, m.sub.L,i,j=[x.sub.L,i,j y.sub.L,i,j 1].sup.T and m.sub.R,i,j=[x.sub.R,i,j y.sub.R,i,j 1].sup.T, are related by the fundamental matrix, F:
m.sub.R,i,j.sup.TFm.sub.L,i,j=0  Equation 14

(60) 0$\begin{array}{cc}\left[\begin{array}{ccc}{x}_{R,i,j}& y_{R,i,j}& 1\end{array}\right]\left[\begin{array}{ccc}{F}_{11}& F_{12}& F_{13}\\ F_{21}& F_{22}& F_{23}\\ F_{31}& F_{32}& F_{33}\end{array}\right]\left[\begin{array}{c}{x}_{L,i,j}\\ y_{L,i,j}\\ 1\end{array}\right]=0& \mathrm{Equation}15\end{array}$
which can be expanded into the epipolar line equation at each i,j [16]:
(F.sub.11x.sub.L,i,j+F.sub.12y.sub.L,i,j+F.sub.13)x.sub.R,i,j+(F.sub.21x.sub.L,i,j+F.sub.22y.sub.L,i,j+F.sub.23)y.sub.R,i,j+(F.sub.31x.sub.L,i,j+F.sub.32y.sub.L,i,j+F.sub.33)=0  Equation 16
and represented in homogeneous coordinates as l.sub.R,i,j=[a.sub.i,j b.sub.i,j c.sub.i,j].sup.T.
The distance between an observed point and the corresponding line is then given:

(61) $\begin{array}{cc}{d}_R\left(m_{R,i,j},l_{R,i,j}\right)=\frac{a_{i,j}{x}_{R,i,j}+b_{i,j}{y}_{R,i,j}+c_{i,j}}{\sqrt{a_{i,j}^2+b_{i,j}^2}}& \mathrm{Equation}17\end{array}$
Where the epipolar error cost function minimizes (for all observations of all features):

(62) $\begin{array}{cc}\mathrm{.Math.}=\underset{j=1}{\overset{J}{.Math.}}\underset{i=1}{\overset{I}{.Math.}}{d}_{L,i,j}^2+d_{R,i,j}^2& \mathrm{Equation}18\end{array}$

(63) Sampson Error—This approach uses the Sampson distance which is a first order approximation of the reprojection error (algebraic error) as the cost term. Sampson distance is based on the fundamental matrix extracted from the current calibration parameters and the 2D locations of a matched feature across multiple views. As shown in FIG. 10, the iterative procedure adjusts the calibration parameters so that the Sampson distance across all image captures and all calibration features is minimized.

(64) Referring next to FIG. 10, a process for determining a cost function based upon Sampson distance is shown. The first step 1010 in FIG. 10 is to acquire multiple multi-view image sets viewing a calibration target. The features of the calibration target, such as the centroids of a circular pattern or the corners of a checkerboard or cuboidal object, are detected in each image at step 1015. The calibration parameters are then initialized to default values, which may be random, derived of hardware tolerances, or estimated from features of the detected features in closed form, at step 1020.

(65) Processing then passes to step 1025 where a fundamental matrix is preferably extracted from the calibration parameters. Sampson error for the observed features is then calculated at step 1035. It is then queried at step 1040 whether the error computed at step 1035 has been minimized. If this inquiry is answered in the affirmative, and it is therefore determined that the variance has been minimized, processing ends at step 1055. If on the other hand it is determined at step 1040 that the variance has not been minimized, then processing passes to step 1045 where it is further queried whether an exit criterium has been reached. If the query at step 1045 is answered in the affirmative, and it is therefore determined that an exit criterium has been reached, processing ends at step 1055. If on the other hand it is determined at step 1045 that an exit criterium has not yet been reached, processing passes to step 1050 where the calibration parameters are adjusted. Processing then returns to step 1025 for further processing.

(66) Mathematically, the Sampson distance is defined as:

(67) $\begin{array}{cc}{s}_{i,j}\left(m_{L,i,j},m_{R,i,j}\right)=\hspace{1em}\hspace{1em}\frac{{\left(m_{R,i,j}^T{\mathrm{Fm}}_{L,i,j}\right)}^2}{\left({\mathrm{Fm}}_{L,i,j}\right)\left(0\right)+\left({\mathrm{Fm}}_{L,i,j}\right)\left(1\right)+\left(F^T{m}_{R,i,j}\right)\left(0\right)+\left(F^T{m}_{R,i,j}\right)\left(1\right)}& \mathrm{Equation}19\end{array}$
where m.sub.L,i,j and m.sub.R,i,j represent the observed image coordinates of a matched feature and F is the fundamental matrix. The Sampson error cost function minimizes (for all observations, i, of all features, j):

(68) $\begin{array}{cc}\mathrm{.Math.}=\underset{j=1}{\overset{J}{.Math.}}\underset{i=1}{\overset{I}{.Math.}}{s}_{i,j}^2& \mathrm{Equation}20\end{array}$

(69) Fundamental Error—In this approach described in, the matched 2D locations of the features are used to estimate the fundamental matrix which can then be used to compute the essential matrix. When decomposed, the essential matrix gives the rotation and translation parameters. These estimated parameters are then compared to the current parameter set and the iterative procedure adjusts the parameters until the difference between the decomposed rotation and translation parameters and the current set of parameters is minimized.

(70) Referring next to FIG. 11, a process for determining a cost function based upon fundamental matrix is shown. The first step 1110 in FIG. 11 is to acquire multiple multi-view image sets viewing a calibration target. The features of the calibration target, such as the centroids of a circular pattern or the corners of a checkerboard or cuboidal object, are detected in each image at step 1115. The calibration parameters are then initialized to default values, which may be random, derived of hardware tolerances, or estimated from features of the detected features in closed form, at step 1120.

(71) Processing then passes to step 1125 where a fundamental matrix is computed from matched 2d features. An essential matrix is then calculated from the fundamental matrix in step 1130, followed by step 1135, where the essential matrix is preferably decomposed into rotation and translation parameters. At step 1138 a difference between the newly-calculated parameters and the initial parameters is determined. It is then queried at step 1140 whether the difference determined at step 1138 has been minimized. If this inquiry is answered in the affirmative, and it is therefore determined that the difference has been minimized, processing ends at step 1155. If on the other hand it is determined at step 1140 that the difference has not been minimized, then processing passes to step 1145 where it is further queried whether an exit criterium has been reached. If the query at step 1145 is answered in the affirmative, and it is therefore determined that an exit criterium has been reached, processing ends at step 1155. If on the other hand it is determined at step 1145 that an exit criterium has not yet been reached, processing passes to step 1150 where the calibration parameters are adjusted. Processing then returns to step 1125 for further processing.

(72) Mathematically, this approach can be described as follows:

(73) Specifically, the fundamental matrix, F, can be recovered from eight or more point correspondences, where F is related to the essential matrix E by the intrinsic camera matrices K.sub.L and K.sub.R:
F=K.sub.R.sup.−TEK.sub.L.sup.−1  Equation 21
The essential matrix E is defined up to scale w by the skew-symmetric translation matrix, [T].sub.x and the rotation matrix R:
E=Rw[T].sub.x  Equation 22
Where the skew-symmetric matrix is defined from the stereo baseline, t=[t.sub.x t.sub.y t.sub.z].sup.T, as:

(74) $\begin{array}{cc}{\left[T\right]}_x=\left[\begin{array}{ccc}0& -t_z& t_y\\ t_z& 0& -t_x\\ -t_y& t_x& 0\end{array}\right]& \mathrm{Equation}23\end{array}$
If the rotation matrix is known, then translation can be extracted from the essential matrix:
R.sup.−1E=w[T.sub.R].sub.x  Equation 24
Allowing the cost to be defines based on the difference in the translation parameters:

(75) $\begin{array}{cc}\mathrm{.Math.}=\frac{1}{3}\left({\mathrm{\Delta t}}_x^2+\Delta t_y^2+\Delta t_Z^2\right)& \mathrm{Equation}25\end{array}$

(76) Combined Measurement and Reprojection Error-based Cost function for Calibration using Checkerboard Patterns—This approach described in, uses a combination of the reprojection error described in Section 4.3 and a measurement error term that is based on the differences in known physical distances between calibration features on a calibration target and 3D distances between the same features obtained as Euclidean distances between triangulated points using the current camera parameter estimate.

(77) Referring next to FIG. 12, a process for determining a cost function based upon dimensions and reprojection is shown. The first step 1210 in FIG. 12 is to acquire multiple multi-view image sets viewing a calibration target. The features of the calibration target, such as the centroids of a circular pattern or the corners of a checkerboard or cuboidal object, are detected in each image at step 1215. The calibration parameters are then initialized to default values, which may be random, derived of hardware tolerances, or estimated from features of the detected features in closed form, at step 1220.

(78) Processing then passes to step 1225 where one or more 2D features are converted into 3D features. A Euclidian distance is then computed between the 3D features in step 1230, followed by step 1232, where a difference between the 3D measurements and a ground truth is computed. At step 1235 a dimensioning error is scaled to pixels using resolution. Processing then preferably passes to step 1237 where a reprojection error is computed, and then to step 1238 where the dimensioning and reprojection errors are aggregated. It is then queried at step 1240 whether the error aggregated at step 1238 has been minimized. If this inquiry is answered in the affirmative, and it is therefore determined that the error has been minimized, processing ends at step 1255. If on the other hand it is determined at step 1240 that the error has not been minimized, then processing passes to step 1245 where it is further queried whether an exit criterium has been reached. If the query at step 1245 is answered in the affirmative, and it is therefore determined that an exit criterium has been reached, processing ends at step 1255. If on the other hand it is determined at step 1245 that an exit criterium has not yet been reached, processing passes to step 1250 where the calibration parameters are adjusted. Processing then returns to step 1225 for further processing.

(79) Mathematically, this approach can be described as follows:

(80) If the J features detected are detected as corners on several images taken of a checkerboard calibration pattern consisting of R row and C column corners, then the total number of features is given as:
J=RC  Equation 26
Additionally, there exists U unique edges formed between every adjacent corner pair, where
U=(R−1)C+(C−1)R=2RC−(R+C)=2J−(R+C)  Equation 27
If a unique edge u in image i is defined by its reconstructed 3D endpoints M.sub.i,u.sub.1 and M.sub.i,u.sub.2, the measurement of the unique edge is given by the Euclidean distance between them,
D.sub.i,u=∥M.sub.i,u.sub.i−M.sub.i,u.sub.2∥  Equation 28
The measurement cost term is then given by the average squared error between the measurements of all unique edges and the true edge length D which are expressed in units of mm:

(81) $\begin{array}{cc}{\mathrm{.Math.}}_{\mathrm{measure}}=\frac{1}{\mathrm{UI}}\underset{i=1}{\overset{I}{.Math.}}\underset{u=1}{\overset{U}{.Math.}}{\left(D_{i,u}-D\right)}^2& \mathrm{Equation}29\end{array}$
If, m.sub.L,i,j and m.sub.R,i,j represent the j.sup.th matched 2D corner for the i.sup.th image for the left and right views, respectively, and M.sub.i,j represents the triangulated 3D co-ordinates that are reprojected onto the image plane as p.sub.L,i,j and p.sub.R,i,j by using Equation 9, the reprojection cost term is given as:

(82) $\begin{array}{cc}{\mathrm{.Math.}}_{\mathrm{reproj}}=\underset{i=1}{\overset{I}{.Math.}}\underset{j=1}{\overset{J}{.Math.}}{.Math.p_{L,i,j}-m_{L,i,j}.Math.}^2+{.Math.p_{R,i,j}-m_{R,i,j}.Math.}^2& \mathrm{Equation}30\end{array}$

(83) The two cost terms are combined by dividing the measurement term by the pixel pitch p.sub.pixel (expressed as mm/pixels):

(84) $\begin{array}{cc}\mathrm{.Math.}=\frac{{\mathrm{.Math.}}_{\mathrm{measure}}}{p_{\mathrm{pixel}}}+{\mathrm{.Math.}}_{\mathrm{reproj}}& \mathrm{Equation}31\end{array}$

(85) FIG. 13 shows a calibration target 1310 with 20 mm squares imprinted thereon in a radial pattern. The calibration target is shown at increasing distances from the camera. Superimposed on each image of calibration target 1310 are lines 1315 that mark the unique edges described in Equation 27 for each of the patterns of squares imprinted on calibration target 1310. The lines are pseudo-colored with the measurement error described in with blue representing low error and red representing high error.

(86) Calibrating the Camera Response Function

(87) The presented calibration targets can also be used in the calibration of the camera response function. The camera response function, ƒ, encapsulates the nonlinear relationship between sensor irradiance, E, and the measured intensity, Z, of a photosensitive element (pixel) over an exposure time, t:
Z.sub.ij=ƒ(E.sub.it.sub.j)  Equation 32
where i represents a spatial pixel index and j represents an image exposure index. Assuming that the function is monotonic, the camera response can be used to convert intensity to irradiance by recovering the inverse of the response, f-1 (Grossberg & Nayar, 2003):
ƒ.sup.−1(Z.sub.ij)=E.sub.it.sub.j  Equation 33
Taking the natural logarithm of both sides, the simplified notation, g, can be used to represent the log inverse function:
g(Z.sub.ij)ln(ƒ.sup.−1(Z.sub.ij))=ln(E.sub.i)+ln(t.sub.j)  Equation 34
Letting Z.sub.min and Z.sub.max define the minimum and maximum pixel intensities of N pixels and P images, g can be solved by minimizing the following quadratic objective function:

(88) 0$\begin{array}{cc}O=\underset{i=1}{\overset{N}{.Math.}}\underset{j=1}{\overset{P}{.Math.}}{\left\{w\left(Z_{\mathrm{ij}}\right)\left[g\left(Z_{\mathrm{ij}}\right)-\ln \left(E_i\right)-\ln \left(t_j\right)\right]\right\}}^2+\lambda \underset{z=Z_{\min }+1}{\overset{Z_{\max }-1}{.Math.}}{\left(w\left(z\right)g^{″}\left(z\right)\right)}^2& \mathrm{Equation}35\end{array}$
Where the second term imposes a smoothness constraint on the second derivative:
g″(z)=g(z−1)+2g(z)−g(z+1)  Equation 36
With respect to the scalar coefficient, λ, and the weighting function, w:

(89) $\begin{array}{cc}w\left(z\right)=\{\begin{array}{c}z-Z_{\min },z\le \frac{1}{2}\left(Z_{\max }+Z_{\min }\right)\\ Z_{\max }-z,z>\frac{1}{2}\left(Z_{\max }+Z_{\min }\right)\end{array}& \mathrm{Equation}37\end{array}$
Once the inverse camera response function is known, map of the scene radiance can be obtained from P images by weighting exposures, which produce intensities closer to the middle of the response function:

(90) $\begin{array}{cc}\ln \left(E_i\right)=\frac{\underset{j=1}{\overset{P}{.Math.}}w\left(Z_{\mathrm{ij}}\right)\left(g\left(Z_{\mathrm{ij}}\right)-\ln \left(t_j\right)\right)}{\underset{i=1}{\overset{P}{.Math.}}w\left(Z_{\mathrm{ij}}\right)}& \mathrm{Equation}38\end{array}$

(91) Multithreading—The iterative nature of the camera calibration process lends several opportunities for parallelization to speed up operations. As a first step, extraction of calibration features for the entire image capture set can be done in a highly parallelized manner by spawning multiple threads on specialized GPU hardware. Additionally, the cost aggregation over multiple image sets and the triangulation and reprojection operations to go from 2D to 3D coordinates and back to 2D coordinates can also be performed parallelly across multiple threads.

(92) It will thus be seen that the objects set forth above, among those made apparent from the preceding descriptions, are efficiently attained and, because certain changes may be made in carrying out the above method and in the construction(s) set forth without departing from the spirit and scope of the invention, it is intended that all matter contained in the above description and shown in the accompanying drawings shall be interpreted as illustrative and not in a limiting sense.

(93) It is also to be understood that this description is intended to coverall of the generic and specific features of the invention herein described and all statements of the scope of the invention which, as a matter of language, might be said to fall there between.