Robust stereo calibration system and method for accurate digital image correlation measurements

This application claims the benefit of U.S. Provisional Patent Application No. 61/190,846 filed Sep. 3, 2008, entitled “ROBUST STEREO CALIBRATION METHOD FOR ACCURATE DIGITAL IMAGE CORRELATION MEASUREMENTS” and having the same inventors as present, which is incorporated herein by reference for all purposes.

Certain aspects of the presently disclosed technology were developed with funding from the National Science Foundation through grant CMS-0556362. Therefore, the government may retain certain rights in this invention.

The accurate measurements of full-field surface shape, displacements and strains for objects that are immersed in a media that is different from where the optical measurement system is provided, is of great interest in a variety of applications. For example, measuring deformation of biological materials over an extended period of time, such as vascular arteries, usually requires that the material be fully submerged in an adventitial bathing media with controlled temperature and pH. According to Snell's law, when light travels through different media, it refracts at the interface. Sutton et al. (Development of a Methodology for Non-Contacting Strain Measurements in Fluid Environments Using Computer Vision, Optics and Lasers in Engineering. (32) 2000; 367-377) calculated and measured the strain errors introduced by refraction when a camera in air measures a planar specimen submerged underwater. This work has clearly shown that strain errors increase rapidly as the angle between the optical axis and the interface changes. When the angle is 12°, the strain is as high as 1%.

Thus, there remains a need for an accurate and effective way to correct distortion induced by refraction and to calibrate measurement systems for multimedia applications.

According to its major aspects and briefly recited, the presently disclosed technology concerns a system and method for determining the shape and position of a multi-media interface and correcting the distortion introduced by refraction that happens at such interface(s). While the present technology applies especially well to optical measurement systems, it is not limited to optical measurement systems.

Certain aspects of the presently disclosed technology concern a system and method of calibrating a measurement system provided in a given media, while making several “readings” for calibration or, in the case of optical systems, images of a physical object in another media. The images are taken when the physical object undergoes in-plane translations, in-plane rotations and out-of-plane tilting with respect to the interface. From these images, the corresponding features on the object are extracted by image registration methods. Then, a mathematical algorithm is applied to these feature points to determine the shape and position of the interface. The resulting interface can then be used to reconstruct the three dimensional distortion-free locations of the object from its images when it is submerged in a media that is different from where the measurement system is and when it deforms. Therefore, the distortion caused by refraction at interfaces is corrected and the measurement system is able to quantify accurate motion and deformation for applications in multi-media.

An advantage of some embodiments of the invention is the preserved accuracy of measurement systems when they are applied to multi-media applications. This advantage allows users to set up measurement systems in their preferred media or environment and use them to measure what is happening in the other media without sacrificing any accuracy.

Another advantage of some embodiments of the invention is that the subject techniques are not based on or derived from a particular shape of interface(s) as long as the shape can be modeled using parametric representations. This permits the determination of interfaces that may be unusual in special applications, including multiply curved interfaces, singly curved interfaces and planar interfaces.

Still another advantage of some embodiments of the invention is that the present method is not based on or derived for applications with only two different media (i.e. one medium where the measurement system is and the other medium where the object is). Between these two media, there could be multiple media as long as these media have distinct interfaces that could be parameterized. This makes the present technology useful in varied applications.

Yet another feature of some embodiments of the present invention is that it allows measurement systems to be used to measure what is happening in the media where special-made measurement systems otherwise have to be employed, e.g., standard optical systems cannot work in an environment filled with erosive chemicals or with high temperatures. Once the harmful environment is restricted by a transparent material, the present method could be applied and make standard optical systems work without loss of accuracy in such a circumstance. This advantage extends regular measurement systems into broader application areas. Also, typically as measurement systems' specialty increases, the cost increases more rapidly. Here, the cost of a specialized measurement system may be obtained at a lower incremental cost.

In one exemplary embodiment of the present subject matter, a method for calibrating a measurement system is based on determination of interface parameters representative of an interface provided between the measurement system and a physical object to be measured. A first plurality of images of a calibration object are obtained, wherein the first plurality of images are obtained by a measurement system comprising two or more image capture devices that are provided in the same medium as the calibration object. A second plurality of images of the calibration object are obtained, wherein the second plurality of images obtained while the measurement system is in a different medium than the calibration object, the different media being separated from one another by the interface, and the second plurality of images obtained by the two of more image capture devices while the calibration object is provided at a respective plurality of different locations and orientations relative to said interface. Data corresponding to measured positions of selected image points from the first and second pluralities of images of the calibration object are extracted and electronically stored. Initial values of interface parameters from data corresponding to measured positions of selected image points from the first plurality of images are electronically calculated. Definitions of modeled positions of image points based on the initial values of the interface parameters are provided. An error function defined by the difference between modeled positions and measured positions of selected image points from the second plurality of images is defined and electronically stored.

In another exemplary embodiment of the present technology, a measurement system includes two or more image capture devices, a first memory, a second memory and a processor. The two or more image capture devices are configured to obtain a first plurality of images of a calibration object when provided in the same medium as the measurement system and a second plurality of images of the calibration object when provided in a different medium as the measurement system and separated by an interface. The first memory stores data corresponding to measured positions of selected image points from the first and second pluralities of images of the calibration object. The second memory stores software in the form of computer-executable instructions. The processor is coupled to the first and second memories and is configured to selectively implement the computer-executable instructions stored in the second memory to process the data stored in the first memory. The processor implements the computer-executable instructions stored in said second memory in order to implement the functions of: calculating initial values of interface parameters from data corresponding to measured positions of selected image points from the first plurality of images; providing definitions of modeled positions of image points based on said initial values of the interface parameters; and establishing an error function defined by the difference between modeled positions and measured positions of selected image points from the second plurality of images.

In a still further exemplary embodiment of the present technology, a method of processing an image taken of a physical object provided in a different physical medium than a measurement system to correct for refraction at an interface between the different media. In accordance with such method, at least one image of a physical object is obtained by a measurement system comprising two or more image capture devices that are provided in a different physical medium than the measurement system and separated therefrom by an interface. Measured image data points associated with the at least one image are stored in an electronic database. The measured image data points are electronically processed to determine actual image data points by applying a determined interface filter function to the measured image data points. Ultimately, the actual image data points are provided as visual output to a user.

Additional aspects and advantages of the present subject matter are set forth in, or will be apparent to, those of ordinary skill in the art from the detailed description herein. Also, it should be further appreciated that modifications and variations to the specifically illustrated, referred and discussed features and elements hereof may be practiced in various embodiments and uses of the present subject matter without departing from the spirit and scope of the subject matter. Variations may include, but are not limited to, substitution of equivalent means, features, or steps for those illustrated, referenced, or discussed, and the functional, operational, or positional reversal of various parts, features, steps, or the like.

Still further, it is to be understood that different embodiments, as well as different presently preferred embodiments, of the present subject matter may include various combinations or configurations of presently disclosed features, steps, or elements, or their equivalents (including combinations of features, parts, or steps or configurations thereof not expressly shown in the figures or stated in the detailed description of such figures). Additional embodiments of the present subject matter, not necessarily expressed in the summarized section, may include and incorporate various combinations of aspects of features, components, or steps referenced in the summarized objects above, and/or other features, components, or steps as otherwise discussed in this application. Those of ordinary skill in the art will better appreciate the features and aspects of such embodiments, and others, upon review of the remainder of the specification, including the Detailed Description of Preferred Embodiments accompanied by the following drawings.

FIG. 1 is schematic representation defining incident and refracted rays and angles and illustrating how light refracts at the interface between two different media;

FIG. 2 is a schematic representation of a stereo camera system demonstrating the source of pseudo strains caused by refraction;

FIG. 3 is a graph showing real and reconstructed grid spacing and corresponding angles for normal and shear strains in grid spacing;

FIG. 4 is a schematic representation of a stereo camera system with an interface between different media, as would be present during an interface calibration process in accordance with the presently disclosed technology;

FIG. 5 is a graphical representation of a transmitted light path from a grid point (Q) in a first medium (W), through two parallel planar interfaces associated with a second medium (G) to a focus point (O) of a measurement system in a third medium (A);

FIG. 6 is a flow chart of exemplary method steps in a process for calibrating a measurement system in accordance with aspects of the presently disclosed technology;

FIG. 7 is a schematic representation of aspects associated with three-dimensional spatial position reconstruction in accordance with the presently disclosed technology;

FIG. 8 a graphical representation of a transmitted light path from a grid point (Q) in a first medium (W), through two parallel planar interfaces associated with a second medium (G) to a focus point (O) of a measurement system in a third medium(A) for three-dimensional spatial position reconstruction of submerged object points;

FIGS. 9A and 9B provide respective graphs showing a calibration grid's image on two sensor planes at a 1 st exemplary view for interface calibration;

FIGS. 10A and 10B provide respective graphs illustrating a calibration grid at 9 th exemplary view in 3-D space for error estimation, with FIG. 10A using the results of interface parameters without errors on sensor planes (as in Table 4), and FIG. 10B using the results of interface parameters with errors on sensor planes (as in Table 5);

FIG. 11 shows the real position and reconstructed position of a measured object based on resulting calibration parameters from simulation, with data being shown without errors for views 18 and 34 in Table 10;

FIG. 12 shows the real position and reconstructed position of a measured object based on calibration parameters from simulation, with data shown with image plane errors of ±0.05 pixels for views 18 and 34 in Table 10;

FIGS. 13A and 13B show pseudo normal strains of reconstructed grid spacing with and without considering refraction;

FIGS. 14A and 14B show pseudo shear strains in reconstructed grid spacing with and without considering the effect of refraction for (a) rotation about the x-axis, θ x , and (b) rotation about the y-axis, θ y ; and

FIG. 15 is a schematic representation of exemplary hardware components for implementing aspects of the subject calibration technology.

It is to be understood by one of ordinary skill in the art that the present discussion is a description of preferred embodiments only, and is not intended as limiting the broader aspects of the present invention. Aspects of the present invention may be described herein in connection with a stereo optical measurement system, partly for convenience and partly because the present invention works well with optical measurement systems, but the present invention is not limited to optical measurement systems. Any measurement system that has distortions caused by refraction can benefit from the present method.

According to Snell's law, when light travels through two different media, it refracts at the interface in such a way that the incident and refracted rays meet two conditions. Defining the incident ray, ī, the refracted ray, t and the interface normal, n , at the point of intersection (see FIG. 1 ), the first condition requires that all three vectors must lie in the same plane. This can be written mathematically in the form

( i × n )· t = 0   (1)

where “×” denotes a cross-product and “·” denotes a scalar product. The second condition relates the angle of refraction to the angle of incidence. Defining the incident angle θ i , the refracted angle θ t , the indices of refraction for two media as η i and η t , the equation is written

η i sin θ i =η t sin θ t .   (2)

The refractive index η for air is 1.0 and for water is 1.33.

In order to demonstrate how spatial distortions, i.e. pseudo strains, are introduced by refraction, numerical simulations have been conducted in a manner that reflects the imaging process for a typical stereo vision system viewing a submerged object, such as represented schematically in FIG. 2 . To perform the simulation, it is assumed that the following parameters are known a priori: (a) all intrinsic and extrinsic camera parameters for both cameras; (b) indices of refraction for all media; (c) interface orientation and location; and (d) orientation and position of a calibration object (e.g., a calibration target with nine grid points).

The simulation procedure is as follows. First, the center of each of nine grid points on a planar grid submerged in water at a pre-specified distance from an air-water interface is imaged onto the sensor plane of both stereo cameras while considering refraction. Second, the matching locations on both sensor planes for each grid point are projected to obtain a three-dimensional spatial position without considering refraction. Finally, by comparing the actual grid shape with the ‘unrefracted’ reconstruction, the normal strains and the shear strains in the grid spacing are obtained. The normal strain is defined for each pair of grid points using the definition of average engineering normal strain, the initial grid spacing, and the reconstructed three-dimensional grid point locations to define the final grid spacing. Similarly, shear strain is defined for each pair of orthogonal grid lines using the initial angle (π/2) and the final angle obtained using the reconstructed three-dimensional grid point locations for all three grid points. FIG. 3 shows the real, reconstructed grid spacing and the corresponding angles.

The process noted above was performed for a range of grid rotation angles, φ, relative to the interface. Here, the angle (φ=0°) corresponds to the grid being placed parallel to the planar interface. Table 1 shows the pseudo normal strains and the pseudo shear strains of grid spacing as a function of the angle, φ. As shown in Table 1, when φ≧5°, the largest pseudo normal strain ε nn >500×10 −6 , errors that should be minimized when attempting to develop advanced measurement methods with high accuracy.

The subject stereo calibration methodology determines the orientation and position of the interface, as well as all stereo system camera parameters, during a two-stage calibration process. Such techniques are based on bundle adjustment and an image-based optimization metric implemented using Levenburg-Marquardt (LM) non-linear algorithm.

The present techniques utilize a method of modeling the interface with a plurality of parameters. In one exemplary embodiment, the interface is modeled as a plane defined by parameters a, b and d as

ax+by +(√{square root over (1 −a 2 −b 2 )}) z+d= 0,   (3)

where x, y, and z are coordinates in a coordinate system and a, b and d are interface parameters. In yet another embodiment, the present method is not limited to planar interfaces and the interface could be of any shape as long as it could be modeled using a parametric representation.

The present method allows multiple interfaces in the optical path from the measurement system to the object. For example, in order to make sure that the interface is constant during calibration and measurement process, a piece of transparent material (as medium G) having thickness, TH, is added between medium A where the measurement system is and medium W where the object is. This additional, parallel planar interface is included in the development since many practical applications will require viewing through a transparent containment vessel. In such a configuration, there would be two interfaces: the interface A/G between medium A and medium G and the interface G/W between medium G and medium W. In such configuration, the two interfaces have the same orientation. The interface A/G is:

ax+by+√ {square root over (1 −a 2 −b 2 )} z+d−TH /√{square root over (1 −a 2 −b 2 )}=0,   (4a)

and the interface G/W is:

ax+by+√ {square root over (1 −a 2 −b 2 )} z+d= 0.   (4b)

The normal n of both interfaces is:

n =( a, b, √ {square root over (1 −a 2 −b 2 )}).   (4c)

After the interface parameters are determined, the accurate 3-D positions of immersed objects from images could be reconstructed by taking account of interfaces like this: (1) Project rays according to the measurement system and the image points of the object. For example, if the measurement system is modeled as a pinhole, the projective ray is from the focal point to points in the sensor plane. The present model is not limited to pinhole-modeled measurement systems, and it can apply to any measurement systems as long as the 3-D location of the object could be reconstructed from its measurement; (2) When the rays hit the interfaces, they refract according to the resulting shape and position of the interface; (3) Pairs of corresponding rays are triangulated and the true 3-D positions of the object points are obtained. After this, the refraction distortion caused by light traveling through media A, media G and media W can be corrected simultaneously.

In one exemplary embodiment, the measurement system is a stereo optical system with two cameras that are modeled as pinhole devices. As noted before, the current invention does not restrict its application to optical systems, and the current invention is valid in any other measurement systems capable of acquiring images through refractive media. In one exemplary embodiment, the method uses a planar grid of known spacing as the calibration object. In another embodiment, calibration can be performed with calibration objects such as cubes, rectilinear parallelepipeds. In yet another embodiment, a calibration object with identifiable features can be translated or rotated by known amounts.

A schematic of a typical setup of a stereo optical system is shown in FIG. 2 . The world coordinate system (WCS) with axes (X w , Y w , Z w ) is coincident with camera 1 's local coordinate system, which is not necessary. Each camera's local coordinate system (camCS) is defined by a rotation matrix [R] and a translation vector {T} to the WCS. Any point p(x cam , y cam , z cam ) defined in a local coordinate system could be transformed to point P(x w , y w , z w ) in WCS using the following form

For camera 1 , [R] is the identity matrix and {T} is a zero vector. For camera 2 , [R] is the rotation matrix and {T} is the translation vector from camCS 2 to WCS. [R] and {T} for both cameras are called extrinsic parameters of the stereo optical system. Similarly, the calibration grid's local coordinate system, GCS, has axes X g , Y g and Z g , with Z g perpendicular to the grid. Any point p(x g , y g , 0) on the grid is transformed to the WCS using equation (6a)

[ x w y w z w ] = [ R g ] ⁡ [ x g y g 0 ] + [ D x D y D z ] ( 6 ⁢ ⁢ a )

where [D] is the translation vector between the WCS and the GCS defined in the WCS, and the rotation matrix can be defined in terms of Eular angles γ x , γ y , and γ z as

[ R g ] = [ ⁢ cos ⁢ ⁢ γ y 0 sin ⁢ ⁢ γ y 0 1 0 - sin ⁢ ⁢ γ y 0 cos ⁢ ⁢ γ y ] ·   [ ⁢ 1 0 0 0 cos ⁢ ⁢ γ x - sin ⁢ ⁢ γ x 0 sin ⁢ ⁢ γ x cos ⁢ ⁢ γ z ] · [ cos ⁢ ⁢ γ z - sin ⁢ ⁢ γ z 0 sin ⁢ ⁢ γ z cos ⁢ ⁢ γ x 0 0 0 1 ] . ( 6 ⁢ b )

For the remaining description of the invention, γ x , γ y , γ z , D x , D y and D z are called the grid's parameters.

A point p(x cam , y cam , z cam ) on the grid is projected to point P sen (Xs, Ys) on the sensor plane by

[ X s Y s 1 ] = [ M p ] ⁡ [ x cam / z cam y cam / z cam 1 ] ( 7 ⁢ a )

where [M p ] is the camera's projection matrix which is a combination of the camera's intrinsic and extrinsic parameters as

[ M p ] = [ f x skew C x 0 f y C y 0 0 1 ] ⁢ ⁢ and ⁢ ⁢ f x = fS x fy = fS y . ( 7 ⁢ b )

In equation (7b), the projection matrix [M p ] is composed of f x , f y , C x , C y , and skew. The lens focal length is f and the scale factors are S x and S y . The center in the sensor plane is (C x , C y ) and skew is the skew factor.

To account for lens distortion, the relationship between the distorted sensor location (X s d , Y s d ) and the undistorted sensor location (X s , Y s ) can be modeled by using a lens distortion factor κ

X s d = X s 1 + κ ⁢ ( X s ) 2 + ( Y s ) 2 ; Y s d = Y s 1 + κ ⁢ ( X s ) 2 + ( Y s ) 2 ( 7 ⁢ c )

and an approximate inverse of equation (7c) can be written as

X s = 2 ⁢ X s d 1 + 1 - 4 ⁢ κ ⁡ [ ( X s d ) 2 + ( Y s d ) 2 ] ; ⁢

⁢ Y s = 2 ⁢ Y s d 1 + 1 - 4 ⁢ κ ⁡ [ ( X s d ) 2 + ( Y s d ) 2 ] . ( 7 ⁢ d )

Lens distortion could also be corrected by using other prior lens distortion removal methods before calibration, for example as disclosed in Shreier et al., Advances in light microscope stereo vision, Expl Mechanics, 2004, 44, 278-288 or Sutton et al., Strain field measurements on mouse carotid arteries using microscopic three-dimensional digital image correlation, J. Biomed. Mat. Res. Part A, July 2007, 178-190, DOI:10.1002/jbm.a.31268, both of which are incorporated herein by reference for all purposes. In this way, all the calibration processes will be implemented on images that are free of lens distortion.

Equations (6a), (6b), and (7a) demonstrate that any grid location defined in GCS by (mΔx, nΔy, 0), where m and n are non-negative integers as grid indices, is transformed into the sensor plane. The parameters involved in the transformation include: (a) six unknown extrinsic parameters (γ x , γ y , γ z , D x , D y , D z ) for the grid; (b) a rotation matrix [R] and a translation vector {T} as the camera's extrinsic parameters; and (c) five of the camera's intrinsic parameters (f, S x , S y , skew, κ).

Suppose I M mo — air (X s mo — air , Y s mo — air ) or I M m mo — air (X s , Y s ) represents an image point on the sensor plane when the calibration grid is in air without refraction. The quantity I M mo — air (X s , Y s ) is obtained by (a) transforming points in the GCS to the WCS using equation (6), transforming points in the WCS to the camCS using equation (5), and (c) transforming points in the camCS into the sensor plane using equations (7a) and (7b). Therefore, I M mo — air (X s , Y s ) is modeled in terms of the unknown parameters: stereo cameras intrinsic and extrinsic parameters and grid parameters, together with other given parameters, such as grid spacing and indices.

As shown in FIG. 4 , an object point Q is imaged onto the sensor planes of both cameras 1 and 2 with coordinates of I M1 and I M2 in the WCS, respectively. FIG. 4 shows the light path followed by the principal ray, along with the refraction that occurs at each interface.

FIG. 5 shows the light path from a grid point Q immersed in fluid to one of the stereo cameras, where all components of the principal ray are in a common plane and the sensor coordinates for I M have been transformed into the WCS using the inverse of equations (7) and (5). Light travelling from Q is refracted at points B X and B S where the ray intersects an optical interface. The refraction angles (α, β, ξ) satisfy equation (2) in Snell's law, with V 3 (QB x ), V 2 (B x B s ) and V 1 (I M O) (representing vectors from Q to B x , B x to B s , and I M to O, respectively).

In order to understand properties of the light path, projection points have been introduced: (a) Q′, O′, and I M ′ are projections of Q, O, and I M onto the interface G/W; and (b) Q″, O″, I M ″ and B X ,″ are projections of Q, O, I M , and B x onto the interface A/G.

According to Snell's law, it is readily shown that (a) I M , O, B s , B x , Q, and their projections onto both interfaces are all on the same plane; (b) O″, I M ″, B s , B x ″, and Q″ are on the same straight line; and (c) O′, I M ′, B x , and Q′ are on another straight line. As developed in the following section, the relationship between Q on the calibration grid and its corresponding image point on the sensor plane is divided into three steps: (a) determine B x as a function of Q; (b) obtain B S as a function of B x ; and (c) define image position (X s , Y s ) in terms of B s .

Since the relationship between the refraction angles is non-linear, expressions in terms of Q and n that may be solved iteratively using estimates for (α, β, ξ) based on points measured on the sensor plane. The approach allows one to obtain numerical solutions without solving non-linear equations. The numerical simulations disclosed herein will show that the proposed optimization process converges well and the results are accurate and reliable, even in the presence of Gaussian random intensity noise.

Establishment of a relationship between Q on the calibration grid and I M (X s , Y s ) on the sensor plane can be divided into three steps: (a) Model of B x in terms of Q; (b) Model of B s in terms of B x ; (c) Model of I M model (X s model , Y s model ) in terms of B s . (Modeled parameters may be referred to herein as “model” or “mo” for short.) Since α, β, ξ are nonlinearly related according to Eq. (2), it is impossible to determine B x , B s or I M model (X s model , Y s model ) in terms of Q and n directly due to no explicit analytical solutions for nonlinear equations. In the present invention, α, β, ξ are parameters to be estimated by points measured on the sensor plane. A unique aspect of this invention is that this approach eliminates the need to solve nonlinear equations, as evidenced by simulations that show the optimization process converges well and the results are accurate and reliable.

Suppose I M measured (X s measured , Y s measured ) is the point measured on the sensor plane and it has been transformed to I M measured (x w , y w , z w ) in WCS by Eq. (3) and Eq. (4). (Measured parameters may be referred to herein as “measured” or “me” for short.) Then the refraction angle α, β, ξ could be expressed as,

β = a ⁢ ⁢ cos ⁡ ( ( OI M measured ) · n ⇀  OI M measured  ) ( 8 ⁢ a ) α = a ⁢ ⁢ sin ⁡ ( η A η W ⁢ sin ⁢ ⁢ β ) ( 8 ⁢ b ) ξ = a ⁢ ⁢ sin ⁡ ( η A η G ⁢ sin ⁢ ⁢ β ) ( 8 ⁢ c )

where η A , η G and η W are indexes of refraction of medium A, medium G and medium W. ∥ ∥ is L 2 matrix norm over magnitude for a vector.

To obtain the relationship between B x (x Bx , y Bx , z Bx ) and Q, it is helpful to first note that B x is on the G/W interface. As shown in FIGS. 4 and 5 , this gives rise to the following:

ax Bx +by Bx +√{square root over (1 −a 2 −b 2 )} z Bx +d= 0   (9a)

V 1 measured , V 3 model and n are on the same plane, such that:

[( I M measured O )×( QB x )]· n = 0   (9b)

The refraction angle relation is:

∥ Q″B x ∥ 2 =[ tan(α)*∥ Q′Q∥] 2   (9c)

As Q′, O′ could be substituted using Q, O, n , the analytical model of B x can be obtained in terms of O, I M me , Q, n by simultaneously solving equations (9a), (9b) and (9c). From Eq. (9c), there would be two possible solutions of B x with the same absolute value but opposite sign in B x 's x or y coordinate. Since O′, B x and Q′ are on the same straight line, and three points have to be in order with B x in the middle, three components of the vector Q′B x have the same sign as those of Q′O′, respectively which serves as a convenience condition to determine the correct sign for B x .

In determining a model for B s , it is helpful to note that B s satisfies the following three conditions:

(a) B s is on the interface A/G such that

ax Bs +by Bs +√{square root over (1 −a 2 −b 2 )} z Bs +d+TH /√{square root over (1 −a 2 −b 2 )}=0   (10a)

(b) V 1 measured , V 2 model and n are on the same plane such that

[( OI M measured )×( B x B s )]· n = 0   (10b)

(3) the refraction angle relation is:

∥ B x ″B s ∥ 2 =[ tan(ξ)* TH] 2   (10c)

As Q″, O″, B x ″ could be substituted using Q, O, n , TH, the analytical model of B s can be obtained in terms of O, I M me , Q, n , TH by solving equations (10a), (10b) and (10c) simultaneously. As O″, B s , and Q″ are on the same straight line, and three points have to be in order with B s in the middle, three components of the vector Q″B s have the same sign as those of the vector Q″O″, respectively, which can also be used as a condition for the program to pick up the correct sign for B s .

In determining a model of I M model (X s model , Y s model ), suppose I M model (x w , y w , z w ) is the modeled image point of Q in camCS. B s obtained from the previous steps is transformed to I M model (X s model , Y s model ) or I M model (X s , Y s ) by Eq. (7c). Therefore I M model (X s model , Y s model ) is modeled in terms of the unknown interface parameters, grid parameters, and stereo cameras' intrinsic and extrinsic parameters, together with measurements on the sensor plane I m me (X s me , Y s me ) and other given parameters, such as refractive indices, vessel wall thickness, grid spacing and indices.

The present method uses an optimization method to obtain the optimal parameters of interfaces. The method does not restrict itself to a particular optimization method and the optimization method could be any kind In one exemplary embodiment, the present method employs the least square optimization method and bundle adjustment. Suppose χ 2 is the error function on two sensor planes and it is defined by the difference between I M mo (X s mo , Y s mo ) and I M me (X s me , Y s me ) as follows:

χ 2 = ∑ k = 1 2 ⁢ ∑ i = 1 NView ⁢ ∑ j = 1 MPNTS ⁢ [ ( X s mo ⁢   i , j , k ⁢ - ⁢ X s ⁢ ⁢ i , j , k me ) 2 + ( Y s mo ⁢   i , j , k ⁢ - ⁢ Y s ⁢ ⁢ i , j , k me ) 2 ] ( 11 )

where the summation is among all the points on the calibration grid, all the views of the grid and two sensor planes from two cameras. The unknown parameters that need to be optimized are the interface parameters a, b, d, the grid parameters for each view, and the stereo's intrinsic and extrinsic parameters. The calibration problem is not reduced to minimization of χ 2 .

Considering the uniqueness and convergence problems of any non-linear optimization process, this work proposes a calibration methodology that separates the process into two successive general parts: (a) stereo system calibration to obtain stereo cameras' intrinsic and extrinsic parameters, and (b) interface calibration to obtain interface parameters that determine orientation and position of each planar interface. The stereo calibration system uses a first plurality of images when the calibration grid is in air and uses I M mo — air (X s , Y s ) as its image point model on sensor planes. The interface calibration uses a second plurality of images when the calibration grid is submerged in fluid and in its I M mo (X s , Y s ), the stereo cameras' intrinsic and extrinsic parameters are taken as given. The interface parameters a, b, and d are to be determined using a method of bundle optimization such as defined in equation (11). Additional details regarding more particular exemplary steps of the calibration process are now provided and are discussed with reference to FIG. 6 .

Referring now to FIG. 6 , step 600 involves obtaining a first plurality of images for interface parameter initialization while cameras in a measurement system and a calibration object (e.g., grid) are in the same medium (e.g., Air-A). Initial values of the interface parameters are then calculated in step 606 based on this first plurality of images. For example, in one embodiment, the system's intrinsic and extrinsic parameters R, T and M p are determined in the single media A by any standard calibration procedures. The interface thickness TH and the grid spacing Δx & Δy are measured. They are all known when calibrating the interface parameters. The present invention does not restrict itself to this procedure. All the unknown parameters could be optimized at a single time or step by step as preferred. Additional details of the first step (a) in the subject calibration may be implemented in accordance with the calibration process disclosed in U.S. Pat. No. 7,133,570 (Schreier et al.), which is hereby incorporated by reference herein for all purposes.

In accordance with a second step 602 , two image capture devices (e.g., digital cameras) take high resolution images simultaneously of a common calibration object (e.g., a calibration grid) at varied locations and orientations relative to the interface. The grid may experience as much variation in orientation as possible including in-plane translations, in-plane rotations and out-of-plane tilting. The step of obtaining high resolution images is repeated several times and a series of image pairs corresponding to a second plurality of images are acquired.

Exemplary step 604 involves extracting measured image point positions I M measured (X s measured , Y s measured ) on the sensor plane. Additional processing of such extracted measured image point positions, as indicated in step 606 involves transforming the measured image point positions (in pixels) to the world coordinate system (in mm) by Eq. (5) and Eq. (7), thus calculating initial values of interface parameters based on the the first plurality of images (e.g., the images obtained when the measurement system and the calibration object are in the same medium). Step 606 may also involve removing lens distortion. In step 608 , the image point position I M model (X s model , Y s model ) for all of the known grid positions is predicted using a combination of Eq. (8), (9) and (10) and transformed to sensor planes as I M mo (X s , Y s ) by equation (7a). In step 610 , the error function χ 2 is established by Eq. (11). The error function may be minimized to obtain the optimal values of a, b and d (the interface parameters) together with grid parameters D x i , D y i , D z i , γ x i , γ y i , γ z i for each view i. The results of the optimization are used to obtain optimized interface parameters, as indicated in step 612 .

A modified LM non-linear optimization algorithm may be employed to determine the interface parameters a, b and d together with the grid parameters for all views. LM is an iterative procedure that is used widely to solve non-linear least-square problems and has become a standard technique used extensively in many computer vision problems. Letting ω be a vector composed of all the unknowns to optimized, ω=[D x 1 , D y 1 , D z 1 , γ x 1 , γ y 1 , γ z 1 , . . . , D x i , D y i , D z i , γ x i , γ y i , γ z i , . . . , a, b, d] T . H is the Hessian matrix and

H ⁡ ( i , j ) = ∂ 2 ⁢ χ 2 ∂ ϖ i ⁢ ∂ ϖ j ⁢ | ϖ .

J is the gradient of χ 2 and

To ensure that the non-linear optimization process converges to the global minimum for χ 2 , reasonable initial values for ω may be determined. For example, the interface parameters a, b and d can be roughly estimated by taking images of the calibration grid or a speckled planar specimen located near the interface and approximately parallel to the interface. Based on the calibrated values for the stereo camera's intrinsic and extrinsic parameters, the calibration grid or the planar specimen's spectral position can be reconstructed and used to estimate the interface parameters. The grid's parameters are estimated based on the rough estimation of the interface parameters by taking refraction into consideration.

It should be appreciated that the calibration method discussed above can be implemented by a measurement system including a special-purpose computer control system for implementing the subject calibration features and steps. An example of such a measurement system 500 is shown in schematic form in FIG. 15 . More particularly, the exemplary measurement system 500 includes two or more image capture devices (shown in FIG. 15 as first and second cameras 502 and 504 ). Computer control system 512 may generally include such components as at least one memory/media element or database for storing data and software instructions as well as at least one processor. In the particular example of FIG. 15 , a processor(s) 522 and associated memory/media elements 524 a, 524 b and 524 c are configured to perform a variety of computer-implemented functions (i.e., software-based data services). At least one memory/media element (e.g., element 524 b in FIG. 15 ) is dedicated to electronically storing software and/or firmware in the form of computer-readable and executable instructions that will be implemented by the one or more processor(s) 522 . Other memory/media elements (e.g., memory/media elements 524 a, 524 c ) are used to store data which will also be accessible by the processor(s) 522 and which will be acted on per the software instructions stored in memory/media element 524 b . The various memory/media elements of FIG. 15 may be provided as a single or multiple portions of one or more varieties of computer-readable media, such as but not limited to any combination of volatile memory (e.g., random access memory (RAM, such as DRAM, SRAM, etc.) and nonvolatile memory (e.g., ROM, flash, hard drives, magnetic tapes, CD-ROM, DVD-ROM, etc.) or any other memory devices including diskettes, drives, other magnetic-based storage media, optical storage media and others.

Although FIG. 15 shows three separate memory/media elements 524 a , 524 b and 524 c , the content dedicated to such devices may actually be stored in one memory/media element or in multiple elements, any such possible variations and other variations of data storage will be appreciated by one of ordinary skill in the art.

In one particular embodiment of the present subject matter, a first portion of memory/media 524 a is configured to store input data for the subject calibration system and related methods. Input data stored in memory/media element 524 a may include captured images, including the first and second pluralities of images captured when the measurement system and a calibration object are in the same media or in different media, extracted data associated with image points defining the images, known parameters such as refraction indices, calibration object parameters, interface parameters and the like. Some of such input data may be captured by the cameras 502 and 504 , or may be alternatively or additionally be pre-programmed into memory/media element 524 a . Still further, such information may also be entered as input data from a user accessing an input device 526 , which may correspond to one or more peripheral devices configured to operate as a user interface with computer control system 512 . Exemplary input devices may include but are not limited to a keyboard, touch-screen monitor, microphone, mouse and the like.

Second memory element 524 b includes computer-executable software instructions that can be read and executed by processor(s) 522 to act on the data stored in memory/media element 524 a to create new output data (e.g., processed image data, calibration constants, interface parameters, new error functions, reconstructed image points, etc.) for storage in a third memory/media element 524 c . Such output data may be provided to a peripheral output device 528 , such as monitor, printer or other device for visually depicting the output data, or as control signals to still further components. Computing/processing device(s) 522 may be adapted to operate as a special-purpose machine by executing the software instructions rendered in a computer-readable form stored in memory/media element 524 b . When software is used, any suitable programming, scripting, or other type of language or combinations of languages may be used to implement the teachings contained herein. In other embodiments, the methods disclosed herein may alternatively be implemented by hard-wired logic or other circuitry, including, but not limited to application-specific circuits.

Referring now to FIG. 7 , aspects of three-dimensional spatial position reconstruction are presented. Suppose I M1 me (X s , Y s ) and I M2 me (X s , Y s ) are the measured corresponding points on the stereo vision sensor planes that have been processed to remove distortion. The corresponding three-dimensional spatial point is Q(x w , y w , z w ). Two rays projected from I M1 me (X s , Y s ) and I M2 me (X s , Y s ) into three-dimensional space are not necessarily on the same plane, and hence, may not intersect due to errors in the calibration and measurement processes, as shown in FIG. 7 . In order to estimate the optimal three-dimensional spatial position Q(x w , y w , z w ), I M reco (X s , Y s ) and its corresponding model location I M mo — reco (X s , Y s ) are introduced. First I M reco (X s , Y s ) is the image point for Q(x w , y w , z w ) on the sensor plane. Second, I M reco (X s , Y s ) is assumed to be in close proximity to the measured location I M me (X s , Y s ). Third, the general least-square error function defined in equation (12) below is minimized to determine the optimal three-dimensional point Q(x w , y w , z w ). Specifics for I M mo — reco (X s , Y s ) are now described.

For object points in air, three-dimensional spatial position reconstruction without refraction is considered. I m mo — reco (X s , Y s ) is obtained by first assuming that Q(x w , y w , z w ) is unknown, and then transforming it from the WCS to the sensor plane by equation (7a). Therefore the reconstruction of Q without refraction is reduced to minimizing χ reco 2 .

For objects that are submerged in a fluid, three-dimensional spatial position reconstruction with refraction is considered, and I m mo — reco (X s , Y s ) is established by two steps. First, the refraction angle bta is introduced as an extra unknown and equations are set up for I M mo — reco (X s , Y s ) similar to those shown in modeling relationships set forth above. Second, a constraint for bta is set up as the sum of distances from B s to O″, from B x ″ to B s , and from Q″ to B x ″ as being equal to the distance from Q″ to O″, as shown in FIG. 8 .

It is noted that I m mo — reco (X s , Y s ) is different from I m mo (X s , Y s ) above in several ways. First, for I M mo — reco (X s , Y s ), the interface parameters a, b and d are known but for I m mo (X s , Y s ), they are not known. Second, there are no grid parameters for I M mo — reco (X s , Y s ), and each point on the object for reconstruction is treated individually. Third, the refraction angle bta for I M mo — reco (X s , Y s ) in FIG. 8 is physically the same as β for I M mo (X s , Y s ) in FIG. 5 , but bta is unknown and needs to be optimized while β is known and calculated based on I M me (X s , Y s ). Fourth, the unknown parameters for I M mo — reco (X s , Y s ) are Q(x w , y w , z w ), bta and λ (λ is used to define a constraint for bta), while those for I M mo (X s , Y s ) are interface parameters and grid parameters.

The least-square error function is set up as the sum of difference between I M me (X s , Y s ) and I M mo — reco (X s , Y s ) plus the constraint for the refraction angle bta as follows:

The unknown parameters for each camera, and λ for each bta. Reconstruction of Q(x w , y w , z w ) with refraction is reduced to minimizing χ reco 2 in equation (13).