Jacques Feldmar

AND Nicholas Ayache

Fabienne Betting

3D-2D Projective Registration of Free-Form Curves and Surfaces

projection) which maps a 3D object onto a 2D image of this object, the relative positions of the 3D object and the Some medical interventions require knowing the correspondence between an MRI/CT image and the actual position of 2D sensor being unknown.

the patient. Examples occur in neurosurgery and radiotherapy,

This problem can be solved with artificial markers visible but also in video surgery (laparoscopy). We present in this in both images, 1 but this produces unacceptable constraints paper three new techniques for performing this task without (e.g., preoperative images must be taken on the same day artificial markers. To do this, we find the 3D-2D projective as the operation and stereotactic frames are very painful transformation (composition of a rigid displacement and a perand can prevent a free access by the surgeon). spective projection) which maps a 3D object onto a 2D image Recently, computer vision techniques have been proof this object. Depending on the object model (curve or surface), posed to solve this registration problem without artificial and on the 2D image acquisition system (X-Ray, video), the markers. Grimson [17] uses an intermediate laser range techniques are different but the framework is common:

finder, which provides a 3D description of the patient's • We first find an estimate of the transformation using bisurface. This surface is then matched against the surface tangent lines or bitangent planes. These are first order semiof the segmented corresponding surface in the volumetric differential invariants. medical image. As the laser range finder is calibrated with • Then, introducing the normal or tangent, we define a disrespect to the camera, the medical image can be fused tance between the 3D object and the 2D image, and we minimize with the video image. Colchester [9] and Kanade [29] also it using extensions of the Iterative Closest Point algorithm. developed reconstruction/rigid registration frameworks.

• We deal with the critical problem of outliers by computing

The approach adopted by Schweikard [28] is a 2D correla-Mahalanobis distances and performing generalized 2 tests.

tion scheme between the radiography and precomputed Results are presented on a variety of real medical data to simulated radiographies of a MRI image. Another correlademonstrate the validity of our approach.

INTRODUCTION

in the 2D image to perform the 3D-2D registration task. This paper is a contribution to this new and exciting Medical images are commonly used to help establish a research field. We introduce three different techniques for correct diagnosis. As they contain spatial information, both 3D-2D registration, depending on the particular clinical anatomical and functional, they can also be used to plan problem. We first present in Section 2 an approach which therapy and even in some cases to control the therapy. A makes use of passive stereo to reconstruct the surface. The recent overview of such research can be found in [START_REF] Ayache | Analysis of three-dimensional medical images-Results[END_REF] and problem is in this case 3D-3D surface rigid registration. in [START_REF] Taylor | Computer Accepted for publication in Int. J. Comput. Vision[END_REF], a spectacular use of planning and control of therapy Then, in Section 3 we present a framework for finding the using medical images and robots can be found in [7, 21, geometric transformation between a 3D image of vessels 31] for surgery and [28] for radiotherapy. Some basic techand a 2D radiography of those vessels. This is the 3D-2D niques involved in these problems are presented in [16].

curve registration problem. Finally, in Section 4, we present One of the most difficult tasks is to register, if possible a technique which allows us to find the transformation in real time, a video or X-ray image of the patient (intraopbetween a 3D surface and an image of the surface, by using erative image) with a preoperative image (MRI or CT).

the silhouette. This is the 3D-2D surface registration Basically, one must find the 3D-2D projective transformaproblem. tion (composition of a rigid displacement and a perspective For the three registration problems addressed in this paper, we have developed a unified framework. The subsections correspond to each step of this approach:

• First, we find an initial estimate of the geometric transformation. To do that, we use bitangent planes or bitangent lines, which are semidifferential invariants 2 [START_REF] Van Gool | Semi-differ-REFERENCES[END_REF] involving only first-order derivatives. Hence, our framework does not depend on the initial relative positions of the objects.

• Then, we find an accurate transformation by defining a distance between the 3D and the 2D object (which is a compromise between the spatial distance and the differ-FIG. 1. A surface and two bitangent points M 1 and M 2 . Let n 1 and ence in normal or tangent orientation) and by minimizing n 2 be the normals at these points. M 1 and M 2 are bitangent if the plane defined by (M 1 , n 1 ) and the plane defined by (M 2 , n 2 ) are the same. it using two-step minimization algorithms which are exten-tracted from the MRI image S 2 . So, the actual algorithm where Match is a function from the 3D space to S 2 . The iterative algorithm minimizes this energy E. In step 1, the is slightly more complex, but is basically as stated. The details can be found in Appendix B.1, and an analysis of variables (R, t) are fixed and E is minimized with respect to Match. The function which minimizes E in this case is the complexity of the search for an initial estimate can be found in Appendix B.2. We simply note that the complexity ClosestPoint. In step 2, the variable Match is fixed and E is minimized with respect to (R, t). Hence, at each step, of the algorithm is quasi-linear in the number of points on the surfaces and that the risk of stopping the algorithm E decreases. Because E is positive, convergence is guaranteed. with a wrong initial estimate decreases extremely quickly with the number of points (when the two surfaces actually This ICP algorithm is efficient and finds the correct solution when the initial estimate (R 0 , t 0 ) of the rigid displace-show some overlapping regions up to a rigid displacement). ment is ''not too bad'' and when each point on S 1 has a 2.2. The Iterative Closest Point Algorithm correspondent on S 2 . But in practice, this is often not the case. For example, in our application, as explained in the Using the pairs of bitangent points as described in the previous section, the reconstructed surface usually only previous subsection, we get an estimate (R 0 , t 0 ) of the rigid partially describes the patient's surface and often includes displacement to superpose S 1 on S 2 . In order to find an a description of the patient's environment. The next two accurate rigid displacement we have developed an extensubsections explain how we deal with these two problems. sion of an algorithm called ''the Iterative Closest Point algorithm'' which was introduced by several researchers [4-6, 26, 36]. 4 We sketch the original ICP algorithm, which 2.3. Working with Incomplete Surfaces searches for the rigid displacement (R, t) which minimizes

In step 1 of the iterative algorithm, we map each point the energy of S 1 to a ''closest point'' on S 2 . But when the two surfaces are partially reconstructed, some points on S 1 do not have

E(R, t) ϭ M i ʦS i ʈRM i ϩ t Ϫ ClosestPoint(RM i ϩ t)ʈ 2 ,
any homologous point on S 2 . Thus, given a point M on S 1 ,( R i Ϫ 1 ,t i Ϫ 1 ), and ClosestPoint(R iϪ1 M ϩ t iϪ1 ), we have where ClosestPoint is the function which associates to a to decide whether (M, ClosestPoint(R iϪ1 M ϩ t iϪ1 )) is a space point its closest point on S 2 .

plausible match. This is very important because, if we ac-The algorithm consists of two iterated steps, each iteracept incorrect matches, the found rigid displacement will be tion i computing a new estimation (R i , t i ) of the rigid disbiased (and therefore inaccurate), and if we reject correct placement.

matches, the algorithm may not converge toward the best solution. 1. The first step builds a set Match i of pairs of points.

As proposed in [1], we make use of the extended Kalman The construction is the following: for each point M on S 1 , filter (EKF). This allows us to associate to the six paramea pair (M, N) is added to Match i , where N is the closest ters of (R i , t i ) a covariance matrix S i and to compute a point on S 2 to the point R iϪ1 M ϩ t iϪ1 . To compute the generalized Mahalanobis distance ͳ for each pair of closest point, several methods have been proposed, for matched points (M, N). This generalized Mahalanobis disexample the distance map method [START_REF] Danielsson | Euclidean distance mapping[END_REF].

tance, under some assumptions on the noise distributions 2. The second step is the least-squares evaluation of the and some first-order approximations, is a random variable rigid displacement (R i , t i ) to superpose the pairs of Match i with a 2 probability distribution. By consulting a table of (see for example [12] who use a quaternion-based method). values of the 2 distribution, it is easy to determine a confidence level for ͳ corresponding to, for example, a The termination criterion depends on the approach used: 95% probability of having the distance ͳ less than .I n the algorithm stops either when (a) the distance between this case, we can consider the match (M, N)a slikely or the two surfaces is below a fixed threshold, (b) the variation plausible when the inequality ͳ Ͻis verified and consider of the distance between the two surfaces at two successive any others as unlikely or unplausible. iterations is below a fixed threshold, or (c) a maximum This distinction between plausible and unplausible number of iterations is reached.

matches implies a change in the second step of the iterative The convergence of the algorithm is clearly demonalgorithm. Given Match i , instead of computing the rigid strated in [START_REF] Cohen | Use of auxiliary variables in computer vision problems[END_REF]. Let us define the energy displacement (R, t) which minimizes the least-squares criterion

E(R, t, Match) ϭ M i ʦS 1 ʈRM i ϩ t Ϫ Match(M i )ʈ 2 , (M,N)ʦMatch i ʈR i M ϩ t i Ϫ Nʈ 2 ,
we recursively estimate the six parameters of (R, t), and the where n 1 (M) is the normal on S 1 at point M and ClosestPoint 6D is the new 6D closest point function. associated covariance matrix which minimizes the criterion In step one, the closest point now has to be computed in 6D space. Because we cannot extend the technique (M,N)ʦMatch i and (M,N)isplausible (R i M described in [START_REF] Danielsson | Euclidean distance mapping[END_REF] (the distance map image would be much too large), we use the kd-tree technique first proposed by Zhang ([36]) for the 3D case. The second step must also be modified: the criterion which defines the best rigid displacement must use the new 6D distance. Otherwise, it is where W j is a covariance matrix associated with each pair not possible to prove the convergence of our new ICP (M, N) which allows us, for example, to increase the imporalgorithm as in Subsection 2.2. Hence, the rigid displacetance of high curvature points. ment (R i , t i ) is now defined as the minimum of the function More details about the meaning of ''plausible or not'' and about the EKF can be found in Appendix C. In practice, we use extended Kalman filters to minimize minimum. To reduce this problem, we propose in this subthis new criterion at step 2. Even though it is in nonlinear, section to make use of the normal information and to the minimization works very well. We have carried out define a new criterion to minimize. In our formulation, numerous experiments on synthetic data, and the global surface points are no longer 3D points: they become 6D minimum was always found even with rather large data points. Coordinates of a point M on the surface S are (x, noise and rather crude initial estimates. Note that this use y, z, n x , n y , n z ), where (n x , n y , n z ) is the normal to S at M. of extended Kalman filters allows us to compute Mahala-For two points M(x, y, z, n x , n y , n z ) and N(xЈ, yЈ, zЈ, nЈ x , nobis distances and to determine if a match is plausible or nЈ y , nЈ z ) we define the distance not as explained in the previous subsection.

ϩ t i Ϫ N) t W Ϫ1 j (R i M ϩ t i Ϫ N),
f (R, t) ϭ (M,N)ʦMatch i d(RM ϩ t, N)
In order to try to demonstrate that the ICP algorithm, using the 6D distance, converges more often to the global d(M, N) ϭ (Ͱ 1 (x Ϫ xЈ) 2 ϩ Ͱ 2 (y Ϫ yЈ) 2 minimum than the standard ICP algorithm, we conducted ϩ Ͱ 3 (z Ϫ zЈ) 2 ϩ Ͱ 4 (n x Ϫ nЈ x ) 2 the following experiment. We chose S 1 and S 2 to be the same surface. Hence, the resulting transformation should ϩ Ͱ 5 (n y Ϫ nЈ y ) 2 ϩ Ͱ 6 (n z Ϫ nЈ z ) 2 ) 1/2 , be the identity. We run both the original and the modified algorithm choosing different initial rigid displacements (R 0 , t 0 ), at an increasing distance from the identity. The where Ͱ i is the inverse of the difference between the maxiresults are reported in Table 1. This shows that our modimal and minimal value of the ith coordinate of points in fied algorithm is in practice much less sensitive to the initial S 2 . Using this definition of distance, the closest point to P estimate (R 0 , t 0 ), and more robust to local minima. on S 2 is a compromise between the 3D distance and the difference in normal orientation. 52.5. Results This new definition of the distance between points naturally implies modifications to steps one and two of the ICP We now present an example application of the framealgorithm in order to minimize the new energy work presented in this section. First we compute on both the stereo reconstructed surface, and on the MRI surface (Fig. 2), the pairs of bitangent

E(R, t) ϭ MʦS 1 d((RM ϩ t, Rn 1 (M)),
points. We find 598 pairs on the stereo surface and 5000 pairs on the MRI surface (obviously, the desired density ClosestPoint 6D((RM ϩ t, Rn 1 (M))) 2 , of the bitangent pairs is a parameter of the bitangent extraction algorithm). Hence, there are more pairs on the MRI surface than on the stereo surface which makes the Note. We chosen for S 1 and S 2 the same surface: 100 random points chosen on a mannequin head surface. Then we randomly chose rigid displacements, ran both the original and the modified algorithm, and verified if the result was the identity or not. The laws for choosing the rigid displacement are Gaussian laws centered on the identity. For each covariance matrix, we made 1000 trials and we report in the table the number of successes. For example, for a standard deviation of 0.1 radians for the rotation and 500 mm for the translation, the standard ICP algorithm found the correct solution 640 times out of 1000 trials, whereas the result was 906 using the normal information. Another example is for a standard deviation of 1.0 radians for the rotation and 10 mm for the translation: the number of successes are respectively 151 and 360. For more information on Gaussian laws for choosing rotations, see [27]. We note that the standard deviation for rotations is a number between 0 and ȏ.

Using these pairs of bitangent points, we estimate the have their closest point at a distance lower than 3 mm.

The average distance between matched points is 1.6 mm. rigid displacement in about 30 s. Applying this estimate, 80% of the points on the stereo surface have their closest The result is presented in Fig. 3, left.

Because we know the point-to-point correspondences point at a distance lower than 8 mm. This error must be compared with the size of the voxel in the MRI image: between the MRI head surface and the stereo face surface (this is the result of the registration), and because of each 4mm ϫ 4m mϫ2 mm. Moreover, recall that there are points on the stereo surface which do not a homologous point of the stereo surface we know the gray level from the video image, we can map the video image onto the point on the MRI surface.

Using this estimate of the rigid displacement, we run MRI surface (Fig. 3, right). The fact that the points on the MRI surface have the right gray levels qualitatively the modified iterative closest point algorithm. The MRI head surface is described by 15000 points and the stereo demonstrates that the MRI/stereo matching is correct.

Once the MRI head surface has been colorized, we can surface by 10000 points. It takes 20 s. Applying this new rigid displacement, 85% of the points on the stereo surface apply the projective transformation computed as explained FIG. 2. The bitangent lines computed on the MRI (left) and stereo (right) face surfaces. Note that we selected the lines whose length varies between 2 cm and 10 cm. in this paper and get a gray-level 2D image. This image is (for example in the brain or kidney) and a 2D angiography of those vessels (Fig. 8). Two techniques are used to get very similar to the video image and this demonstrates that the 3D images of the vessels. It can either be a MRI image we actually recover the right projective transformation where all the vessels are visible, or it can be a CT image. (Fig. 4).

In the latter case, opaque liquid is injected in a vessel and Finally, we projected the brain onto the video image only those vessels linked to the one injected are visible in (see Fig. 5) using the computed projective transformation.

the 3D image. The 2D image is a standard radiography In fact, we now have enough geometric and textural paramand some opaque liquid is also injected in order to make eters to produce a stereo pair of realistic images from a the vessels visible. continuous range of viewpoints and provide the surgeon An interesting application of the 3D-2D vessel registrathe feeling of seeing inside the patient's head and guide tion is in interventional radiography. A common operation him/her during the operation as explained in [17].

for the radiologist is to introduce a catheter into an artery and guide it toward a lesion within a vessel. The radiologist

3D-2D PROJECTIVE CURVE REGISTRATION

injects opaque fluid and takes radiographies throughout the operation (which can be quite long) in order to see In this section, we consider the problem of finding the projective transformation between a 3D image of vessels the vessels and the catheter. The opaque liquid can be

FIG. 5.

The projection of the brain in the two video images using the projective transformation computed as explained in this paper. A stereoscopic display could provide the surgeon with the feeling of seeing inside the patient's head. The two presented images cannot be visually fused in this position, because the baseline between the two optical centers of the video cameras was vertical in this experiment: one can notice that the camera of left image was bellow the camera of the right one. This will be corrected in further experiments. dangerous for the patient, and the X-radiations are noxious FocalPlane having equation z ϭ 1. This means that a 3D point M ϭ (X, Y, Z) is transformed to the point for both the radiologist and the patient. Also, understanding the 3D shape of the vessels from two radiographies is m ϭ (X/Z, Y/Z, 1) which is the intersection between the line a very hard task.

OM and the plane FocalPlane (see Fig. 6). For clarity, Knowledge of the projective transformation which maps we make this assumption in the sequel of this article. the 3D vessels onto the 2D ones enables us to determine

We consider that a point m ϭ (x, y, 1) belonging to the position of the catheter (which is visible in the radiograthe FocalPlane is a 2D point with coordinates (x, y). phy) with respect to the 3D image. Hence, the radiologist We call the perspective transformation which transforms could visualize his catheter in the 3D image during the the 3D point M to the 2D point m, the transformation operation. This is a significant step toward finding the right Proj, and we denote this by m ϭ Proj(M). The problem path within the vessels. Since this reduces the risk of an now is to find the rigid displacement (R, t) such that, incorrect interpretation, the operation becomes safer.

for all points M i on the 3D curve corresponding to the Moreover, by reducing the interpretation time, the quanpoint m i on the 2D curve, m i ϭ Proj(Rm i ϩ t). The tity of injected liquid and radiations can be lowered.

projective transformation search for is Proj o (R, t). Some recent work ( [START_REF] Malandain | Topological segmentation of discrete surfaces[END_REF]) enables us to represent the vessels as curves described by points in both the 3D and the 2D images. Moreover, spline approximation yields the tangents to these curves [18]. Hence, from a computer vision point of view, the problem is to find the best projective transformation which maps a given 3D curve onto a given 2D curve, the two curves being described by points and their associated tangents. This problem has already been addressed for special curves such as circles or ellipses but not for free-form curves to our knowledge (a hint can be found in [19] for nonplanar parameterized curves).

Note that we know the camera calibration matrix T corresponding to the 2D image processing. 7 So, it is possible to assume without loss of generality (see Appendix E) that the camera image formation corresponds to perspective projection with center O ϭ (0, 0, 0) and a focal plane FIG. 6. Each point M i on the 3D curve is first transformed into RM i ϩ t and then transformed into Proj (RM i ϩ t). Finally, an affine 7 We make the pinhole camera assumption in this paper. This is well justified for our applications.

transformation of FocalPlane provides the image coordinates.

C which lie in the plane defined by (O, m 1 , m 2 ). This proves that T 1 and T 2 are coplanar.

SECOND BITANGENT LINE PROPERTY. Let t be the common tangent line on c at m 1 and m 2 . Let m 3 be an intersection point between t and c. m 3 is the projection of a point M 3 which is an intersection point between the plane defined by (M 1 , M 2 , T 1 , T 2 ) and the curveC .

Proof. O, m 1 , m 2 , m 3 , M 1 , M 2 , T 1 , T 2 all lie in the same plane P. Because M 3 is on the line Om 3 , M 3 also lies in P.

We use these two properties to find an initial estimate of the projective transformation which maps a given 3D curve C onto a given 2D curve c as follows: 3. If no solution was found, return to point 1.

More details on this algorithm can be found in Appendix D.1. We note that the complexity is quasi-linear in the 3.1. Finding an Initial Estimate of the number of points describing the curve. It allows us to find Projective Transformation an estimate of the projective transformation which maps A well-known result is that knowledge of three points the 3D curve onto the 2D curve. But, because we need M 1 , M 2 , M 3 on the 3D object and of the three corresponda more accurate transformation, we use this estimate to ing 2D points m 1 , m 2 , m 3 in the 2D image yields four rigid initialize a two-step minimization algorithm. displacements (R, t) which satisfy m i ϭ Proj(RM i ϩ t) 3.2. Extension of the ICP Algorithm to 3D-2D (see for example [START_REF] Linnainmaa | Pose determination of a three-dimensional object using triangle pairs[END_REF]). Note that there are six equations Curve Registration and six unknowns. Of course, it is not possible to select all the 3D and 2D triplets and to explore the possible 3.2.1. First Extension. The transformation of a 3D point combinations in order to find an estimate of the transfor-M by the estimate of the projective transformation found mation. We must select from the 3D and 2D curves two by the method explained in the previous subsection corresubsets of triplets which can match each other. In order sponds to: (i) the application of a rigid displacement (R 0 , to reduce the number of potential matches, we use two t 0 )toMand, (ii) the application of the projection Proj. In projective properties of curves (see Fig. 7):

order to get a more accurate projective transformation, we propose to minimize the energy FIRST BITANGENT LINES PROPERTY. Let c be a 2D curve which is the perspective projection of a 3D curveC( cϭ

E ( Ͱ ,ͱ ,Ͳ ,t x ,t y ,t z )ϭ M i ʦ C
ʈ Proj(RM i ϩ t) Proj(C)). Two points on c which share the same tangent lines are projections of two points on C such that the tangent lines are coplanar.

Ϫ ClosetPoint(Proj(RM i ϩ t))ʈ 2 , Proof. Let M 1 and M 2 be two points on C. Let T 1 and T 2 where (R, t) is the rigid displacement corresponding to the be the attached tangents. Let m 1 and m 2 be the projection of M 1 and M 2 .Ifm 1 and m 2 share the same tangent line, these tangent lines are projections of two 3D tangent lines on three Euler angles Ͱ, ͱ, Ͳ, t is the translaction vector (t x , t y , where T i is the normalized tangent vector to C at point M i and Proj tangent is the function which associates with t z ), and where ClosestPoint is the function which associates with each point m in FocalPlane its closest point on c. the normalized tangent vector RT i at point RM i ϩ t, the normalized tangent vector to the projected curve at point As in the previous section, the algorithm consists of two iterated steps, each iteration j computing a new estimation Proj(RM i ϩ t). The modifications to the previous extended ICP algorithm to minimize this new criterion are straight-(R j , t j ) of the rigid displacement: forward and similar to the modifications in the 3D-3D 1. Build a set of pair of points Match j by attaching to rigid case. In step 1, the closest points are computed in 4D each point M i in C the closest point m i on c to the point space using the kd-tree technique, and in step 2 the 4D Proj(R jϪ1 M i ϩ t jϪ1 ).

criterion corresponding to the new 4D distance is 2. Compute the best rigid displacement which minimizes

(M i ,m i )ʦMatch j d(Proj(RM i ϩ t), m i ) 2 . M i ʦC ʈProj(RM i ϩ t) Ϫ m i ʈ 2 .
As in the 3D-3D case, experiments demonstrate that the Using the same kind of demonstration as the converuse of tangent information makes the minimization algogence demonstration of the ICP algorithm presented in rithm less sensitive to its initialization. Subsection 2.2, it is straightforward to prove that the algorithm converges and minimizes the energy E. Moreover, 3.3. Results because we also use extended Kalman filters to perform the minimization at step 2 of the algorithm, we get the 3.3.1. 3D Morphometer Data. We present an applicacovariance matrix associated to the rigid displacement tion of this framework to the data shown in Fig. 8. The (R jϪ1 , t jϪ1 ). Hence, as in the 3D-3D rigid case, we can 3D data have been acquired by the ''3D morphometer''. compute for each pair in Match j a Mahalanobis distance Figure 9 (left) shows the initial estimate. One can oband decide whether a pair in Match j is a plausible match serve that it is quite different from the one shown in Fig.

or not. This allows us to deal efficiently with the outliers 8 (left). Indeed, different projections of the same 3D curve problem. Note that this problem is for example crucial for can produce very different images. This demonstrates that the 3D-2D brain vessel registration problem because all the search for the initial transformation (Subsection 3.1) the vessels are visible in the MRI image, whereas only the is crucial. Figure 9 (right) shows the result of the minimizainjected ones are visible in the radiography.

tion process presented in Subsection 3.2 from this initial estimate. The computation time is 20 s on a DEC alpha 3.2.2. Extension Using the Tangent Information. As for workstation. For matched points, the average distance is the 3D-3D rigid registration problem, the exploitation of 6.8 pixels for the initial estimate and 0.83 pixel for the tangent information enhances the convergence of the final registration. 3D-2D ICP algorithm. In this case, the points on the 2D It is difficult to estimate the accuracy of the registration curves are no longer 2D points. They are 4D points (x, y, with this measure of the error. An improved quantitative tan x, tan y), where (x, y) are the classical spatial coordimeasure would require two radiographies taken with a 90Њ nates and (tan x, tan y) are the coordinates of the normalangle. This would allow us to measure the z error. We ized tangent vector to the curve. Hence, given two points plan to conduct such experiments in the near future. m 1 and m 2 with coordinates (x 1 , y 1 , tan x 1 , tan y 1 ) and (x 2 , y 2 , tan x 2 , tan y 2 ), we can define the distance 3.3.2. 3D MR Aniography Data. We also conducted some experiments to match 3D MR angiographies onto

d(m 1 , m 2 ) ϭ (Ͱ 1 (x 1 Ϫ x 2 ) 2 ϩ Ͱ 2 (y 1 Ϫ y 2 ) 2
2D X-ray angiographies. This scheme has a large potential use in practice because the acquisition of MR angiography

ϩ Ͱ 3 (tan x 1 Ϫ tan x 2 ) 2 ϩ Ͱ 4 (tan y 1 Ϫ tan y 2 ) 2 ) 1/2 ,
is noninvasive. This is the reason why we treat the vessels as 3D curves. Indeed, the vessels generally appear thinner where Ͱ i is the inverse of the difference between the maxiin MR images than in 2D X-rays. Even if the skeleton is mal and minimal value of the ith coordinate of points in c. not invariant under projective transformation, this approx-In this case, the criterion corresponding to this new disimation works well in practice because the vessels are tance is thin enough.

A stereotactic frame was attached to the head of the

E(Ͱ, ͱ, Ͳ, t x , t y , t z ) ϭ ͚ M i ʦC d(RM i ϩ t), Proj tangent(RT i ))
, patient during the acquisition of the images. This allowed us to compare quantitatively the projective transformation ClosestPoint 4D((Proj(RM i ϩ t), Proj tangent(RT i )))) 2 , T frame found using the stereotactic frame and the transfor- mation T frameless found without the frame as described in frame technique is probably more robust. We believe that using two 2D angiographies, a frontal one and a lateral this section, using only the vessels visible in the images.

In order to compare T frame and T frameless , we define an one, the frameless technique could be robust enough for clinical use. Note that when several radiographies are avail-error measure as follows. For each point M in the 3D volume of interest, we apply T frame to M. We get a 2D able, and assuming that the relative acquisition positions are known, it is then straightforward to define a new crite-point m. The inverse of m for the transformation T frameless is a set of points MЈ forming a line L. We define the error rion being the sum of the individual criterions define above for each image and therefore to take into account both associated with the point M as the distance between M and the line L.

radiographies. Typically, we get an average error of 0.82 mm, a maximal error of 2.04 mm and a minimal error of .0001 mm. This

3D-2D PROJECTIVE SURFACE REGISTRATION

error is close to the uncertainty attached to the estimation of T frame . Hence, from an accuracy point of view the Even if a vessel is not visible both in the 3D image and in the 2D radiography, it can be useful to find the projective ''frame'' and ''frameless'' techniques are similar. The enormous advantage of the frameless technique is that it is transformation which maps the 3D image onto the 2D radiography. This is for example the case in radiotherapy noninvasive and permits comparison between images taken with a long time interval. On the other hand, the (as explained in [28]), where one problem is to link the

MRI image to the actual position of the patient in order

We call Normal 2D the function which associates the normal n to the point m and we note n ϭ Normal 2D (m). to guide the radiation sources. We recall that the usual practice to complete this task is to screw a stereotactic 4.2. Finding an Initial Estimate of the frame in the skull of the patient! Projective Transformation The approach described in this section resembles the approach introduced by Lavallee in [START_REF] Lavalle ´e | Matching 3-d smooth surfaces 1. N. Ayache. Artificial Vision for Mobile Robots-Stereo-Vision and with their 2-d projections using 3-d distance maps[END_REF]. We find the posi-As in the previous section, we use two properties of tion of the 3D object so that its projection corresponds surface projection in order to select independently on the to the occluding contours detected in the 2D image. But 3D surface, and on the occluding contour, triplets of points Lavalle ´e's approach assumes that it is possible to determine which can correspond. By exploring possible associations globally the interior and the exterior of the occluding conof 2D and 3D triplets and verifying them, we will find an tour. This is often a hard task in practice. The use of the estimate of the required projective transformation. These vector normals allows us to avoid this problem. Moreover, properties are: we do not make any assumption about the initial position PROPERTY 1.9 Two 2D points m 1 and m 2 on the occludof the 3D object. Ponce also used the occluding contour ing contour which share the same tangent line are projections in [20] to find the pose of 3D objects, but his scheme of two points M 1 and M 2 on the surface which share the requires second-order derivatives, which can significantly same tangent plane. reduce the robustness of the approach.

Proof. All tangent lines to the surface at points M 1 and 4.1. Fundamental Property of Occluding Contour M 2 lie in the plane (O, m 1 , m 2 ). The basic idea of our approach is to use a PROPERTY 2. 10 Let m 1 and m 2 be two 2D points on the FUNDAMENTAL PROPERTY OF OCCLUDING CONTOUR (see occluding contour and let M 1 and M 2 be two points on the Fig. 10). If a point M on a 3D surface S is so that projection 3D surface such that m 1 ϭ Proj(M 1 ) and m 2 ϭ Proj(M 2 ). point m ϭ Proj(M) lies on the occluding contour c, then The angle between Normal 2D (m 1 ) and Normal 2D (m 2 ) is the normal vector Normal 3D (M) to S at point M is equal to the normal vector n of the plane P defined by (m, O, t), equal to the angle between the normals Normal 3D (M 1 ) and x-and y-coordinates of points on c. The energy corresponding to this new distance is Normal 3D (M 2 ) on the surface at points M 1 and M 2 .

Proof. The two planes defined respectively by (m i , E(Ͱ, ͱ, Ͳ, t x , t y , t z , Match) ϭ

m i ʦc d(m i , RMatch(m i ) ϩ t) 2 , Normal 2D (m i )) and (M i , Normal 3D (M i )) are identical.
Hence, to find the initial estimate (R 0 , t 0 ) of the rigid where Match is a function which associates a point on S displacement so that Proj o (R 0 , t 0 ) is the projective transwith each point on c. This energy measures how well the formation which transforms S with the 2D occluding con-''fundamental property of occluding contour'' (Subsection tour c, we proceed as follows:

4.1) is verified. 1. Choose (m 1 , m 2 , m 3 ) on the 2D curve c such that:

As in the two previous sections, we use a two-step algo-(a) m 1 and m 2 share the same tangent line; and (b) the rithm to minimize this energy. This algorithm can be underangle Ͱ between Normal 2D (m 1 ) and Normal 2D (m 3 )i sa s stood as an extension of the ICP algorithm. The first step close as possible to ȏ.

is the minimization of E with respect to Match, assuming that Ͱ, ͱ, Ͳ, t x , t y , t z are fixed. The second step is the 2. For each triplet (M 1 , M 2 , M 3 ) on the 3D surface S minimization of E with respect to Ͱ, ͱ, Ͳ, t x , t y , t z assuming such that (a) M 1 and M 2 share the same tangent plane; and that Match is fixed. We now describe the algorithm which (b) the angle between Normal 3D (M 1 ) and Normal 3D (M 3 )is consists of two iterated steps, each iteration j computing Ͱ, compute the projective transformations which maps a new rigid displacement (R j , t j ): (M 1 , M 2 , M 3 ) onto (m 1 , m 2 , m 3 ). Stop if a computed projective transformation is a correct initial estimate. 111. We associate with each point m i on c a point M i .To 3. If no solution was found, return to point 1.

do so, we first compute the set More detail about this algorithm (in particular about the CloseNormal(m i ) way in which we deal with discretization error) can be found in Appendix D.2. We note that the complexity is ϭ ͕M ʦ S/Angle(R t jϪ1 Normal 2D (m i ), Normal 3D (M)) Ͻ͖, quasi-linear in the number of points describing the surface. As in the previous section, this algorithm yields an estimate where is a given threshold. 12 Next we compute the point of the searched projective transformation and we use it N which minimizes d(m i , R jϪ1 M ϩ t jϪ1 ) for all M in to initialize an energy minimization algorithm and get an CloseNormal. In order to ensure that E decreases at this accurate transformation.

step 1 (this is the condition required to prove the convergence of the algorithm), we compare d(m i , R jϪ1 N ϩ t jϪ1 ) 4.

Extension of the ICP Algorithm to 3D-2D

and d(m i , R jϪ1 M i Prev ϩ t jϪ1 ), where M i Prev is the Surface Registration point associated with m i in the previous iteration. 13 If d(m i , R jϪ1 N ϩ t jϪ1 ) is smaller than d(m i , R jϪ1 M i Prev ϩ The previous estimate of the projective transformation t jϪ1 ), the matched point M i becomes N, otherwise it remains corresponds to, first, the application of a rigid displacement as M i Prev. (R 0 , t 0 ) and, then, the application of the projection Proj.

2. Compute the rigid displacement (R j , t j ) which mini-In order to get a more accurate projective transformation, mizes we propose a distance between the points on the surface and the points on the occluding contour and minimize an 2 . appropriate energy function.

m i ʦc d(m i , RM i ϩ t)
Let m ϭ (x, y) be a point on c and let M ϭ (X, Y, Z) be a point on S. We define the 5D distance

The demonstration of the convergence of this algorithm is straightforward, since at each step E decreases. More- 2 over, because we use extended Kalman filters to perform ϩ (Normal 2D (m) Ϫ Normal 3D (M)) 2 ) 1/2 , the minimization of step 2, we get the covariance matrix corresponding to the six parameters of (R j , t j ). Hence, as in the two previous sections, we can compute a Mahalanobis where Ͱ 1 and Ͱ 2 are the inverse of the difference between distance and deal with points on c which have no correthe maximum and the minimum values of respectively the spondent on S by performing 2 tests (see Appendix C for more details). This allows us to find an accurate projective cause they are considered as outliers (they will be considered later during the iterative process), this shows that transformation which transforms the surface S with occludthe initial estimate is already quite good. For the final ing contour c, discarding outliers.

d(m, M) ϭ (Ͱ 1 (x Ϫ X/Z) 2 ϩ Ͱ 2 (y Ϫ Y/Z)
transformation the average error is 0.76 pixels and 0.17 4. [START_REF] Besl | A method for registration of 3-D shapes[END_REF]

. Results

degrees which demonstrates that it is accurate. We plan to attach a stereotactic frame to the mannequin head and We now present an application of the framework preperform the experiment again in order to have a reference sented in this section for surface registration using occludmeasure of the error and then try to better demonstrate ing contours.

the accuracy of the recovered transformation. Note that Figure 11 presents the data to register: a 3D model of when the relative positions of the cameras are known, it a mannequin head extracted from a 3D CT-scan image is possible to integrate into our framework several camera (left), and a 2D video image of this mannequin head (right). views. The criterion is, in this case, the sum for each image Figure 12 (left) shows the initial estimate, and Fig. 12 of the criterion defined in the previous subsection. (right) shows the result of the minimization stage. The CPU time is 10 s on a DEC alpha workstation.

THE NONCALIBRATED CASE

For the initial estimate, the average error is 2.1 pixels for spatial distance and 2 degrees for the difference in

In the two previous sections, we assumed that we knew the calibration matrix attached to the camera. Even if this normal orientation. Even if some points are rejected be- seems to be a reasonable assumption because in practice surface and the X-ray image. This yields the estimate of the projective transformation (Fig. 14,left). The average it is possible to calibrate a camera (see for example [1]), or because it is quite often possible to know the calibration error is 1.5 pixels for spatial distance and 1.9 degrees for the difference in normal orientation. Then we run the parameters for a radiology system, we now explain how to proceed when we do not know the calibration matrix. 14 It minimization algorithm to recover accurately the 11 parameters of the projective transformation. The average is obvious that the problem is harder. Since five calibration parameters are unknown, the problem is to find all the error for the resulting transformation is 0.79 pixels and 0.7 degrees which demonstrates that it is accurate. eleven parameters of the projective transformation.

CONCLUSION

Finding an Initial Estimate of the Projective Transformation

We presented a unified framework to perform 3D-2D We do not (yet) have a satisfactory automatic solution projective registration. We first showed how to use the for this problem. In practice, we find the estimate manually. bitangent lines or the bitangent planes to find an estimate Using an appropriate interface, it is possible to select interof the transformation. Then we proposed definitions of the actively pairs of corresponding points on the 3D object and distance between the 3D object and the 2D image and we the 2D object. Selecting at least six pairs of corresponding presented extensions of the iterative closest point algopoints, it is possible, with a least-squares minimization, to rithm to minimize this distance dealing with occlusion in determine an initial estimate.

an rigorous way.

Finding an Accurate Projective Transformation

Validation

In practice, there are two cases where we are provided Although we presented results on real data, a lot of work with an estimate of the best projective transformation, and still needs to be done. Because the techniques described in we want to find the best one with respect to the eleven this article are related to medical operations, it is especially parameters. This is the case when the estimate has been important to perform a rigorous validation. Two properties found manually as explained in the previous subsection. must be established precisely: robustness and accuracy. In But this problem also arises when the camera calibration the presented framework, we first segment the data, then is not very accurate, a common occurrence in practice.

find an initial estimate of the projective transformation In such cases, we minimize the same criteria as the ones and finally we perform an energy minimization. described in the two previous sections. But instead of min-

The segmentation stage is obviously very important. imizing them with respect to six parameters, we minimize Both robustness and accuracy depend on it. We did not them with respect to the eleven parameters of the projecdetail this stage in this paper, more details can be found tive transformation. The minimization algorithms are alin [START_REF] Malandain | Filtrage, topologie et mise en correspondance d'images[END_REF][START_REF] Malandain | Topological segmentation of discrete surfaces[END_REF]. In practice, it works well and the registration most the same. We observe that these 11D minimizations remains robust thanks to our treatment of the outliers. are more sensitive to the initialization conditions than the Our search for an initial estimate only involves first-6D case. But, because our initial estimates are not too order derivatives and has reasonable algorithmic complexcrude, it works well in practice. This allows us to find ity. We could have chosen to compute higher order derivaan accurate projective transformation even if the camera tives to achieve a lower computational complexity or to calibration is not perfect, or when the initial estimate has compute lower order derivatives to achieve a better robeen found interactively.

bustness. We made the choice presented in this article because it seems to us, after experiments, that it is a good 5.

Results

trade-off between robustness and complexity (see [START_REF] Van Gool | Semi-differ-REFERENCES[END_REF] for an interesting discussion on these issues). In practice, the We have a 3D CT-scan image of a skull and a 2D X-Ray different thresholds have been easily chosen and the image of it (Fig. 13, right). We do not know the calibration method is not sensitive to their precise values. But we will parameters of the radiography system. We first extract the not discuss this point here. In practice, an appropriate skull surface in the 3D image (Fig. 13, left). Then, we user interface allows the user to design interactively some design six pairs of corresponding points on the 3D skull corresponding points on the objects to register. Hence, even if an initial estimate is proposed automatically as 14 We do note here about the 3D-3D registration problem addressed described in this paper, it can be corrected interactively if in Section 2. But note that it is possible to perform the stereo with it is wrong. Note that the initial estimate does not need noncalibrated cameras. In this case, the reconstructed surface is known to be accurate because it is just an initialization of the up to a 3D-3D projective transformation. Then the corresponding regisminimization stage. tration problem would be to find the 3D-3D projective transformation (15 parameters) to superpose the two surfaces.

The robustness of the minimization stage is not easy to check. It is possible to perform statistics. One can add apy, it would be useful to combine the techniques described in Sections 2 and 4. Indeed, the first one is well adapted some noise to the objects to register, perturb the correct solution, run the algorithm, and check if it converges to-for frontal views and the second one for lateral views. Note that it is not difficult to combine different viewpoints with ward the correct solution. Such experiments were presented in Table 1. The next stage is to perform the experi-the different techniques described in this paper since the defined criterions can easily be added together, as noted ments on a larger data set and it will take a long time. Checking the accuracy is difficult but important too. The in Subsection 3.3.2. best way to do it is probably to compare the automatic 6.2. Future Work registration with the one obtained with a stereotactic frame. Such an experiment was presented in Subsection

We believe that the techniques described in this paper offer different possible and interesting extensions.

3.3.2.

Even if the experiments presented in this paper seem First, we would like to compare from a practical point of view the techniques described in this paper with respect to demonstrate the robustness and accuracy of the techniques (at least in the calibrated case), we believe that for to the techniques presented by others [9,17,[START_REF] Lavalle ´e | Matching 3-d smooth surfaces 1. N. Ayache. Artificial Vision for Mobile Robots-Stereo-Vision and with their 2-d projections using 3-d distance maps[END_REF]22,28,29,35]. For example, we plan to develop new procedures to maximal security, it is better to use more than one sensor when possible. For neuradiology, it would be better to use validate and compare rigorously the accuracy of these techniques. two X-ray views as noted in Subsection 3. are not randomly distributed (they belong to a smooth we choose W k corresponding to the covariance of the estimated error after the correct registration, and S 0 as the surface). But this can help to understand why, in practice, for our problems, the choice of the three parameters ͳ, ⌬, diagonal matrix corresponding to the covariance of the error on the initial estimation of the rigid displacement pa-and is not difficult and why a good solution is found after a very small number of iterations (two or three).

rameters.

The problem of outliers is to distinguish between likely matches and unlikely matches in Match j . For each pair

FIG. 3 .

 3 FIG. 3. (left) The result of the rigid registration. The MRI surface is clearer and the stereo surface is darker. The alternation dark/clear shows that the registration is quite accurate. (right) The colorized MRI head surface obtained by attaching to each matched point the gray level of its corresponding point on the stereo surface. Thanks to the transparency effect, one can observe the brain.

FIG. 4 .

 4 FIG. 4.The left image shows one of the video images. The right image is the gray-level image obtained by applying to the colorized MRI surface the recovered projective transformation. Let us call this image the ''simulated image.'' Note that the white areas come from unmatched points. In order to demonstrate that the video and the simulated images are very similar, we computed the contours C in the video image: C is visible in the left picture. Then we superimposed C int the simulated image. One can observe that C fits very well with the simulated image. This demonstrates that the recovered projective transformation is correct.

1. 7 .

 7 Choose a triplet (m 1 , m 2 , m 3 ) of points on the 2D curve c such that (a) the tangent lines at point m 1 and m 2 are the same line l; and (b) m 2 is an intersection point between c and l. 2. For each triplet (M 1 , M 2 , M 3 ) of points on the 3D curve C such that (a) the tangent lines at point M 1 and M 2 lie in the same plane P; and (b) M 3 is an intersection point FIG. Illustration of the two properties of curves projection used between C and P, 8 compute the projective transformations in this paper. Two points m and n which share the same tangent lines l which maps (M 1 , M 2 , M 3 )o n( m 1 ,m 2 ,m 3 ). Stop if a are projections of two 3D points M and N such that the tangent lines at computed projective transformation is a correct initial esthese points lie in the same plane Plane. Moreover, if p is an intersection point between l and the 2D curve, then it is the projection of a 3D point timate. P which is an intersection point between the 3D curve and Plane.

FIG. 8 .

 8 FIG. 8. (left) The 3D model of the vessels. It was extracted from a 3D image acquired by the ''3D morphometer'' of General Electric Medical System Europe. (right) A radiography of the vessels (an angiogram). Our goal is to find the point to point correspondence between the left and right objects.

FIG. 9 .

 9 FIG. 9. (left)The initial estimate of the transformation. (right) The result of the minimization process. One can observe that the radiography vessels are no longer visible which shows that the 3D model is quasi-perfectly projected on them. This demonstrates the accuracy of the found projective transformation.

FIG. 10 .

 10 FIG. 10. The illustration of the fundamental property of occluding controur formation. The projection of the point M is m. m lies on the occluding contour and the tangent to this contour is t. The normal vector to S at point M is equal to the normal to the plane defined by (O, m, t).

FIG. 11 .

 11 FIG. 11. (left) The 3D model of the mannequin head. It was extracted from a 3D CT-scan image. Thanks are due to Cannes Hospital. (right) A video picture of the mannequin head. Our goal is to find the point-to-point correspondence between the left and the right objects.

FIG. 12 .

 12 FIG.12. (left) The initial estimate of the transformation: we projected the 3D model using the initial transformation. The contour superimposed in this image is the contour extracted in the video image. (right) The final transformation: we projected the 3D model using the transformation found after minimization. You can observe that the contour which comes from the video image exactly corresponds to the contours of the 3D model. This demonstrates that we recovered the correct projective transformation.

FIG. 13 .

 13 FIG. 13. (left) The skull surface extracted from the 3D CT-scan image. (right) The skull X-ray image. The goal is to recover the 3D-2D projective transformation (11 parameters) mapping the 3D surface (left) onto the 2D image (right), without any calibration information.

  FIG. 14. (left) The transformation found designing six pairs of corresponding points on the 3D surface and the 2D image. (right) The result of the minimization. You can observe that the contour which comes from the X-ray image exactly corresponds to the contours of the 3D model (except small areas where the contour extraction is not perfect and which are considered as outliers). This demonstrates that we recovered the right projective transformation.

TABLE 1 Superiority of the Algorithm Using the Surface Normal Information: Quantitative Experiments

 1

	r (radians)-t (mm)	0.1-0	0.25-0	0.5-0	0.75-0	1.0-0	1.5-0
	points	1000	963	598	293	157	61
	points ϩ normals	1000	999	881	596	370	170
	r (radians)-t (mm)	0.1

-10 0.25-10 0.5-10 0.75-10 1.0-10 1.5-10 points

		1000	964	582	286	154	68
	points ϩ normals	1000	997	858	567	360	168
	r (radian)-t (mm)	0.1-100	0.25-100	0.5-100	0.75-100	1.0-100	1.5-100
	points	997	921	512	247	139	52
	points ϩ normals	1000	995	834	549	348	152
	r (radian)-t (mm)	0.1-500	0.25-500	0.5-500	0.75-500	1.0-500	1.5-500
	points	640	576	312	158	85	33
	points ϩ normals	906	870	695	456	295	147

The problem is then very similar to the camera calibration problem. * E-mail: Jacques.Feldmar@sophia.inria.fr.

Note that an extension of this algorithm to the nonrigid case can be found in [13] and[START_REF] Szeliski | Matching 3-d anatomical[END_REF].

Of course, only two parameters are necessary to describe the orientaalgorithm more efficient. The extraction requires about tion of the normal (for example the two Euler angles). But we use (n x , n y , n z ) because the Euclidean distance better reflects the difference of 30

s.[START_REF] Chen | Object modeling by registration of multiple range images[END_REF] orientation between the normals (that is not the case with the Euler angles because of the modulo problem) and we can use kd-trees to find the closest point as explained later.[START_REF] Chen | Object modeling by registration of multiple range images[END_REF] CPU times are given for a DEC-ALPHA workstation.

The algorithms to compute the set of triplets are given in Appendix A.

Note that this property has already been used by Zisserman[START_REF] Van Gool | Semi-differ-REFERENCES[END_REF] to where t is the tangent vector to the occluding contour at find the pose of symmetric objects.

Note that property 1 is a consequence of property 2. point m.

A technique for verifying the hypothesis is proposed in[START_REF] Basri | The alignment of objects with smooth surfaces[END_REF], using the surface curvature, and making the affine camera assumption. We

Note that the computation of this set is fast because we compute, as a preprocessing, a 3D kd-tree on the normals to the surface S. note that we use only the surface normals and that our camera model is the pinhole one (see Appendix D.2).

The convention is that d(m i , R 0 M i Prev ϩ t 0 ) is infinity.

A rigid displacement (R, t) maps each point M to RM ϩ t where R i sa3ϫ3 rotation matrix and t a translation vector.to set a value for p and q, and the points on the surfaces

This is the same kind of assumption as in Appendix A.1.Working with the new set of 2D points, because A is an

Most often in our applications, the contours are closed and the affine transformation of FocalPlane, everything happens distance from the object to the camera is large with respect to the size of the object, therefore it is possible to choose such a point.as if the camera image formation process was a perspective

ACKNOWLEDGMENTS

ential invariants, in Geometric Invariance in Computer Vision (J. L. Mundy and A. Zisserman, Eds.), MIT Press, Cambridge, MA, We thank the members of the Epidaure group for their practical help 1992. and stimulating discussions. We also thank the members of the Robotvis 16. W. Grimson, Object Recognition by Computer: The Role of Geometric group at INRIA and especially Frederic Devernay and Cyril Zeller. Some Constraints, MIT Press, Cambridge, MA, 1990. very interesting discussions with Mike Brady (currently invited professor 17. W. E. L. Grimson, T. Lozano-Perez, W. M. Wells III, G. J. Ettinger, in our group at INRIA) and Gerard Medioni also influenced our work. S. J. White, and R. Kikinis, An automatic registration method for Thanks are also due to General Electric Medical System Europe frameless stereotaxy, image guided surgery, and enhanced reality (GEMS) (Catherine Picard and Yves Trousset) for providing us with visualization, in IEEE Proceedings of Computer Vision and Pattern some of the data presented here. Thanks to Dr. Michel Royon (Cannes Recognition 1994 (CVPR '94), Seattle, June 1994. Hospital) and Dr. Jean-Jacques Baudet (La Timone Hospital) for provid-18. A. Gue ´ziec and N. Ayache, Smoothing and matching of 3-D-space ing us data and for fruitful discussions. Finally, thanks to Drs. David curves, Int. J. Comput. Vision 12(1), 1994. Dean and Court Cutting for skull data. This work was supported in part by a grant from Digital Equipment 19. R. M. Haralick and Yu Hong Chu, Solving camera parameters from Corporation and by the European Basic Research Action VIVA (Esprit the perspective projection of a parameterized curve, Pattern Recognit. project). Future work will be also partially supported by a grant from 17(6), 1984, 637-645. GEMS. 20. T. Joshi and J. Ponce, Hot curves for modelling and recognition of smooth curved 3d objects, in IEEE Proceedings of Computer Vision and Pattern Recognition 1994 (CVPR '94), Seattle, June 1994.

Then, we wish to validate our approach on a larger scale. angle must depend on the local curvature of the surface and on the local density of the discretization. In practice, At least for the vessel registration problem, we plan a clinical validation in collaboration with General Electric we choose M as the average angle between the normal at point M and the normal at the four closest points in S to Medical Systems Europe and a hospital. We will use a frontal and a lateral view for these experiments.

M. In order to control the final density of the set Bitangent(S), we introduce a parameter ͳ: we add in Bi-We believe that it should be possible to perform realtime tracking of the patient and enable the surgeon to move tangent(S) a new bitangent pair (M, MЈ) only if there is no pair (Q, QЈ)i nBitangent(S) so that (d(M, Q) Ͻ ͳ) either the patient or the 2D sensor. Indeed, for tracking we need only to correct a projective transformation which is and d(MЈ, QЈ) Ͻ ͳ)o r( d ( M ,Q Ј )Ͻͳ )and d(MЈ, Q) Ͻ ͳ). We use a hash table for performing this test. quite close to the right solution as explained in Subsection 5.2. Hence, because the initialization would be good, it

The complexity of the technique is O(p ϫ n), where p is the maximum number of points in CloseNormal(M) and should be possible to use the extensions of the iterative closest point algorithm presented in this paper with just a n is the number of points describing the surface. In practice, because of the way we compute M , the more the discretiza-few points without local minimum problem and to get fast convergence for the eleven parameters of the projective tion of S is dense, the more M is small. Hence, in practice, p depends only on the ''complexity'' of the surface 15 and transformation.

Finally, we would like to develop a new 3D-2D registra-the technique is quasi-linear. Of course, we do not pretend that this analysis is rigorous: we just try to explain why the tion scheme, for example extending some ideas presented by Fua in [START_REF] Fua | Registration without correspondences, in projection of center O ϭ (0, 0, 0[END_REF] and try to better understand the 3D-2D algorithm is efficient in practice. registration problem when the 3D object is deformable.

A.2. Computing the Coplanar Tangents on a 3D Curve This problem arises, for example, in the processing of X-ray mammographies.

The problem is here to compute, given a 3D curve described by points and attached tangents, the pairs of points

APPENDIX A

such that the tangent lines are coplanar. The technique is very simple. A set of m planes is associated with each Computing Special Triplets of Point tangent line. It is just a uniform sampling of the set of planes containing the tangent line.

A.1. Computing the Bitangent Planes on a

Hence, the problem is the same as in Appendix A.1: Surface and the Bitangent Line on a 2D Curve computing the bitangent planes on a surface. Because m is a constant number (typically 20), the complexity is also We present in this appendix our technique for computquasi-linear. ing, given a surface S described by points and attached normals, the pairs of points sharing the same tangent plane.

A.3. Computing the Intersections between a Set of 2D

We do not present the details of the technique for comput-Lines and a 2D Curve and between a Set of ing the bitangent lines on a 2D curve because it is the same Planes and a 3D Curve problem in 2D space.

As a preprocessing, we first compute a kd-tree Tree The problem is, given a set of 2D lines and a 2D curve based on the Euclidean coordinates of the normal vectors. described by points, to compute the intersection points Then for each point M, we compute the set between the lines and the curve. We first tessellate the plane using squares of edge size d. Hence, each point of the curve belong to a square. For each nonempty square CloseNormal(M) ϭ ͕N ʦ S/Angel(Normal 3D(M), Square, we compute the set CloseLine of lines which are Normal 3D(N)) Ͻ M ͖, at a distance below d ϫ ͙2 ϩ ͳ from the square center.

Then, exploring all the curve points in Square and all the lines in CloseLine, we simply retain the set of pairs (M, l) where Normal 3D(M) is the normal to S at point M, angle such that M is a point in Square and l is a line in CloseLine is the function which returns the angle between two norand the distance between M and l is smaller then ͳ.W e malized vectors, and M is a threshold automatically comchoose ͳ as the average distance between two consecutive puted as explained later. For each point MЈ in CloseNorpoints on the curve. mal(M), let P (resp. PЈ) be the tangent plane at point M The complexity of the technique is in the worse case (resp. MЈ) and let m (resp. mЈ) be the orthogonal projection O(p ϫ n), where p is the number of lines and n the number of M on P (resp. PЈ). We decide that M and MЈ are bitangent points if the (M ٙ MЈ, MmЈ) and (MЈ ٙ M, MЈm) are both smaller than M . 15 Clearly, we assume that the surface is not flat but smooth. This is a reasonable assumption for anatomical surfaces.

The key parameter of the algorithm is M . This threshold of points describing the curve. So, in the worst case the The parameter ͳ is the estimation of the maximal distance between one point on the surface S 1 and its closest plane tessellation is unuseful. But in practice, it makes the computation much faster and the calculation time is ac-point on S 2 for a good estimate of the best rigid registration (ͳ depend on the noise level). is the estimation of the ceptable.

We do not present the details for computing the intersec-ratio of the number of points on S 1 which have a correspondent on S 2 , divided by the number of points on S 1 ( tions between a set of planes and a 3D curve. Indeed, it is the same problem as the previous one in 3D space.

depends on the occlusion level).

B.2. Complexity Analysis

APPENDIX B

The goal of this subsection is not to provide a very Finding the Initial Estimate of the 3D-3D rigorous analysis of the algorithm complexity but just to Rigid Registration try to explain ''why it works.'' First, if the surface does not present too many flat areas, B.1. The Algorithm the complexity of the bitangent calculation and the number In practice, two corresponding pairs of bitangent points of bitangents is quasi-linear in the number of points on cannot be exactly superposed because of the point discretithe surface (see Appendix A.1). zation error and because of the error in the normal compu-Moreover, the number of iterations of the algorithm is tation. Moreover, only a part of the reconstructed surface relatively small because we just have to try a few pairs P 1 S 1 may be superposed on the patient's surface extracted in the set Bitangent(S 1 ) (step 1 of the algorithm). Indeed, from the MRI image S 2 . Indeed, in S 1 the patient's surface assume that the probability for a pair P 1 in Bitangent(S 1 ) is not complete, and there can be a part with corresponds to have no corresponding pair in Bitangent(S 2 )isq. After to the reconstruction of the patient's environment. So, in n iterations the probability that we didn't verify a right practice, we first compute on S 1 and S 2 the two sets of hypothesis is q n . Then, if the verification stage is correct, pairs of bitangent points Bitangent(S 1 ) and Bitangent(S 2 ). the number of iterations will be very small. Then, we sort the set Bitangent(S 2 ) on the distance between

We now focus on the complexity of each iteration. It the two bitangent points. This way, given a pair P 1 in is simply ͉CloseDistance(P 1 )͉ ϫ Verification Complexity, Bitangent(S 1 ), we can very quickly compute the set where Verification Complexity is the complexity of the verification stage. The number of pairs in CloseDistance(P 1 ) CloseDistance(P 1 ) depends on ⌬. But even if ⌬ can be relatively small in practice, ͉CloseDistance(P 1 )͉ is porportional to ͉Bitangenϭ ͕P 2 ʦ Bitangent(S 2 )/͉d(P 1 ) Ϫ d(P 2 )͉ Ͻ⌬͖, t(S 2 )͉, that is quasi-proportional to the number of points where ⌬ is a given threshold. on S 2 . Finally, Verification Complexity is constant. Indeed, Then we apply the following algorithm:

we can a priori fix the number nЈ of points in SЈ 1 which are used to verify. Moreover, we can choose nЈ relatively small: 1. We randomly choose a pair P 1 in Bitangent(S 1 ) the risk to accept a bad hypothesis is still very small. In-2. We compute the set CloseDistance(P 1 ) deed, assume that (R, t) is a rigid displacement which 3. For each pair P 2 in CloseDistance(P 1 ), we compute corresponds to a wrong hypothesis. If the points are ranthe two rigid displacements D and DЈ which correspond domly distributed within a volume and if for a point M ʦ to the superposition of P 1 on P 2 . For each of these two RS 1 ϩ t the probability that there exists a rigid displacements, we now estimate the number of points point N ʦ S 2 such that ʈMN UUǞ ʈ Ͻ ͳ is p, then the probability on the transformed surface RS 1 ϩ t 16 which have their to step the algorithm with this bad rigid displacement is closest point on S 2 below a given distance ͳ to check whether D or DЈ reaches our termination criterion. First, iϭn iϭk nϪi , we randomly choose a subset SЈ 1 of points on S 1 . Then, for each point P in SЈ 1 , we compute the closest point Q to RP ϩ t on S 2 . Let m be the number of points P such that where n ϭ ͉SЈ 1 ͉ and k ϭ ϫ n. For instance, if p ϭ 0.4, ʈPQʈ Ͻ ͳ. If the ratio m/͉SЈ 1 ͉ is larger than (0 Ͻ Ͻ 1), ϭ 0.7 and n ϭ 500 this probability is approximatively we decide that (R, t) is a good estimate of the rigid displace-10 Ϫ41 . ments which superposes S 1 on S 2 and we stop the algorithm.

We tried to explain why the complexity of each iteration 4. If no solution has been found, we return to point 1. is quasi-linear in the number of points on the surfaces and why only a few iterations are necessary in practice. We do not pretend that this analysis is rigorous. It is not possible

APPENDIX C

(M k , M k ), we consider that we have an estimate a ˆkϪ1 of a

Using the Generalized Mahalanobis Distance and 2 with associated covariance matrix S kϪ1 and a noisy mea-Tests to Deal with Outliers surement (x ˆk , ⌳ k )o fx i , and we would like to then test the validity of the Eq. ( 1). If we consider the truncated Assume that Match j is perfect and that the computed expansion of f k (x k , a) geometric transformation perfectly superposes the two sets of points. For each pair (M k , N k )i nMatch j , we can write the equation

and assuming that w k and v k are zero-mean independent where x k is the vector obtained by concatenation of the Gaussian noise, it is clear that f k (x ˆk , a ˆkϪ1 ) is also-to within M k and N k coordinates and a is the vector of the parameters a first-order approximation-a Gaussian process whose of the geometric transformation. 17 mean and covariance are given by In reality, we cannot write Eq. ( 1) for all pairs (M k , N k ) in Match j . But we can assume that this is because each A(f k (x ˆk , a ˆkϪ1 )) ϭ 0 pair corresponds to a measurement x ˆk of x k corrupted with a random error v k :

As a result, if the rank of Q i is q, the generalized Mahalaϭ ⌳ k Ն 0, A(v i v t j ) ϭ 0 ᭙i ϶ j. nobis distance

Given the measurements x ˆk , the EKF allows us to recurͳ(x ˆk , a ˆkϪ1 )

(5) sively compute an estimate a ˆk of a and the associated covariance matrix S k (a ˆk ϭ a ϩ w k , where w k is a random error: is a random variable with a 2 probability distribution with A(w k ) ϭ 0 and A(w k w t k ) ϭ S k Ն 0). Only an initial estimate q degrees of freedom. 18 (a ˆ0 , S 0 ) is required and the minimized criterion is closely By consulting a table of values of the 2 distribution, it is related (up to a first order approximation) to easy to determine a confidence level for the Mahalanobis distance corresponding to, for example a 95% probability C ϭ (a ˆ0 Ϫ a) t S Ϫ1 0 (a ˆ0 Ϫ a) of having the distance (5) less than . In this case, we can ϩ

consider as likely or plausible the measurements x ˆk for which the inequality ͳ(x ˆk , a ˆkϪ1 ) Ͻis verified, and consider any others as unlikely or unplausible. where

Taking S Ϫ1 0 Ȃ 0, C is the criterion we want to minimize. See Ayache [1] for more de-

APPENDIX D tails.

Hence, in practice, the difficulty is just to choose the Finding the Initial Estimate of the 3D-2D Curve covariance matrix S 0 associated with the initial estimate Registration and of the 3D-2D Curve Registration (R 0 , t 0 ) of the rigid displacement and to choose the covariance matrix W k associated with each point M k . In practice,

D.1. The Algorithm for 3D-2D Curve Registration

Let c be the 2D curve and C be the 3D curve. We first compute on c the set Bitangent(c) of pair of points on c 17 For example, for the 3D-3D rigid registration problem (Section 2.2), if M k ϭ (x 1 , y 1 , z 1 ) and N k ϭ (x 2 , y 2 , z 2 ) and if (Ͱ, ͱ, Ͳ, t x , t y , t z ) are the three Euler angles and translation parameters of the searched rigid displacement, we can have x k ϭ (x 1 , y 1 , z 1 , x 2 , y 2 , z 2 ) and a ϭ (Ͱ, ͱ, Ͳ, 18 If q Ͻ p ϭ the size of the measurement vector

which share the same tangent lines. For each pair pair 2D 1. We choose a pair (m 1 , m 2 )i nBitangent 2D (c) and a point m 3 on c. in Bitangent(c) which corresponds to a tangent line t,w e compute the set of points Inter 2D (pair 2D ) which are inter-2. For each pair (M 1 , M 2 )i nBitangent 3D (S), we consections points between t and c. Then, we compute the struct the set set CoplanarTangent(C) of pairs of points on C sharing coplanar tangent lines. For each Pair 3D in CoplanarTan-Triplet 3D (M 1 , M 2 ) ϭ ͕(M 1 , M 2 , M 3 )/͉ gent(C) with attached plane P, we compute the set of points Angle(Normal 2D (m 1 ), Normal 2D (m 3 ))-Inter 3D (Pair 3D ) which are intersection points between the plane P and the curve C. The algorithms for computing Angle(Normal 3D (M 1 ), Normal 3D (M 3 ))͉ Ͻ͖, all these sets are given in Appendix A.

Let us define the two sets where Angle is the function which returns the angle between two normalized vectors and is a given threshold.

We use the kd-tree Tree to compute efficiently this set. and i ʦ Inter 2D (pair 2D )͖ 3. For each triplet (M 1 , M 2 , M 3 )i nTriplet 3D (M 1 , M 2 ), we compute the set Hypotheses of the four rigid displaceand ment (R, t) such that m i ϭ Proj(RM i ϩ t).

4. For each rigid displacement (R, t)inHypotheses,we Triplet 3D ϭ ͕(Pair 3D , I)/Pair 3D ʦ CoplanarTangent(C) and check that the transformed surface verify the ''fundamental I ʦ Inter 3D (Pair 3D )͖. property of occluding contour'' (Subsection 4.1). To do this, we compute the ratio MatchedPoint of the number We can now describe the algorithm: of points on c which have a correspondent on RS ϩ t,o n the number of points on c. To determine if a point m has 1. We choose a triplet (m 1 , m 2 , m 3 )i nTriplet 2D a correspondent on S, we proceed as follows. We first 2. For each triplet (M 1 , M 2 , M 3 )i nTriplet 3D , we comcompute the set pute the set Hypotheses of the four projective transformations which map the 3D and the 2D triplets.

CloseAngle(m) ϭ ͕M ʦ S/Angle(R t Normal 2D (m), 3. For each projective transformation in Hypotheses,we Normal 3D (M)) Ͻ͖. first verify that the projection of the tangent line attached to M 3 is the tangent line attached to m 3 . If the test is Again, to compute efficiently the set CloseAngle, we use positive, we apply the same verification scheme as dethe kd-tree Tree. If there exists a point M in CloseAngle(m) scribed in Appendix B.1 for the rigid 3D-3D case: we such that the 2D distance between m and Proj(RM ϩ t) compute the ratio MatchedPoints of the number of points is bellow a given threshold ͳ, we consider that m has a M on C such that the distance between the projection of correspondent. Otherwise, we consider that m has no cor-M and its closest point on c is bellow a given distance ͳ, respondent for the rigid displacement (R, t). If the ratio on the number of points on C.I fMatchedPoint is larger MatchedPoint is larger than a given threshold (0 Ͻ Ͻ than (0 Ͻ Ͻ1) then we stop the algorithm and the result 1), we stop the algorithm and the initial estimate is the is the current projective transformation.

projective transformation Proj o (R, t). 4. If no solution was found, we return to point 1.

5. If no solution was found, we return to point 1.

D.3. Complexity Analysis of the Initial Estimate Search D.2. The Algorithm for 3D-2D Surface Registration

for the 3D-2D Curve Registration Problem and for the To find the initial estimate (R 0 , t 0 ) of the rigid displace-3D-2D Surface Registration Problem ment such that Proj o (R 0 , t 0 ) is the projective transformation which transforms a surface S with the 2D occluding Basically, the complexity analysis for these two problems is not different from the analysis presented in Appen-contour c, we first compute the set Bitangent 2D (c) of pair of points which share the same tangent line on c and the dix B.2. Indeed the algorithms are very similar: instead of two set of pairs of points, we have two sets of triplets of set Bitangent 3D (S) of pair of points which share the same tangent plane on S. We also compute a kd-tree Tree based points. The discussion for explaining why the number of global iterations is small, why the risk to stop the algorithm on the Euclidean coordinates of the normal vectors on the surface S which will allow us, given a normalized vector with a bad hypothesis is small, and why the verification stage can be performed in constant time would be very m, to find quickly the set of normals n on S such that the angle between m and n is bellow a given threshold. Then similar to the previous one. In fact, the complexity of the algorithms depends on the number of hypotheses to verify we apply the following algorithm: at each iteration. This is the number N of 3D triplets which

The result of the calibration process is a 3 ϫ 4 matrix T. This is the matrix of the 3D-2D projective transformation can correspond to a given 2D triplet.

For the 3D-2D curve registration, we have which maps a reference object of known shape (the calibration grid) onto an image of this object. This matrix T can be decomposed into T ϭ AD. D is the 3 ϫ 4 matrix of N ϭ ͉Triplet 3D ͉ Ͻ ͉CoplanarTangent(C)͉ the extrinsic parameters and it is a matrix of change of ϫ max

͉Inter 3D (Pair 3D )͉. frame from the calibration grid frame to the camera frame. A i sa3ϫ3 matrix such that

We explained in Appendix A.2 that ͉CoplanarTangent(C)͉ is quasi-linear in the number of points in C.I fC is not linear but smooth,

It is called the matrix of the intrinsic parameters and deis small. 19 Finally, the complexity of the algorithm is scribe the camera image formation process. Basically, it is O( p ϫ n), where p is a small integer and n is the number an affine transformation of FocalPlane. of points describing C.

Let us recall why this decomposition is possible (the For the 3D-2D surface registration, we have proof is constructive). We note T ϭ [H͉t]( His a 3 ϫ 3 matrix and t a3ϫ1 matrix). For any matrix H, HH T

͉Angle (M 1 )͉, clearly a symmetric matrix. If, moreover, H is a square nonsingular matrix, then HH T is also positive definite. Indeed, for each vector v, different from 0, where v T HH T v ϭ ʈH T vʈ 2 Ͼ 0, Angle (M 1 ) ϭ ͕M 3 ʦ S/͉ Angle(Normal 3D (M 1 ), Normal 3D (M 3 ))since the kernel of H T reduced to 0 because of the nonsingularity of H. Consequently, a nonsingular upper triangular Angle(Normal 2D (m 1 ), Normal 2D (m 3 ))͉ Յ͖ matrix U exists such that ( is a given threshold). HH T ϭ UU T We explained in Appendix A.1 that ͉Bitangent 3D (S)͉ is quasi-linear in the number of points in S. p ϭ (it is one form of the Cholesky decomposition). If we max (M 1 ,M 2 )ʦBitangent 3D (S) ͉Angle (M 1 )͉ is more difficult to estidenote by O the matrix U Ϫ1 H, we see that mate. In practice, assuming that the normals on the surface S are approximately uniformly distributed on the unit OO T ϭ U Ϫ1 HH T U Ϫ1T ϭ I sphere, we try to choose m 3 such that ͱ ϭ ȏ. 20 Hence, p is small. Finally, the complexity is roughly O( p ϫ n), where n is the number of points on S. Because p is small, the which shows that O is an orthogonal matrix. Hence, we algorithm is efficient in practice. In practice, before starting the registration process, we We explain in this appendix why it is possible to assume first decompose T as explained upper. Then, we replace that the camera image formation corresponds to a perspeceach point (x, y) on the 2D curve by (xЈ, yЈ), such that tive projection of center O ϭ (0, 0, 0) on a focal plane FocalPlane of equation z ϭ 1 without loss of generality.

(xЈ, yЈ,1 ) t ϭA Ϫ 1 ( x ,y ,1 ) t .


Jacques Feldmar

Nicholas Ayache

Rigid, Affine and Locally Affine Registration of Free-Form Surfaces

published or not. The documents may come


Nicholas Ayache

Keywords:

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.


published or not. The documents may come