System and method for 3D measurement and surface reconstruction

[0001] This continuation application claims priority to U.S. patent application Ser. No. 10/891,632 entitled A SYSTEM AND METHOD FOR 3D MEASUREMENT AND SURFACE RECONSTRUCTION, filed on Jul. 15, 2004 now abandoned, the contents of which are hereby incorporated by reference in their entirety.

[0002] The present invention relates to a system and method for 3D measurement and surface reconstruction of an image reconfigurable vision, and in particular to a reconfigurable vision system and method.

[0003] In many practical applications, such as reverse engineering, robotic exploration/navigation in clustered environments, model construction for virtual reality, human body measurements, and advanced product inspection and manipulation by robots, the automatic measurement and reconstruction of 3D shapes with high speed and accuracy is of critical importance. Currently, the devices widely used in industry for obtaining 3D measurements involve the mechanical scanning of a scene, for example in a laser scanning digitizer, which inevitably makes the measurement a slow process. Some advanced active vision systems using structured lighting have been explored and built. However, the existing systems lack the ability to change their settings, to calibrate by themselves and to reconstruct the 3D scene automatically.

[0004] To reconstruct a complete and accurate 3D model of an unknown object, two fundamental issues must be addressed. The first issue is how to acquire the 3D data for reconstructing the object surface. Currently, a laser range finder/scanner [1] is widely used for 3D surface data acquisition in industry. However, due to the mechanical scanning involved, the acquisition speed is limited. To increase the efficiency in the 3D imaging, pattern projections can be employed [2]. Portable 3D imaging systems based a similar principle have also been designed recently.

[0005] The second issue is how to determine the next viewpoint for each view so that all the information about the object surface can be acquired in an optimal way. This is also known as the NBV (Next Best View) problem, which determines the sensor direction (or pose) in the reconstruction process. The problem of viewpoint planning [3] for digitalization of 3D objects can be treated in different ways depending on whether or not the object's geometry is known beforehand [4,5]. For an unknown object, since the number of viewpoints and their viewing direction are unknown or cannot be determined prior to data acquisition, conventional 3D reconstruction processes typically involve an incremental iterative cycle of viewpoint planning, digitizing, registration and view integration and is conventionally based on a partial model reconstructed thus far. Based on a partial model reconstructed, the NBV algorithm then provides quantitative evaluations on the suitability of the remaining viewpoints. The evaluation for each viewpoint is based on all visible surface elements of the object that can be observed. The viewpoint with the highest visibility (evaluation score) is selected as the NBV.

[0006] In general, there are two fundamental problems to be solved when determining the Next Best View. The first problem is to determine the areas of the object which need to be sensed next and the second is to determine how to position the sensor to sample those areas. As there is no prior knowledge about the object, it is impossible to obtain a complete description of an object when occlusion occurs. Therefore, it is not generally possible to obtain precisely the invisible portions from either the current viewpoint or the acquired partial description of the object, so only an estimation of the Next Best View may be derived.

[0007] Various Next Best View algorithms have been proposed to date, for example Connolly [6] uses octree to represent object space, and the regions that have been scanned are labeled as seen, regions between the sensor and the surface are labeled as empty and all other regions are labeled as unseen. A set of candidate viewpoints is enumerated at fixed increments around the object. The Next Best View is calculated based on the evaluation of the visibility of each candidate viewpoint. This algorithm is computationally expensive and it does not incorporate the sensor geometry.

[0008] Maver and Bajesy [7] presented a solution to the NBV problem for a specific scanning setup consisting of an active optical range scanner and a turntable. In this document, unseen regions of the objects are represented as polygons. Visibility constraints for the sensor to view the unseen region are computed from the polygon boundaries. However, this solution is limited to a particular sensor configuration.

[0009] Pito [8] proposes an approach based on an intermediate position space representation of both sensor visibility constraints and unseen portions of the viewing volume. The NBV is determined as the sensor position that maximized the unseen portion of the object volume. This approach has been demonstrated to have achieved automatic viewpoint planning for a range sensor constrained to move on a cylindrical path around the object.

[0010] Whaite and Ferrie [9] use the superellipsoid model to represent an object and define a shell of uncertainty. The Next Best View is selected at the sensor position where the uncertainty of the current model fitted to the partial data points is the largest. This algorithm enables uncertainty-driven exploration of an object to build a model. However, the superellipsoid cannot accurately represent objects with a complex surface shape. Furthermore, surface visibility constraints were not incorporated in the viewpoint planning process.

[0011] Reed and Allen [10] propose a target-driven viewpoint planning method. The volume model is used to represent the object by extrusion and intersection operations. The constraints, such as sensor imaging constraints, model occlusion constraints and sensor placement constraints, are also represented as solid modeling volumes and are incorporated into the viewpoint planning. The algorithm involves expensive computation on the solid modeling and intersection operation.

[0012] Scott [11] considers viewpoint planning as integer programming. However, in this system the object must be scanned before viewpoint planning to obtain prior knowledge about an unknown object. Given a rough model of an unknown object, a sequential set of viewpoints is calculated to cover all surface patches of the object with registration constraint. However, the object must be scanned before viewpoint planning to obtain the prior knowledge about unknown objects.

[0013] In many applications, a vision sensor often needs to move from one place to another and change its configuration for perception of different object features. A dynamic reconfigurable vision sensor is useful in such applications to provide an active view of the features.

[0014] Active robot vision, in which a vision sensor can move from one place to another for performing a multi-view vision task, is an active research area. A traditional vision sensor with fixed structure is often inadequate for the robot to perceive the object's features in an uncertain environment as the object distance and size are unknown before the robot sees the object. A dynamically reconfigurable sensor may assist the robot in controlling the configuration and gaze at the object surfaces. For example, with a structured light system, the camera needs to see the object surface illuminated by the projector, to perform the 3D measurement and reconstruction task.

[0015] The system must be calibrated and traditionally, the calibration task is accomplished statically by manual operations. A calibration target/device is conventionally designed with a precision calibration fixture to provide a number of points whose world coordinates are precisely known [12]-[14]. With a planar calibration pattern, the target needs to be placed at several accurately known positions in front of the vision sensor. For dynamically reconfigurable vision systems, the vision system needs to have the ability of self-recalibration without requiring external 3D data provided by a precision calibration device.

[0016] Self-calibration of vision sensors has been actively researched in the last decade. However, most of the conventionally available methods were developed for calibration of passive vision systems such as stereo vision and depth-from-motion [15]-[22]. Conventionally these systems require dedicated devices for calibrating the intrinsic and extrinsic parameters of the cameras. Due to the special calibration target needed, such a calibration is normally carried out off-line before a task begins. In many practical applications, on-line calibration during the execution of a task is needed. Over the years, efforts have been made in research to achieve efficient on-line calibrations.

[0017] Maybank and Faugeras [23] suggested the calibration of a camera using image correspondences in a sequence of images from a moving camera. The kinds of constructions that could be achieved from a binocular stereo rig were further addressed in [24]. It was found that a unique projective representation of the scene up to an arbitrary projective transformation could be constructed if five arbitrary correspondences were chosen and an affine representation of the scene up to an arbitrary affine transformation could be constructed if four arbitrary correspondences were adopted.

[0018] Hartly [25] gave a practical algorithm for Euclidean reconstruction from several views with the same camera based on Levenberg-Marquardt Minimization. A new approach based on stratification was introduced in [26].

[0019] In this context, much work has been conducted in Euclidean reconstruction up to a transformation. Pollefeys et al [27] proposed a method to obtain a Euclidean reconstruction from images taken with an uncalibrated camera with variable focal lengths. This method is based on an assumption that although the focal length is varied, the principal point of the camera remains unchanged. This assumption limits the range of applications of this method. A similar assumption was also made in the investigations in [28,29]. In practice, when the focal length is changed (e.g. by zooming), the principal point may vary as well. In the work by Heyden and Astrom [30], they proved that it is possible to obtain Euclidean reconstruction up to a scale using an uncalibrated camera with known aspect ratio and skew parameters of the camera. A special case of a camera with Euclidean image plane was used for their study. A crucial step in the algorithm is the initialization which will affect the convergence. How to obtain a suitable initialization was still an issue to solve [31]. Kahl [32] presented an approach to self-calibration and Euclidean reconstruction of a scene, assuming an affine model with zero skew for the camera. Other parameters such as the intrinsic parameters could be unknown or varied. The reconstruction which needed a minimum of three images was an approximation and was up to a scale. Pollefeys et al gave the minimum number of images needed for achieving metric reconstruction, i.e. to restrict the projective ambiguity to a metric one according to the set of constraints available from each view [31].

[0020] The above-mentioned reconstruction methods are based on passive vision systems. As a result, they suffer from the ambiguity of correspondences between the camera images, which is a difficult problem to solve especially when free-form surfaces [33] are involved in the scene. However, to avoid this problem, active vision may be adopted. Structured light or pattern projection systems have been used for this purpose. To reconstruct precisely a 3D shape with such a system, the active vision system consisting of a projector and a camera needs to be carefully calibrated [34, 35]. The traditional calibration procedure normally involves two separate stages: camera calibration and projector calibration. These individual calibrations are carried out off-line and they have to be repeated each time the setting is changed. As a result, the applications of active vision systems are limited, since the system configuration and parameters must be kept unchanged during the entire measurement process.

[0021] For active vision systems using structured-light, the existing calibration methods are mostly based on static and manual operations. The available camera self-calibration methods cannot be applied directly to structured-light systems as they need more than two views for the calibration. Recently, there has been some work on self-calibration [36]-[40] of structured-light systems. Fofi et al. [36] investigated the self-calibration of structured-light systems, but the work was based on the assumption that a square projected onto a planar surface will most generally give a quadrilateral shape in the form of a parallelogram”.

[0022] Jokinen [37] studied a self-calibration method based on multiple views, where the object is moved by steps. Several maps were acquired for the registration and calibration. The limitation of this method is that the object must be placed on a special device so that it can be precisely moved.

[0023] Using a cube frame, Chu et al. [38] proposed a calibration free approach for recovering unified world coordinates.

[0024] Chen and Li [39, 40] recently proposed a self-recalibration method for a structured-light system allowing changes in the system configuration in two degrees of freedom.

[0025] In some applications, such as seabed metric reconstruction with an underwater robot, when the size or distance of the scene changes, the configuration and parameters of the vision system need to be changed to optimize the measurement. In such applications, uncalibrated reconstruction is needed. In this regard, efforts have been made in recent research. Fofi et al [41] studied the Euclidean reconstruction by means of an uncalibrated structured light system with a colour-coded grid pattern. They modeled the pattern projector as a pseudo camera and then the whole system as a two-camera system. Uncalibrated Euclidean reconstruction was performed with varying focus, zoom and aperture of the camera. The parameters of the structured light sensor were computed according to the stratified algorithm [26], [42]. However, it was not clear how many of the parameters of the camera and projector could be self-determined in the uncalibrated reconstruction process.

[0026] Thus, there is a need for a reconfigurable vision system and method for 3D measurement and reconstruction in which recalibration may be conducted without having to use special calibration apparatus as required by traditional calibration methods.

[0027] In general terms, the present invention provides a method and system for the measurement and surface reconstruction of a 3D image of an object comprising projecting a pattern onto the surface to be imaged, examining distortion produced in the pattern by the surface, converting for example by a triangulation process the distortions produced in the pattern by the surface to a distance representation representative of the shape of the surface. The surface shape of the 3D image may then be reconstructed, for example electronically such as digitally, for further processing.

[0028] In a preferred embodiment, the object is firstly sliced into a number of cross section curves, with each cross-section to be reconstructed by a closed B-spline curve. Then, a Bayesian information criterion (BIC) is applied for selecting the control point number of B-spline models. Based on the selected model, entropy is used as the measurement of uncertainly of the B-spline model to predict the information gain for each cross section curve. After obtaining the predicted information gain of all the B-spline models, the information gain of the B-spline models may be mapped into a view space. The viewpoint that contains maximal information gain for the object is then selected as the Next Best View. A 3D surface reconstruction may then be carried out.

[0029] An advantage of one or more preferred embodiments of the invention system is that the 3D information of a scene may be acquired at high speed by taking a single picture of the scene.

[0030] With this method, a complex 3D shape may be divided into a series of cross section curves, each of which represents the local geometrical feature of the object. These cross section curves may be described by a set of parametric equations. For reconstruction purposes using parametric equations, the most common methods include spline function (e.g. B-spline) [43], implicit polynomial [44], [45] and superquadric (e.g. superellipsoid). [46]. Compared with implicit polynomial and superquadric, B-spline has the following main advantages:

1. Smoothness and continuity, which allows a curve to consist of a concatenation of curve segments, yet be treated as a single unit;

2. Built-in boundedness, a property which is lacking in implicit or explicit polynomial representations whose zero set can shoot to infinity;

3. Parameterized representation. which decouples the x, y coordinates enabling them to be treated separately;

[0045] FIG. 1 shows an active vision system according to a preferred embodiment of the invention. The system comprises an LCD projector 1 adapted to cast a pattern of light onto an object 2 which is then viewed by a camera and processor unit 3. The relative position between the projector 1 and the camera in the camera and processing unit 3 has six degrees of freedom (DOF). When a beam of light is cast from the projector 1 and viewed obliquely by the camera, the distortions in the beam line may be translated into height variations via triangulation if the system is calibrated including the relative position between the projector 1 and camera. The vision system may be self-recalibrated automatically if and when this relative position is changed. The camera and processor unit 3 preferably includes a processor stage, as well as the camera, for processing the observed distortions in the projected pattern caused by the object 2 and associated data and for enabling and carrying out reconstruction. In a further preferred embodiment, the processor stage may be remotely located from the camera and may be connectable thereto to receive the data for processing and carrying out the reconstruction process.

[0046] FIG. 2 shows the geometrical relationship between the projector 1, the object 2 and the camera of the system of FIG. 1, and, in particular, the pattern projected by the projector 1 onto the object 2 and viewed by the camera 3.

[0047] For the camera,
x_c=P_cw_c,
where
x_c=[λx_c λy_c λ]⁷ are the coordinates on the image sensor plane,
λεR is an uncertain scalar
w_c=[X_c Y_c Z_c 1]^T are the 3D coordinates of an object point from the view of the camera (FIG. 2), and
P_c is the 3×4 perspective matrix

[0048] $P_c={\left[\begin{array}{cccc}{v}_c& k_{\mathrm{xy}}& x_c^0& 0\\ 0& s_{\mathrm{xy}}{v}_c& y_c^0& 0\\ 0& 0& 1& 0\end{array}\right]}_{3\times 4}.$
where
v_c is the distance between image plane and camera optical center,
s_xy is the ratio between the horizontal and vertical pixel cell sizes,
k_xy represents the placement perpendicularity of the cell grids,
and (x_c⁰,y_c⁰) is the center offset on the camera sensor.

[0049] Similarly, for the projector,
x_p=P_pw_p,
where
x_p=[κy_p κy_p κ]^T are the coordinates on the projector plane,
κεR is also an uncertain scalar,
w_p=[x_p Y_p Z_p 1]^T are the 3-D coordinates of the object point based on the view of projector (see FIG. 2), and
P_p is the inverse perspective matrix of the projector

[0050] $P_p={\left[\begin{array}{cccc}{v}_p& 0& x_p^0& 0\\ 0& v_p& y_p^0& 0\\ 0& 0& 1& 0\end{array}\right]}_{3\times 4}.$

[0051] The relationship between the camera coordinate system and projector coordinate system is

[0052] $w_p={\mathrm{Mw}}_c=\left[\begin{array}{cc}R& t\\ 0^T& 1\end{array}\right]w_c,$
where M is a 4×4 matrix and t is the translation vector,
t=s[t_xt_yt_z]^T.

[0053] Here s is a scaling factor to normalize t_x²+t_y²+t_z²=1.
x_p=P_pw_p=P_pMw_c.

[0054] Let

[0055] $H=P_{p}M={\left[\begin{array}{c}{h}_1\\ h_2\\ h_3\end{array}\right]}_{\left(3\times 4\right)},$
where h₁, h₂, and h₃ are 4-dimensional vectors. We have
κx_p=h₁w_c,κy_p=h₂w_c and κ=h₃w_c.
So
(x_ph₃−h₁)w_c=0.

[0056] Then the following can be derived:

[0057] $\left[\begin{array}{c}{P}_c\\ x_p{h}_3-h_1\end{array}\right]w_c=\left[\begin{array}{c}{x}_c\\ 0\end{array}\right],$

[0058] Denote x_c+=[x_c^T,0]^T and

[0059] $Q=\left[\begin{array}{c}{P}_c\\ x_p{h}_3-h_1\end{array}\right]=\left[\begin{array}{cccc}{q}_{11}& q_{12}& q_{13}& 0\\ 0& q_{22}& q_{23}& 0\\ 0& 0& 1& 0\\ q_{41}& q_{42}& q_{43}& q_{44}\end{array}\right].$

[0060] The 3-dimensional world position of a point on the object surface can be determined by
w_c=Q⁻¹x_c+.

[0061] As mentioned above, the relative positions of the camera 3 and the projector 1 may be changed dynamically during run-time of the system. As the camera (which is acting as a sensor) is reconfigurable during the run-time, it should be automatically recalibrated for 3-D perception tasks. Here the recalibration means that the camera (sensor) has been calibrated before installation in the system, but it should require calibrated again as the relative configuration changes. It is assumed for present purposes that the intrinsic parameters such as the focal lengths, scale factors, distortion coefficients will remain unchanged whereas the extrinsic parameters of the positions and orientations between the camera and projector have to be determined during the run-time of the system.

System Reconfiguration and Automatic Recalibration

[0062] The whole calibration of the structured light system of FIG. 1 may be divided into two parts. The first part concerns the calibration of the intrinsic parameters including the focal lengths and optical centers, this is called static calibration and may be performed off-line in a static manner. The second part deals with calibration of the extrinsic parameters of the relative position of the camera 3 and the projector 1, and this is hereinafter referred to as self-recalibration. The static calibration needs to be performed only once. The self-recalibration is thus more important and needs to be performed online whenever the system configuration is changed during a measurement task.

[0063] Once the perspective projection matrices of the camera 3 and the projector 1 relative to a global coordinate frame are computed from the static calibration, it is possible to obtain P_c and P_p which are the perspective projection matrices of the camera and the projector respectively relative to a global coordinate frame. The dynamic self-recalibration task requires the determination of the relative position M between the camera 3 and the projector 1. There are 6 unknown parameters, three for 3-axis rotation and three for 3-dimensional translation (as shown in FIG. 1).

[0064] For a point on an object surface, it is known that its coordinates on the camera's sensor plane x_c=[λx_c λy_c λ]^T and on the projector's source plane x_p=[κx_p κy_p κ]^T are related via the following:
x_p^TFx_c=0,
where F is a 3×3 essential matrix:

[0065] $F=\mathrm{sRS}=s\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],$

[0066] Here R is a 3-axis rotation matrix and S is a skew-symmetric matrix

[0067] $s=\left[\begin{array}{ccc}0& -t_z& t_y\\ t_z& 0& -t_x\\ -t_y& t_x& 0\end{array}\right],$
based on the translation vector t.

[0068] The recalibration task is to determine the 6 independent parameters in R and t (FIG. 2). For each surface point, x_p^TFx_c=0 may be expressed as:
a_i^Tf=0.
Here
f=[F₁₁F₂₁F₃₁F₁₂F₂₂F₃₂F₁₃F₂₃F₃₃]^T,
a_i=[x_cx_px_cy_px_cy_cx_py_cy_py_cx_py_p1]^T.
where (x_c, y_c) is the coordinates on the camera's sensor and (x_p, y_p) is the coordinates on the projector's LCD/DMD.

[0069] The projected patterns can be in black/white (b/w) or in colors. In either case, a coding method is in general needed. For b/w projections, gray codes can be used with the stripe light planes, which allows robust identifications of the stripe index (FIG. 3).

[0070] In a preferred embodiment of the invention, if an illumination pattern with colour-encoded grids is used, a cell's coordinates on the projector 1 can be immediately determined by the colours of adjacent neighbouring cells in addition to its own when projecting a source pattern. Via a table look-up, each cell's position can be uniquely identified. An example of such coded color pattern is shown in FIG. 4. The method of computing the values of blur diameters from an image has been proposed in [39,40] which is incorporated herein by reference. In a preferred embodiment, the system may comprise a Color Coded Pattern Projection system, a CCD camera, and a mini-platform for housing the components and providing the relative motion in 6DOF between the projector 1 and the camera 3. The method for automatic system recalibration and uncalibrated 3D reconstruction according to one or more preferred embodiments of the present invention may be implemented using such a system. With the adaptively adjustable sensor settings, this system will provide enhanced automation and performance for the measurement and surface reconstruction of 3D objects.

[0071] For n points observed, an n×9 matrix A can be obtained as the calibration data:
A=[a₁,a₂, . . . ,a_n]^T,
Af=0,

[0072] If it is assumed that the structured light vision system has 6DOF in its relative pose, i.e. three position parameters and three orientation parameters, between the camera 3 and the projector 1. The focal lengths of the projector 1 and the camera 3 are assumed to have been obtained in a previous static calibration stage. The optical centers are fixed or can be described by a function of the focal length. The projector 1 generates grid patterns with horizontal and vertical coordinates so that the projector's LCD/DMD can be considered an image of the scene. The relative position between the camera 3 and the projector 1 may be described by

[0073] $\left[\begin{array}{c}{X}_p\\ Y_p\\ Z_p\end{array}\right]=R\left[\begin{array}{c}{X}_c\\ Y_c\\ Z_c\end{array}\right]-\mathrm{Rt}.$

[0074] If there are n (n>5) points observed on a plane and they do not lie on the same line, we have proved that the rank of the calibration matrix A is six, i.e.
Rank(A)=6.

[0075] If the following 6-by-6 matrix is considered which is a sub-matrix of matrix A,
A₆=[r_a1,r_a2,r_a3,r_a4,r_a5,r_a6,]^T
where r_ai=[1x_ci y_ci x_pi y_pi x_ciy_ci,]^T and x_ci is the x value of the ith point projected on the camera coordinate system.

[0076] The matrix A₆ can be diagonalized by basic row-operations:
D(A₆)=diag(1,x_c,2−x_c,1, . . . )

[0077] Since x_c,i≠x_c,j, y_c,i≠y_c,j, x_p,i≠x_p,j, y_p,i≠y_p,j,

it can be proved that every element in D(A₆) is non-zero if no four points of the sampled data lie on the same line. Therefore det(A₆)≠0 and,
Rank(A)≧Rank(A₆)=6.

[0078] On the other hand, based on the projection model of the camera 3 and projector model 1, the coordinates of a surface point projected on the camera (sensor) may be given by (X/Z, Y/Z). For a point, (x_c,y_c) and (x_p,y_p) are related by:

[0079] $Z_p\left[\begin{array}{c}{x}_p\\ y_p\\ 1\end{array}\right]=Z_{c}R\left[\begin{array}{c}{x}_c\\ y_c\\ 1\end{array}\right]-\mathrm{Rt}.$
where Z_c and Z_p are the depth values based on the view of the camera 3 and projector 1, respectively.

[0080] For the camera 3 and projector 1, the scene plane may be defined as follows:

[0081] $Z=\frac{C_3}{1-C_{1}x-C_{2}y}$

[0082] Let r₁, r₂, and r₃ be the three rows in R, and then

[0083] $\frac{C_{3c}}{1-C_{1c}{x}_p-C_{2c}{y}_p}\left[\begin{array}{c}{x}_p\\ y_p\\ 1\end{array}\right]=\frac{C_{3c}}{1-C_{1c}{x}_c-C_{2c}{y}_c}\left[\begin{array}{c}{r}_1\\ r_2\\ r_3\end{array}\right]\left[\begin{array}{c}{x}_c\\ y_c\\ 1\end{array}\right]-\left[\begin{array}{c}{r}_1\\ r_2\\ r_3\end{array}\right]t,$
which contains three equations and then, this equation may be considered to be equivalent to the following system:

[0084] $\{\begin{array}{c}{\tau }_{11}{x}_c{x}_p+{\tau }_{12}{x}_c{y}_p+{\tau }_{13}{x}_p{y}_c+{\tau }_{14}{y}_c{y}_p+\\ {\tau }_{15}{x}_c+{\tau }_{16}{x}_p+{\tau }_{17}{y}_c+{\tau }_{18}{y}_p+{\tau }_{19}=0\\ {\tau }_{21}{x}_c{x}_p+{\tau }_{22}{x}_c{y}_p+{\tau }_{23}{x}_p{y}_c+{\tau }_{24}{y}_c{y}_p+\\ {\tau }_{25}{x}_c+{\tau }_{26}{x}_p+{\tau }_{27}{y}_c+{\tau }_{28}{y}_p+{\tau }_{29}=0\\ {\tau }_{31}{x}_c{x}_p+{\tau }_{32}{x}_c{y}_p+{\tau }_{33}{x}_p{y}_c+{\tau }_{34}{y}_c{y}_p+\\ {\tau }_{35}{x}_c+{\tau }_{36}{x}_p+{\tau }_{37}{y}_c+{\tau }_{38}{y}_p+{\tau }_{39}=0\end{array}\mathrm{or}\Gamma a=0,$
where Γ is a 3-by-9 matrix, a is described above, and {τ_ij} are constants.

[0085] It can be proved that there is no linear relationship among the above three equations, i.e. rank(Γ)=3.

[0086] Considering 9 points as the calibration data, then matrix A is 9-by-9 in size. Since it is constrained by Γa=0, the maximum rank of A is 6.

[0087] Therefore the rank of matrix A must be 6.

[0088] The general solution of the equation a_i^Tf=0 has the form of
f=ξ₁f₁+ξ₁f₂+ξ₁f₃,
where ξ₁, ξ₂, and ξ₃ are real numbers, f_i is a 9-dimensional vector, and [f₁; f₂; f₃] is the null-basis of A.

[0089] Using singular value decomposition (SVD), we have
B=svd(A^TA)=UDV^T,
where A^TA is a 9×9 matrix, D is a non-decreasing diagonal matrix, and U and V are orthogonal matrices.

[0090] Then, f₁, f₂, and f₃ are the three vectors in V corresponding to the least eigenvalues. Theoretically, if there is no noise, matrix B is exactly as just described, i.e. of rank 6. In such a case, there would be three vanishing/singular values in the diagonal matrix D and the sum and or mean of squared errors (SSE and/or MSE) would be zero, since the vector f lies exactly in the null-space of B.

[0091] However, in a practical system, there maybe fewer than 3 singular values in the matrix D, as the matrix B can be perturbed. Since the data are from real measurements, B may have a rank of 7 or even higher. In such a case it is still possible to take the three column vectors from the matrix V as the basis vectors corresponding to the least values in D. It will still be the best in the sense B will have been tailored (with rank B>6) to some other matrix C with rank C=6 in such a way that C is the “nearest” to B among all the n-by-9 matrices with rank 6 in terms of the spectral norm and Frobenius norm.

[0092] Define
f=Hk,
where k=[ξ₁ ξ₂ ξ₃]^T and H=[f₁ f₂ f₃]=[H_u^T, H_m^T, H_l^T]^T. H_u, H_m, and H_l are 3×3 matrices and each for three rows in H.

[0093] The above can be written as
H_uk=F_c1,H_mk=F_c2, and H_lk=F_c3,
where F_c1, F_c2, and F_c3 are the three columns in F. Therefore,

[0094] $G=F^{T}F=\left[\begin{array}{ccc}{k}^T{H}_u^T{H}_{u}k& k^T{H}_u^T{H}_{m}k& k^T{H}_u^T{H}_{l}k\\ k^T{H}_m^T{H}_{u}k& k^T{H}_m^T{H}_{m}k& k^T{H}_m^T{H}_{l}k\\ k^T{H}_l^T{H}_{u}k& k^T{H}_l^T{H}_{m}k& k^T{H}_l^T{H}_{l}k\end{array}\right]$

[0095] As R is orthogonal, F^TF can also be expressed as

[0096] $G=S^T{R}^T\mathrm{RS}=S^{T}S=\left[\begin{array}{ccc}1-t_x^2& -t_x{t}_y& -t_x{t}_z\\ -t_x{t}_y& 1-t_y^2& -t_y{t}_z\\ -t_x{t}_z& -t_y{t}_z& 1-t_z^2\end{array}\right].$

[0097] The three unknowns of k=[ξ₁ ξ₂ ξ₃]^T can be determined. The normalized relative position t_n=[t_x t_y t_z]^T can then be solved:

[0098] $\hspace{1em}\{\begin{array}{c}{t}_x=\pm \sqrt{1-k^T{H}_u^T{H}_{u}k}\\ t_y=\pm \sqrt{1-k^T{H}_m^T{H}_{m}k}\\ t_z=\pm \sqrt{1-k^T{H}_l^T{H}_{l}k}\end{array}$

[0099] It should be noted that multiple solutions exist. In fact, if [k t]^T is a solution of the system, [±k ±t]^T must also be the solutions. One of these solutions is correct for a real system setup. To find this, the re-projection method can be used.

[0100] When k and t are known, the rotation matrix R can be determined by

[0101] $R=\left[\begin{array}{c}{r}_1\\ r_2\\ r_3\end{array}\right]=\left[\begin{array}{c}{F}_{c1}\times t+\left(F_{c2}\times t\right)\times \left(F_{c3}\times t\right)\\ F_{c2}\times t+\left(F_{c3}\times t\right)\times \left(F_{c1}\times t\right)\\ F_{c3}\times t+\left(F_{c1}\times t\right)\times \left(F_{c2}\times t\right)\end{array}\right]$
where “x” is the cross product of two vectors.

[0102] Among the six unknown parameters, the five in R and t_n have been determined so far and reconstruction can be performed but up to a scaling factor. The last unknown, s in t, may be determined here by a method using a constraint derived from the best-focused location (BFL). This is based on the fact that for a lens with a specified focal length, the object surface point is perfectly focused only at a special distance.

[0103] For an imaging system, the mathematical model for standard linear degradation caused by blurring and additive noise is usually described by

[0104] $g\left(i,j\right)=\sum _{k=1}^m\sum _{l=1}^{n}h\left(i-k,j-l\right)f\left(k,l\right)+n\left(i,j\right),\mathrm{or}$ $g=h\otimes f+n$
where f is the original image, h is the point spread function, n is the additive noise, m×n is the image size. The operation represents two-dimensional convolution. The blur can be used as a cue to find the perfectly focused distance.

[0105] For the projector in such a system, the most significant blur is the out-of-focus blur. This results from the fact that for a lens with a specific focal length the illumination pattern will be blurred on the object surface unless it is projected on the perfectly focused distance. Since the noise n only affects the accuracy of the result, it will not be considered in the following deduction. The illumination pattern on the source plane (LCD or DMD) to be projected is described as

[0106] $I_s\left(x\right)=\{\begin{array}{cc}{L}_a,& \left(-\frac{T}{2}<x-2\mathrm{nT}<\frac{T}{2}\right),n\in N\\ 0,& \mathrm{otherwise}\end{array},$
where T is the stripe width of the source pattern.

[0107] The scene irradiance caused by a light source is inversely proportional to the square of the distance from the light source. On the other hand, the image intensity on the camera sensor is independent of the scene distance and is only proportional to the scene irradiance. Therefore, the image intensity can be described as

[0108] $l=\frac{C_c{C}_I}{l^2},$
where C_c is the sensing constant of the camera and C_l is the irradiant constant of light projection. They are related to many factors, such as the diameter of sensor's aperture, the focal length, and properties of surface materials.

[0109] Assume that the intensity at the point where the best-focused plane intersects the principal axis of the lens is l₀ (x=0, z=z₀), where z₀ is the best-focused distance. That is

[0110] $l_0=\frac{C_c{C}_I}{{\left(z_0+v_p\right)}^2},$
where v_p is the focal length, i.e. the distance between the source plane and optical center of the projector lens.

[0111] Consider a straight-line in the scene projected on the X-Z coordinate system:
z=c₁x+c₀.

[0112] For an arbitrary point, we have l=√{square root over (x²+(z+v_p)²)} and thus

[0113] $I=\frac{{\left(z_0+v_p\right)}^2{I}_0}{l^2}.$

[0114] In the view of the projector, when the illumination pattern casts on a line, the intensity distribution becomes nonlinear and is given by

[0115] $\begin{array}{cc}{I}_i\left(x\right)=\frac{{\left(z_0+v_p\right)}^2}{{\left(z+v_p\right)}^2+{x^2\left(1+\frac{v_p}{z}\right)}^2}{I}_0,\left(-\frac{T}{2}<\frac{v_p}{z}x-2\mathrm{nT}<\frac{T}{2}\right),n\in N.& \left(8\right)\end{array}$

[0116] Transforming the x-axis to align with the observed line, we have:
x_l=√{square root over (1+c₁²)}x.

[0117] The above gives

[0118] $I_i\left(x_I\right)=\frac{{\left(z_0+v_p\right)}^2}{{\left[I+{\left(\frac{c_2{x}_l}{c_2{x}_l+c_0}\right)}^2\right]\left[\left(c_2+I\right)x_l+c_0\right]}^2}{I}_0,\mathrm{where}{c}_2=\frac{c_1}{\sqrt{1+c_1^2}}.$

[0119] The illumination will be blurred unless it is projected on a plane at the perfectly focused distance:

[0120] $z_0=\frac{v_p{f}_p}{v_p-f_p},$
where f_p is the intrinsic focus length of the projector. For all other points in the scene, the z-displacement of a surface point to the best-focused location is:

[0121] $\Delta z=z-z_0=z-\frac{v_p{f}_p}{v_p-f_p}.$
where v_p is the distance to from the image plane to the optical center.

[0122] The corresponding blur radius is proportional to Δz:

[0123] $\sigma =\frac{v_p-f_p}{v_p{F}_{\mathrm{num}}}\Delta z,\mathrm{where}{F}_{\mathrm{num}}=\frac{f_p}{r}$
is the f-number of the lens setting.

[0124] For out-of-focus blur, the effect of blurring can be described via a point spread function (PSF) to account for the diffraction effect of light wave. A Gaussian model is normally used.

[0125] With our light stripes, the one-dimensional PSF is

[0126] $h_{\sigma }\left(x\right)=\frac{1}{\sqrt{2\pi }\sigma }{e}^{-\frac{x^2}{2{\sigma }^2}}.$

[0127] The brightness of the illumination from the projector is the convolution of the ideal illumination intensity curve with the PSF blur model:

[0128] $I\left(x\right)=I_i\left(x\right)\otimes h_{\sigma }\left(x\right)={\int }_{-\infty }^{+\infty }{I}_i\left(u\right)h_{\sigma }\left(x-u\right)du.$

[0129] The Fourier transform of the above is
I_F(ω)=I_Fⁱ(ω)H_σ(ω),  (18)
where H_σ(ω) is the Fourier transform of the Gaussian function

[0130] $H_{\sigma }\left(\omega \right)={\int }_{-\infty }^{+\infty }\frac{1}{\sqrt{2\pi }\sigma }{e}^{-\frac{x^2}{2{\sigma }^2}}{e}^{-\mathrm{j\omega }x}dx=e^{-\frac{{\sigma }^2{\omega }^2}{2}}.$

[0131] Without significant loss of accuracy, l_i(x) may be approximated by averaging the intensity on a light stripe to simplify the Fourier transform, I_i(x)=Ī(x). If a coordinate system with its origin at the center of the bright stripe is used, this intensity can be written as

[0132] $\stackrel{\_}{I}\left(x\right)=I_0\left[\varepsilon \left(x+\frac{T}{2}\right)-\varepsilon \left(x-\frac{T}{2}\right)\right],$
where ε is a unit step function.

[0133] The Fourier transform of the above is

[0134] $I_F^i\left(\omega \right)=I_{0}T\frac{\sin \left(\frac{\omega T}{2}\right)}{\frac{\omega T}{2}}=I_0{\mathrm{TS}}_a\left(\frac{\omega T}{2}\right).$

[0135] Since l(x) is measured by the camera, its Fourier transform l_F(ω) can be calculated. Using integration, we have

[0136] $\int e^{-\frac{{\sigma }^2{\omega }^2}{2}}d\omega =I_{0}T\int \frac{I_F\left(\omega \right)}{S_a\left(\frac{\omega T}{2}\right)}d\omega .$

[0137] The left side is found to be

[0138] $\mathrm{Left}=\frac{\sqrt{2}}{\sigma }{\int }_{-\infty }^{+\infty }{e}^{-{\left(\frac{\sigma \omega }{\sqrt{2}}\right)}^2}d\left(\frac{\mathrm{\sigma \omega }}{\sqrt{2}}\right)=\frac{\sqrt{2}\sqrt{\pi }}{\sigma }.$

[0139] Therefore the blur radius can be computed by

[0140] $\sigma ={\left[\frac{I_{0}T}{\sqrt{2\pi }}\int \frac{I_F\left(\omega \right)}{S_a\left(\frac{\omega T}{2}\right)}d\omega \right]}^{-1}.$

[0141] Neglecting the effect of blurring caused by multiple illumination stripes, we have the following theorem to determine the blur radius with low computational cost and high precision.

[0142] Theorem 1. With the projection of a step illumination on the object surface, the blur radius is proportional to the time rate flow of irradiant light energy in the blurred area:

[0143] $\sigma =\frac{\sqrt{2\pi }}{I_0}S,$
where l₀ is the ideal intensity and S is the area size as illustrated in FIG. 5b.

[0144] This means that the blur radius σ[m] is proportional to the area size under the blurring curve: σ=√{square root over (2π)}S/I₀. The time rate flow of radiant light energy, i.e. the irradiant power or irradiant flux Φ [watt], is also the area size S [watt] under the blurring curve (or surface, in the case of two-dimensional analysis) illustrated in FIG. 5b.

[0145] Therefore, we only need to compute the area size S for every stripe to determine the blur radius. In a simple way, the edge position (x=0) can be detected by a gradient method, and S is then determined by summating the intensity function from 0 to x₁. However, even using a sub-pixel method for the edge detection, errors are still considerable since l(x) changes sharply near the origin.

[0146] To solve this problem, we propose an accurate method in the sense of energy minimization. As illustrated in FIG. 6, we have

[0147] $F_s\left(x_o\right)=S_1\left(x_o\right)+S_2\left(x_o\right)=\frac{I_0}{\sqrt{\pi }}{\int }_{\frac{x_o}{\sqrt{2}\sigma }}^{+\infty }{e}^{-y^2}dy+I_0-\frac{I_0}{\sqrt{\pi }}{\int }_{-\infty }^{\frac{x_o}{\sqrt{2}\sigma }}{e}^{-y^2}dy.$

[0148] It can be proved that the derivative of the above function is
F_s′≧0, if x_o≧0,
where “=” holds if and only if x_o=0.

[0149] The same situation occurs when x_o≦0. Therefore, at x_o=0, we have
S=min(F_s)/2.

[0150] This means that the same quantity of light energy flows from S₂ to S₁. This method for computing S is more stable than traditional methods and it yields high accuracy.

[0151] Now the best-focused location can be computed by analyzing the blur information in an image. With Theorem 1, by integrating the blurring curve on each stripe edge, the blur radius σ_i can be calculated. These blur radiuses are recorded as a set
D={(i,σ_i)|iεN},
where i is the stripe index on the projector's source pattern.

[0152] The blur size is proportional to the displacement of a scene point from the BFL, σ=k_fΔz. Since the blur diameters are unsigned, a minimum value σ_min in the data set D can be found. For a line in the scene, in order to obtain a straight line corresponding to the linearly changing depth in the scene, we separate D into two parts and apply linear best-fit to obtain two straight lines:
σ_l(x)=k₁x+k₂ and σ_r(x)=k₃x+k₄.

[0153] Finding the intersection of the two lines gives a best-focused location (as shown in FIG. 7)

[0154] $x_b=\frac{k_4-k_2}{k_1-k_3},$
which corresponds to Δz=0 or

[0155] $z\left(x_b\right)=\frac{v_p{f}_p}{v_p-f_p}.$
The corresponding coordinates on the image are (x_b,y_b), where y_b is determined according to the scanning line on the image.

[0156] From the above analysis, there exists a point which represents the best-focused location for a line in the scene. Now consider the blur distribution on a plane in the scene. This happens when we analyze multiple scan lines crossing the light stripes or when the source illumination is a grid pattern.

[0157] For a plane in scene, z=c₁x+c₂y+c₃, the blur radius is

[0158] $\begin{array}{c}\sigma \left(x,y\right)=\frac{v_p-f_p}{v_p{F}_{\mathrm{num}}}z-\frac{v_p{f}_p}{v_p-f_p}\\ =\frac{\left(v_p-f_p\right)\left(c_{1}x+c_{2}y+c_3\right)-v_p{f}_p}{v_p{F}_{\mathrm{num}}}.\end{array}$

[0159] The best focused locations form a straight line which is the intersection of two planes. A valley line can be found since the blur radius is unsigned.

[0160] For a freeform surface in the scene, the best-focused location can be determined by extending the above method making some minor modifications. For each light stripe, we can also compute its blur diameter and obtain a pair (Z_ri, σ_i), where iεN is the stripe index and z_ri iε[0, 1] is its relative depth in the camera view. Plotting these pairs in a coordinate system with σ as the vertical axis and z_r as the horizontal axis, we can also find a valley to be the best focused distance, z_rb.

[0161] For the point with minimum blur value, i.e. the best-focused location (Δz=0), is constrained by

[0162] $Z_c\left[\begin{array}{c}{x}_c\\ y_c\\ 1\end{array}\right]-s\left[\begin{array}{c}{t}_x\\ t_y\\ t_z\end{array}\right]=\frac{v_p{f}_p}{v_p-f_p}{R}^{-1}\left[\begin{array}{c}{x}_p\\ y_p\\ 1\end{array}\right]=\left[\begin{array}{c}{b}_1\\ b_2\\ b_3\end{array}\right],$

[0163] The scaling factor s can thus be determined:

[0164] $s=\frac{x_c{b}_2-y_c{b}_1}{y_c{t}_x-x_c{t}_y}.$

[0165] The procedures for self-recalibration of a structured light vision sensor are given as follows:

• • [0166] Step 1: projecting grid encoded patterns onto the scene;
  • [0167] Step 2: determining t_n and R for 5 unknowns;
  • [0168] Step 3: computing the blur distribution and determining the best-focused location;
  • [0169] Step 4: determining the scaling factor s; and
  • [0170] Step 5: combining the relative matrix for 3-D reconstruction.

[0171] The method presented here automatically resolves six parameters of a color-encoded structured light system. When the internal parameters of the camera 3 and the projector 1 are known via pre-calibration, the 6 DOF relative placement of the active vision system can be automatically recalibrated with neither manual operations nor assistance of a calibration device. This feature is very important for many situations when the vision sensor needs to be reconfigured online during a measurement or reconstruction task.

[0172] The method itself does not require the six parameters all to be the external ones. In fact, only if the number of unknown parameters of the system does not exceed six, the method can still be applied. For example, if the two focal lengths, v_c and v_p, are variable during the system reconfiguration and the relative structure has 4DOFs, we may solve them in a similar way, by replacing the matrix F by F=(P_c⁻¹)^TRSP_p⁻¹ and modifying the decomposition method accordingly.

[0173] When the unknown parameters exceed six, the above-described method may not solve them directly. However, a 2-step method may be used to solve this problem. That is, before the internal parameters are reconfigured, we take an image firstly to obtain the 6DOF external parameters. Then after changing the sensors' internal parameters, we take an image again to recalibrate them.

[0174] Special features of a method according to a preferred embodiment of the invention the invented method include:

• • [0175] The single image based recalibration allows measurement or reconstruction to be performed immediately after reconfiguration in the software, without requiring any extra requirement.
  • [0176] Metric measurement of the absolute geometry of the 3D shape may be obtained by replacing r(x_b, y_b) with z_rb. This is different from most of the currently available conventional methods where the 3D reconstruction supported is up to a certain transformation.
    Automatic Viewpoint Planning

[0177] To obtain a complete 3D model of an object 2 with the vision system as shown in FIG. 1, multiple views may be required, e.g. via robot vision as shown in FIG. 8 which illustrates the components of FIG. 1 housed in a robot apparatus. The viewpoint planning is charged with the task of determining the position and orientation and system configuration parameters for each view to be taken. It is assumed for the purposes of this description that we are dealing with an unknown object, i.e. assuming no prior knowledge about the object model. It is also assumed here that the object is general free-form surfaced. The approach is preferably to model a 3D object via a series of cross section curves. These cross section curves can be described by a set of parametric equations of B-spline curves. A criterion is proposed to select the optimal model structure for the available data points on the cross section curve.

[0178] For object reconstructions, two conventional approaches are volume reconstruction and surface reconstruction. The volume based technique is concerned with the manipulation of the volumetric objects stored in a volume raster of voxels. Surface reconstruction may be approached in one of the two ways: 1) representing the surface with a mosaic of flat polygon titles, usually triangles: and 2) representing the surface with a series of curved patches joined with some order of continuity. However, as mentioned above, a preferred reconstruction method for embodiments of the present invention is to model a 3D object via a series of cross section curves.

Model Selection for Reconstruction

[0179] The object 2 is sliced into a number of cross section curves, each of which represents the local geometrical features of the object. These cross section curves may be described by a set of parametric equations. For reconstruction of cross section curves, compared with implicit polynomial [47] and superquadric, B-spline has the following main advantages:

1) Smoothness and continuity, which allows any curve to be composed of a concatenation of curve segments and yet be treated as a single unit,

2) Built-in boundedness. a property which is lacking in implicit or explicit polynomial representation whose zero set can shoot to infinity; and

3) Parameterized representation, which decouples the x, z coordinates to be treated separately.

2.1 Closed B-Spline Curve Approximation

[0180] A closed cubic B-spline curve consisting of n+1 curve segments may be defined by

[0181] $\begin{array}{cc}p\left(t\right)=\sum _{j=0}^{n+3}{B}_{j,4}\left(t\right)·{\Phi }_j& \left(1\right)\end{array}$
where p(t)=[x(t), z(t)] is a point on B-spline curve with location parameter t. In this section we use the chord length method for parameterization. In (1), B_j,4(t) is the jth normalized cubic B-spline basis function, which is defined over the uniform knots vector
[u₋₃,u₋₂, . . . ,U_n+4,]=[−3,−2, . . . ,n+4]

[0182] In addition, the amplitude of B_j,4(t) is in the range of (0.0, 1.0), and the support region of B_j,4(t) is compact and nonzero for tε[u_j, u_j+4]. The (Φ_j)_j=0ⁿ⁺³ are cyclical control points satisfying the following conditions:
Φ_n+1=Φ₀,Φ_n+2=Φ₁,Φ_n+3=Φ₂

[0183] By factorization of the B-spline model, the parameters of B-spline model can be represented as:
[Φ_x^τ,Φ_z^τ]^τ=[Φ_x0, . . . Φ_xn,Φ_z0, . . . Φ_zn]^τ

[0184] For a set of m data points r=(r_i)_i=1^m=([x_i, y_i])_i=1^m, let d² be the sum of the squared residual errors between the data points and their corresponding points on the B-spline curve, ie

[0185] $\begin{array}{c}{d}^2=\sum _{i=1}^m{r_i-p\left(t_i\right)}^2\\ =\sum _{i=1}^m{\left[x_i-\sum _{j=0}^{n+3}{B}_{j,4}\left(t_i\right)·{\Phi }_{\mathrm{xj}}\right]}^2+\sum _{i=1}^m{\left[y_i-\sum _{j=0}^{n+3}{B}_{j,4}\left(t_i\right)·{\Phi }_{\mathrm{yj}}\right]}^2\end{array}$

[0186] From the cyclical condition of control points in the equation Φ_n+1=Φ₀, Φ_n+2=Φ₁, Φ_n+3=Φ₂, there are only n+1 control points to be estimated. The LS estimation method of the n+1 control points may be obtained from the curve points by minimizing d² above with respect to Φ=[Φ_x^T, Φ_y^T]^T=[Φ_x0, . . . Φ_xn, Φ_y0, . . . , Φ_yn]^T.

[0187] The following estimation of Φ may then be obtained by factorization of the B-spline:

[0188] $\begin{array}{cc}\begin{array}{c}{\Phi }_x={\left[B^{\tau }B\right]}^{-1}{B}^{\tau }x\\ {\Phi }_y={\left[B^{\tau }B\right]}^{-1}{B}^{\tau }y\end{array}\}\mathrm{where}x={\left[x_1,\dots ,x_m\right]}^{\tau },y={\left[y_1,\dots ,y_m\right]}^{\tau }B=\left[\begin{array}{cccc}{\stackrel{\_}{B}}_{0,4}^1+{\stackrel{\_}{B}}_{n+1,4}^1& {\stackrel{\_}{B}}_{1,4}^1+{\stackrel{\_}{B}}_{n+2,4}^1& {\stackrel{\_}{B}}_{2,4}^1+{\stackrel{\_}{B}}_{n+3,4}^1\dots & {\stackrel{\_}{B}}_{n,4}^1\\ {\stackrel{\_}{B}}_{0,4}^2+{\stackrel{\_}{B}}_{n+1,4}^2& {\stackrel{\_}{B}}_{1,4}^2+{\stackrel{\_}{B}}_{n+2,4}^2& {\stackrel{\_}{B}}_{2,4}^2+{\stackrel{\_}{B}}_{n+3,4}^2\dots & {\stackrel{\_}{B}}_{n,4}^2\\ ⋮& ⋮& ⋮& ⋮\\ {\stackrel{\_}{B}}_{0,4}^m+{\stackrel{\_}{B}}_{n+1,4}^m& {\stackrel{\_}{B}}_{1,4}^m+{\stackrel{\_}{B}}_{n+2,4}^m& {\stackrel{\_}{B}}_{2,4}^m+{\stackrel{\_}{B}}_{n+3,4}^m\dots & {\stackrel{\_}{B}}_{n,4}^m\end{array}\right]& \left(2\right)\end{array}$
and B_j,4ⁱ=B_j,4(τ_i).

[0189] The chord length method may preferably be used for the parameterization of the B-spline. The chord length L of a curve may be calculated as follows:

[0190] $L=\sum _{i=2}^{m+1}r_i-r_{i-1}$
where r_m+1=r₁ for a closed curve. The t_i associated with point q_i may be given as:

[0191] $t_i=t_{i-1}+\frac{r_i-r_{i-1}}{L}·t_{\max }$
where t₁=0 and t_max=n+l
Model Selection with Improved BIC Criterion

[0192] For a given set of measurement data, there exists a model of optimal complexity corresponding to the smallest prediction (generalization) error for further data. The complexity of a B-spline model of a surface is related to its control point (parameter) number [43],[48]. If the B-spline model is too complicated, the approximated B-spline surface tends to over-fit noisy measurement data. If the model is too simple, then it is not capable of fitting the measurement data, making the approximation results under-fitted. In general, both over- and under-fitted surfaces have poor generalization capability. Therefore, the problem of finding an appropriate model, referred to as model selection, is important for achieving a high level generalization capability.

[0193] Model selection has been studied from various standpoints in the field of statistics. Examples include information statistics [49]-[51] Bayesian statistics [52]-[54], and structural risk minimization [55]. The Bayesian approach is a preferred model selection method. Based on posterior model probabilities, the Bayesian approach estimates a probability distribution over an ensemble of models. The prediction is accomplished by averaging over the ensemble of models. Accordingly, the uncertainty of the models is taken into account, and complex models with more degrees of freedom are penalized.

[0194] For a given set of models {M_k, k=1,2, . . . } and data r, there exists a model of optimal model structure corresponding to smallest generalization error for further data and the Bayesian approach may be used to select the model with the largest (maximum) posterior probability to account for the data acquired so far.

[0195] In a first preferred method, the model M may be denoted by:
M=arg_{Mk,k=1 . . . k max}max{p(r|M_K)}
where the posterior probability of model M_k may be denoted by

[0196] $\begin{array}{c}p\left(r❘M_k\right)={\int }_{{\Phi }_k}p\left(r❘{\Phi }_k,M_k\right)p\left({\Phi }_k❘M_k\right)d{\Phi }_k\\ \cong {\left(2\pi \right)}^{d_k/2}{H\left({\hat{\Phi }}_k\right)}^{-1/2}p\left(r❘{\hat{\Phi }}_k,M_k\right)p\left({\hat{\Phi }}_k❘M_k\right)\end{array}$

[0197] Neglecting the term p({circumflex over (Φ)}_k|M_k), the posterior probability of model M_k becomes [11]:

[0198] $M=\arg \underset{M_k,k=1,\dots k_{\max }}{\max }\left\{\log p\left(r❘{\hat{\Phi }}_k,M_k\right)-\frac{1}{2}\log H\left({\hat{\Phi }}_k^{\prime }\right)\right\}$

• • [0199] where {circumflex over (Φ)}_k is the maximum likelihood estimate of Φ_k, and d_k is the parameter number of model M_k, H({circumflex over (Φ)}_k)
  • [0200] is the Hessian matrix of −log p(r|Φ_k,M_k) evaluated at {circumflex over (Φ)}_k.

[0201] The likelihood function p(r|{circumflex over (Φ)}_k,M_k) of closed B-spline cross section curves can be factored into x and y components as
p(r|{circumflex over (Φ)}_k,M_k)=p(x|{circumflex over (Φ)}_kz,M_k)·p(y|{circumflex over (Φ)}_ky,M_k)
where {circumflex over (Φ)}_kz and {circumflex over (Φ)}_ky can be calculated by

[0202] $\{\begin{array}{c}{\Phi }_x={\left[B^{T}B\right]}^{-1}{B}^{T}x\\ {\Phi }_y={\left[B^{T}B\right]}^{-1}{B}^{T}y\end{array}$

[0203] Consider, for example, the x component. Assuming the residual error sequence to be zero mean and white Gaussian with variance σ_kx²({circumflex over (Φ)}_kx), we have the following likelihood function:

[0204] $p\left(x❘{\hat{\Phi }}_{\mathrm{kx}},M_k\right)={\left(\frac{1}{2\pi {\sigma }_{\mathrm{kx}}^2\left({\hat{\Phi }}_{\mathrm{kx}}\right)}\right)}^{m/2}·\exp \left\{-\frac{1}{2{\sigma }_{\mathrm{kx}}^2\left({\hat{\Phi }}_{\mathrm{kx}}\right)}\sum _{k=0}^{m-1}{\left[x_k-B_k{\hat{\Phi }}_{\mathrm{kx}}\right]}^2\right\}$

• • [0205] and σ_kx²({circumflex over (Φ)}_kx,M_k) is estimated by:

[0206] ${\hat{\sigma }}_{\mathrm{kx}}^2\left({\hat{\Phi }}_{\mathrm{kx}}\right)=\frac{1}{m}\sum _{k=0}^{m-1}{\left[x_k-B_k{\hat{\Phi }}_{\mathrm{kx}}\right]}^2$

[0207] In a similar way, the likelihood function of the y component can also be obtained. The corresponding Hessian matrix insert Ĥ_k of −log p(r|Φ_k,M_k)

• • [0208] evaluated at {circumflex over (Φ)}_k may be denoted by:

[0209] $H\left({\hat{\Phi }}_k\right)=\left[\begin{array}{cc}\frac{B^{T}B}{{\hat{\sigma }}_{\mathrm{kx}}^2\left({\hat{\Phi }}_{\mathrm{kx}}\right)}& 0\\ 0& \frac{B^{T}B}{{\hat{\sigma }}_{\mathrm{ky}}^2\left({\hat{\Phi }}_{\mathrm{ky}}\right)}\end{array}\right]$

[0210] By approximating

[0211] $\frac{1}{2}\log H\left({\hat{\Phi }}_k\right)$
by the asymptotic expected value of Hessian

[0212] $\frac{1}{2}\left(d_{\mathrm{kx}}+d_{\mathrm{ky}}\right)\log \left(m\right),$
we can obtain the BIC criterion for B-spline model selection as follows:

[0213] $M=\arg \underset{M_k,k=1,\dots k_{\max }}{\max }\left\{\begin{array}{c}-\frac{m}{2}\log {\hat{\sigma }}_{\mathrm{kx}}^2\left({\hat{\Phi }}_{\mathrm{kx}}\right)-\frac{m}{2}\log {\hat{\sigma }}_{\mathrm{ky}}^2\left({\hat{\Phi }}_{\mathrm{ky}}\right)\\ -\frac{1}{2}\left(d_{\mathrm{kx}}+d_{\mathrm{ky}}\right)\log \left(m\right)\end{array}\right\}$
where d_kx and d_ky are the number of control points in the x and y directions respectively, and m is the number of data points.

[0214] In the above equation, the first two terms {circumflex over (σ)}_kx² and {circumflex over (σ)}_ky²

measure the prediction accuracy of the B-spline model, which increases with the complexity of the model.

[0215] In contrast, the second term decreases and acts as a penalty for using additional parameters to model the data. However, since the {circumflex over (σ)}_kx² and {circumflex over (σ)}_ky² only depend on the training sample for model estimation, they are insensitive when under fitting or over fitting occurs. In the above equation, only penalty terms prevent the occurrence of over-fitting. In fact, an honest estimate of σ_kx² and σ_ky² should be based on a re-sampling procedure. Here, the available data may be divided into a training sample and a prediction sample. The training sample is used only for model estimation, whereas the prediction sample is used only for estimating the prediction data noise {circumflex over (σ)}_kx² and {circumflex over (σ)}_ky². That is, the training sample is used to estimate the model parameter {circumflex over (Φ)}_k by: Φ_x=[B^TB]⁻¹B^Tx, Φ_y=[B^TB]⁻¹B^T_y while the prediction sample is used to predict data noise σ_k² by

[0216] ${\hat{\sigma }}_{\mathrm{kx}}^2\left({\hat{\Phi }}_{\mathrm{kx}}\right)=\frac{1}{m}\sum _{k=0}^{m-1}{\left[x_k-B_k{\hat{\Phi }}_{\mathrm{kx}}\right]}^2.$
In fact, if the model {circumflex over (Φ)}_k fitted to the training data is valid, then the estimated variance {circumflex over (ρ)}_k from a prediction sample should also be a valid estimate of data noise.

[0217] In another preferred embodiment, for a given a set of models insert p51 and data r, the Bayesian approach selects the model with the largest posterior probability. The posterior probability of model M_k may be denoted by:

[0218] $p\left(M_k❘r\right)=\frac{p\left(r❘M_k\right)p\left(M_k\right)}{\sum _{L=1}^{k_{\max }}p\left(r❘M_L\right)p\left(M_L\right)}$

• [0219] where p(r|M_k) is the integrated likelihood of model M_k and p(M_k) is the prior probability of model M_k. To find the model with the largest posterior probability, evaluate
  p(M_k|r) for k=1,2, . . . ,k_max
  • [0220] and select the model that has the maximum p(M_k|r), that is

[0221] $\begin{array}{c}M=\arg \underset{M_k,k=1,\dots k_{\max }}{\max }\left\{p\left(M_k❘r\right)\right\}\\ =\arg \underset{M_k,k=1,\dots k_{\max }}{\max }\left\{\frac{p\left(r❘M_k\right)p\left(M_k\right)}{\sum _{L=1}^{k_{\max }}p\left(r❘M_L\right)p\left(M_L\right)}\right\}\end{array}$

[0222] Here, we assume that the models have the same likelihood a priori, so that p(Mk)=1/k_max, (k=1 . . . k_max). Therefore, the model selection in

[0223] $p\left(M_k❘r\right)=\frac{p\left(r❘M_k\right)p\left(M_k\right)}{\sum _{L=1}^{k_{\max }}p\left(r❘M_L\right)p\left(M_L\right)}$
will not be affected by p(M_k). This is also the case with

[0224] $\sum _{L=1}^{k_{\max }}p\left(r❘M_L\right)p\left(M_L\right)$
since it is not a function of M_k. Consequently, the factors p(M_k) and

[0225] $\sum _{L=1}^{k_{\max }}p\left(r❘M_L\right)p\left(M_L\right)$
may be ignored in computing the model criteria.

[0226] Equation

[0227] $\begin{array}{c}M=\arg \underset{M_k,k=1,\dots k_{\max }}{\max }\left\{p\left(M_k❘r\right)\right\}\\ =\arg \underset{M_k,k=1,\dots k_{\max }}{\max }\left\{\frac{p\left(r❘M_k\right)p\left(M_k\right)}{\sum _{L=1}^{k_{\max }}p\left(r❘M_L\right)p\left(M_L\right)}\right\}\end{array}$
then becomes

[0228] $M=\arg \underset{M_k,k=1,\dots k_{\max }}{\max }\left\{p\left(r❘M_k\right)\right\}$

[0229] To calculate the posterior probability of model M_k, we need to evaluate the marginal density of data for each model p(r|M_k) which requires multidimensional integration:

[0230] $p\left(r❘M_k\right)={\int }_{{\Phi }_k}^{}p\left(r❘{\Phi }_k,M_k\right)p\left({\Phi }_k❘M_k\right)d{\Phi }_k$
where Φ_k is the parameter vector for model M_k, p(r|φ_k,M_k) is the likelihood and p(Φ_k|M_k) is the prior distribution for model M_k.

[0231] In practice, calculating the multidimensional integration is very difficult, especially for obtaining a closed form analytical solution. The research in this area has resulted in many approximation methods for achieving this. The Laplace's approximation method for the integration appears to be a simple one and has become a standard method for calculating the integration of multi-variable Gaussians [53]. This gives:

[0232] $\begin{array}{c}p\left(r❘M_k\right)={\int }_{{\Phi }_k}p\left(r❘{\Phi }_k,M_k\right)p\left({\Phi }_k❘M_k\right)d{\Phi }_k\\ \cong {\left(2\pi \right)}^{d_k/2}{H\left({\hat{\Phi }}_k\right)}^{-1/2}p\left(r❘{\hat{\Phi }}_k,M_k\right)p\left({\hat{\Phi }}_k❘M_k\right)\end{array}$
where {circumflex over (Φ)}_k is the maximum likelihood estimate of Φ_k, d_k denotes the number of parameters (control points for B-spline model) in model M_K, and H({circumflex over (Φ)}52 is the Hessian matrix of −log p(r|Φ_k,M_k) evaluated at {circumflex over (Φ)}_k,

[0233] $H\left({\hat{\Phi }}_k\right)=\frac{{\partial }^2\log p\left(r❘{\Phi }_k,M_k\right)}{\partial {\Phi }_k\partial {\Phi }_k^T}{❘}_{{\Phi }_k={\hat{\Phi }}_k}$

[0234] This approximation is particularly good when the likelihood function is highly peaked around {circumflex over (Φ)}_k. This is usually the case when the number of data samples is large. Neglecting the terms of p({circumflex over (Φ)}_k|M_k) and using log in the calculation, the posterior probability of model M_k becomes:

[0235] $M=\arg \underset{M_k,k=1,\dots k_{\max }}{\max }\left\{\log p\left(r❘{\hat{\Phi }}_k,M_k\right)-\frac{1}{2}\log H\left({\hat{\Phi }}_k\right)\right\}$

[0236] The likelihood function p(r|{circumflex over (Φ)}_k,M_k) of a closed B-spline cross section curve may be factored into x and y components as
p(r|{circumflex over (Φ)}_k,M_k)=p(x|{circumflex over (Φ)}_kz,M_k)·p(y|{circumflex over (Φ)}_ky,M_k)
where {circumflex over (Φ)}_kx and {circumflex over (Φ)}_ky may be calculated by

[0237] $\{\begin{array}{c}{\Phi }_k={\left[B^{T}B\right]}^{-1}{B}^{T}x\\ {\Phi }_y={\left[B^{T}B\right]}^{-1}{B}^{T}y\end{array}$

[0238] Consider the x component. Assuming that the residual error sequence is zero mean and white Gaussian with a variance σ_kx²({circumflex over (φ)}_kx). The likelihood function may be denoted as follows:

[0239] $p\left(x❘{\hat{\Phi }}_{\mathrm{kx}},M_k\right)={\left(\frac{1}{2\pi {\sigma }_{\mathrm{kx}}^2\left({\hat{\Phi }}_{\mathrm{kx}}\right)}\right)}^{m/2}$ $\exp \left\{-\frac{1}{2{\sigma }_{\mathrm{kx}}^2\left({\hat{\Phi }}_{\mathrm{kx}}\right)}\sum _{k=0}^{m-1}{\left[x_k-B_k{\hat{\Phi }}_{\mathrm{kx}}\right]}^2\right\}$ $\mathrm{with}{\sigma }_{\mathrm{kx}}^2\left({\hat{\Phi }}_{\mathrm{kx}},M_k\right)\mathrm{estimated}\mathrm{by}$ ${\sigma }_{\mathrm{kx}}^2\left({\hat{\Phi }}_{\mathrm{kx}}\right)=\frac{1}{m}\sum _{k=0}^{m-1}{\left[x_k-B_k{\hat{\Phi }}_{\mathrm{kx}}\right]}^2$

[0240] Similarly, the likelihood function of the y component may also be obtained. The corresponding Hessian matrix Ĥ_k of −log p(r|Φ_k,M_k)

[0241] Evaluated at {circumflex over (Φ)}_k is

[0242] $H\left({\hat{\Phi }}_k\right)=\left[\begin{array}{cc}\frac{B^{T}B}{{\hat{\sigma }}_{\mathrm{kx}}^2\left({\hat{\Phi }}_{\mathrm{kx}}\right)}& 0\\ 0& \frac{B^{T}B}{{\hat{\sigma }}_{\mathrm{ky}}^2\left({\hat{\Phi }}_{\mathrm{ky}}\right)}\end{array}\right]$

[0243] Approximating

[0244] $\frac{1}{2}\log H\left({\hat{\Phi }}_k\right)$
by the asymptotic expected value of Hessian insert

[0245] $\frac{1}{2}\left(d_{\mathrm{kx}}+d_{\mathrm{ky}}\right)\log \left(m\right)$
the Bayesian information criterion (BEC) for selecting the structure of a B-spline curve is

[0246] $M=\arg \underset{M_k,k=1,{\mathrm{\dots k}}_{\max }}{\max }\left\{\begin{array}{c}-\frac{m}{2}\log {\hat{\sigma }}_{\mathrm{kx}}^2\left({\hat{\Phi }}_{\mathrm{kx}}\right)-\frac{m}{2}\log {\hat{\sigma }}_{\mathrm{ky}}^2\left({\hat{\Phi }}_{\mathrm{ky}}\right)\\ -\frac{1}{2}\left(d_{\mathrm{kx}}+d_{\mathrm{ky}}\right)\log \left(m\right)\end{array}\right\}$
where d_kx and d_ky are the number of control points in x and y directions respectively, and m is the number of data points.

[0247] In the conventional BIC criterion as shown in the above equation, the first two terms measure the estimation accuracy of the B-spline model. In general, the variance {circumflex over (σ)}_k² estimated from

[0248] ${\hat{\sigma }}_{\mathrm{kx}}^2\left({\hat{\Phi }}_{\mathrm{kx}}\right)=\frac{1}{m}\sum _{k=0}^{m-1}{\left[x_k-B_k{\hat{\Phi }}_{\mathrm{kx}}\right]}^2$
tends to decrease with the increase in the number of control points. The smaller the variance value in {circumflex over (σ)}_k², the bigger the value of the first two terms (as the variance is much smaller than one) and therefore the higher the order (i.e. the more control points) of the model resulting from

[0249] $M=\arg \underset{M_k,k=1,\dots k_{\max }}{\max }\left\{\begin{array}{c}-\frac{m}{2}\log {\hat{\sigma }}_{\mathrm{kx}}^2\left({\hat{\Phi }}_{\mathrm{kx}}\right)-\frac{m}{2}\log {\hat{\sigma }}_{\mathrm{ky}}^2\left({\hat{\Phi }}_{\mathrm{ky}}\right)\\ -\frac{1}{2}\left(d_{\mathrm{kx}}+d_{\mathrm{ky}}\right)\log \left(m\right)\end{array}\right\}$

[0250] However, if too many control points are used, the b-spline model will over-fit noisy data points. An over-fitted B-spline model will have poor generalization capability. Model selection thus should achieve a proper tradeoff between the approximation accuracy and the number of control points of the B-spline model. With a conventional BIC criterion, the same data set is used for estimating both the control points of the B-spline model and the variances. Thus the first two terms in the above equation cannot detect the occurrence of over fitting in the B-spline model selected.

[0251] In theory, the third term in the above equation could penalize over-fitting as it appears directly proportional to the number of control points used. In practice, however, it may be noted that the effect of this penalty term is insignificant compared with that of the first two terms. As a result, the conventional BIC criterion is rather insensitive to the occurrence of over-fitting and tends to select more control points in the B-spline model to approximate the data point, which normally results in a model with poor generalization capability.

[0252] The reason for the occurrence of over-fitting in conventional BIC criterion lies in the way the variances σ_kx² and σ_ky² are obtained. A reliable estimate of σ_kx² and σ_ky² should be based on re-sampling of the data, in other words, the generalization capability of a B-spline model should be validated using another set of data points rather than the same data used in obtaining the model.

[0253] To achieve this, the available data may be divided into two sets: a training sample and a prediction sample. The training sample may be used only for model estimation, whereas the prediction sample may be used only for estimating data noise σ_kx² and σ_ky².

[0254] For a candidate B-spline model M_k with d_kx and d_ky control points in the x and y directions, the BIC may be evaluated via the following steps:

1) Estimate the model parameter {circumflex over (Φ)}_k using the training sample by

[0255] $\hspace{1em}\{\begin{array}{c}{\Phi }_x={\left[B^{T}B\right]}^{-1}{B}^{T}x\\ {\Phi }_y={\left[B^{T}B\right]}^{-1}{B}^{T}y\end{array}$
2) Estimate the data noise σ_k² using the prediction sample by equation

[0256] ${\hat{\sigma }}_{\mathrm{kx}}^2\left({\hat{\Phi }}_{\mathrm{kx}}\right)=\frac{1}{m}\sum _{k=0}^{m-1}{\left[x_k-B_k{\hat{\Phi }}_{\mathrm{kx}}\right]}^2$

[0257] If the model {circumflex over (Φ)}_k fitted to the training data is valid, then the estimated variance {circumflex over (σ)}_k² from the prediction sample should also be a valid estimate of the data noise. It may be seen that the data noise σ_k² estimated from the prediction sample may be more sensitive to the quality of the model than one directly estimated from a training sample, as the variance σ_k² estimated from the prediction sample may also have the capability of detecting the occurrence of over-fitting.

[0258] Thus, in one or more preferred embodiments, a Bayesian based approach may be adopted as the model selection method. Based on the posterior model probabilities, the Bayesian based approach estimates a probability distribution over an ensemble of models. The prediction is accomplished by averaging over the ensemble of models. Accordingly, the uncertainty of the models is taken into account, and complex models with more degrees of freedom are penalized. Given a set of models {M_k,k=1, 2, . . . , k_max} and data r, the Bayesian approach selects the model with the largest posterior probability. To find the model with the largest posterior probability, we evaluate p(M_k|r) for k=1, 2, . . . , k_max and select the model that has the maximum p(M_k|r), that is

[0259] $\begin{array}{c}M=\arg \underset{M_k,k=1,\dots k_{\max }}{\max }\left\{p\left(M_k❘r\right)\right\}\\ =\arg \underset{M_k,k=1,\dots k_{\max }}{\max }\left\{\frac{p\left(r❘M_k\right)p\left(M_k\right)}{\sum _{L=1}^{k_{\max }}p\left(r❘M_L\right)p\left(M_L\right)}\right\}\end{array}$