GB2355789A

GB2355789A - Surface registration

Info

Publication number: GB2355789A
Application number: GB9925343A
Authority: GB
Inventors: Ivan Daniel Meir; Norman Ronald Smith
Original assignee: Tricorder Technology PLC
Current assignee: Tricorder Technology PLC
Priority date: 1999-10-26
Filing date: 1999-10-26
Publication date: 2001-05-02
Also published as: AU1042001A; GB9925343D0; WO2001031581A1

Abstract

To register two surfaces by associating corresponding points on the surfaces, a point on one surface is associated with a point on the second surface by projecting the first point to a slope discontinuity (peak or ridge) on the second surface. A modified Iterative Closest Point (ICP) algorithm is used particularly for matching of pre and post-operative facial scans. The accuracy, speed of convergence and stability of the algorithm is improved by minimising the function S w<SB>i</SB> q<SB>i</SB>'-(Rq<SB>i</SB>+T) <SP>2</SP> where w<SB>i</SB> is a weighting factor which is dependent on the relationship between local surface neighbourhoods of candidate pairs of corresponding points q<SB>i</SB>' and q<SB>i</SB>. Additionally the algorithm detects points lying in regions (R1, R2) above peaks or ridges which would not normally be processed correctly and associates them with the peak or a perpendicular to the ridge as appropriate to define their closest corresponding points on the other surface, thereby avoiding instability or failure of convergence in the algorithm.

Description

2355789 Image processing method and =aratus The present invention relates

to a method and apparatus for processing 3D datasets, particularly but exclusively 3D surfacerepresentations.

A major field of application of the invention is the processing of 3D surface representations of the human body, particularly faces. It is often desirable to register 3D surfaces having at least a partial overlap eg to combine adjacent overlapping surface regions (eg as obtained by an optical scanner) into one overall surface or to compare pre- and post operative images of patients undergoing surgery.

An algorithm which is conventionally used for registration of 3D surfaces is the Iterative Closest Point (ICP) algorithm, originated by Besl and McKay (P. J. Besl, N.

D. McKay. A Method for Registration of 3-D Shapes. IEEE Trans Pattern Analysis and Machine Intelligence pp.239-256, Vol. 14, no.2, February 1992 which is incorporated herein by reference). Before describing this algorithm the general problem of registering 3D surfaces will be outlined with reference to Figures I and 2.

Figure 1 shows a "register" surface SREG which is shown in the form of atriangulated mesh of points ql, q2.... q13.... qm (in general, qi) which is to be rotated in three dimensions by a rotation matrix R and translated in three dimensions by a translation vector T so as to merge optimally with a "reference" surface SREF similarly comprising a simplex mesh of points pl, p2, p3 p13 (in general, pi) and having a region which overlaps with and corresponds to a region of the register surface.

The problem is essentially to find R and T which minimise the aggregate distances between a set of points on SREG and their corresponding closeest points on S REF. Accordingly the expression:

N qi' - (Rqi + T) 11 2 (1) is minimised in the ICP algorithm (N being the number of points in the region of overlap and qi' being the closest point on SREG to qi). Thus the sum of the squares 2 of the distances is minimised to find R and T. This is illustrated in Figure 2 which shows the operation of matrix R and vector T on surface SREG, the resulting t aggregate of distances dl 1, d22, d33.... dii being minimised. In practice the sum of 5 the squares of the distances in the region of overlap will never be zero because of measurement inaccuracies.

The efficacy of a standard ICP approach is limited by the following factors:

a) erroneous data can detrimentally influence the registration; b) the algorithm is susceptible to local minima - hence it requires initial alignment to be very close and it is hard to consistently reproduce similar results, and c) although the speed and performance of the ICP can be optimised by sampling a 15 specific subset of points from the 'register' surface, the optimal selection of points is extremely difficult to determine and to repeat. More sophisticated versions of the ICP algorithm have been devised, eg by Maurer et al (Calvin R Maurer et al, IEEE Trans Medical Imaging Vol 15 No 6, December 20 1996 pp 836 to 848, incorporated herein by reference) who apply a weighting factor approximately proportional to the size of the triangles (of the triangulated reference surface). Also the registration is based on matching of surfaces as well as points. Maurer et al examine the closest vertex on the reference surface and examine eachof its adjacent triangles until they determine the triangle to which the vertex projects. They claim that this approach produces the actual closest point on the surface 25 44approximately 99.5% of the time" and that "the occasional errors are small and do not affect the accuracy of the registration". The nature or cause of the errors is not described. Z. Zhang in "Iterative Point Matching for Registration of Free-Form Curves and 30 Surfaces" International Journal of Computer Vision Vol.1-3, No.2, pp. 119-152 (which is also incorporated herein by reference) discloses an algorithm in which the expression (1) above is modified by introducing a weighting factor which takes into account uncertainty in the positions of the points. 35 It is known that point to triangle matching is important when point to point matching is inadequate due to a low sampling density of the reference surface. We have found that an implementation of the ICP algorithm, using point to triangle matching only, which neglects the exceptional cases of ridges and peaks leads to unstable behaviour, inaccuracy, or lack of convergence.

3 An object of the present invention is to overcome or alleviate the above disadvantages of the prior art.

In one aspect the invention provides a method of registering two surfaces with corresponding surface regions by associating corresponding points of the surfaces wherein a first point on one surface facing a region of the other surface is associated with a corresponding point in that region by projecting the first point to a slope discontinuity (eg a ridge or peak) of that region.

In a preferred embodiment, if the first point lies over a ridge, an association between the first point and a perpendicular extending from the ridge to the first point is established.

In a preferred embodiment, if the first point lies over a peak, an association between the first point and the peak is established.

In another aspect the invention provides a method of iteratively registering two spatial regions with corresponding parts by repeatedly:

i) finding corresponding sets of points in the respective spatial regions, ii) relatively moving the spatial regions to minimise an aggregate of the distances between the corresponding pairs of points, and iii) finding new sets of corresponding points, wherein a weighting factor which is dependent on the relationship between local surface neighbourhoods of candidate pairs of corresponding points is applied during the iteration to accelerate the convergence or enhance the stability of the iteration.

In one embodiment the spatial regions are surfaces and the local spatial neighbourhoods are surface properties. For example the weighting factor can be inversely related to the angle between the surface normals of the candidate corresponding points, or could be inversely related to the distance between the pairs of candidate points relative to the average distance between pairs of candidate corresponding points. Also the weighting factor could be dependent upon the relationship between respective surface curvatures or textures for example.

In one embodiment the following function is minimised:

I 4 N wil I qi' - (Rqi + T) 11 2 (2) where wi is a weighting factor which is dependent on the relationship between local surface neighbourhoods of candidate pairs of corresponding points qi' and qi.

It will be noted that the weighting factor utifised by Maurer et al was a non-relative weighting factor, ie dependent only on local surface neighbourhood of a point, not 10 on the relationship between local neighbourhoods of points on both surfaces.

In a related aspect the invention relates to a method and apparatus for registering semi-rigid or partially rigid surfaces, particularly but not exclusively the human face.

It is often useful to align to surface which are not identical throughout, but share a common rigid baseline. For example, surface data of the face captured prior to maxillo-facial surgery and data capture postoperatively might share the common rigid base-line of the forehead, as opposed to other regions of the face which might have changed significantly.

Accordingly in a further aspect the invention provides a method of registering a partially rigid or semi-rigid first surface representation of a subject with a second surface representation of that subject wherein a relatively rigid portion of the first surface representation is selected and used as input to a registration algorithm so as to register with a corresponding portion of that surface representation. For example 25 the algorithm can be the modified ICP algorithm described above.

We have found that the relatively rigid portions tend to be registered preferentially by the algorithm, particularly when the algorithm is enhanced by utilising weighting factors such as those disclosed in accordance with the previously mentioned aspect 30 of the invention. Hence differences due eg to surgery can be revealed.

In one embodiment the method for semi-rigid registration comprises:

i) selecting (either manually or automatically) a region of the 'register' surface which 35 appears to be rigid, ii) sampling vertices (possibly randon-fly) from the selected region, and iii) feeding the sampled 'register' vertices together with the entire reference surface to an ICP based registration algorithm (preferably an algorithm in accordance with a previously-mentioned aspect of the invention).

The corresponding rigid points of the 'reference' surface will be identified automatically as the algorithm converges.

It is not necessary to select the rigid region(s) manually or indeed on the basis of any prior knowledge of the surfaces. Instead the algorithm can be tailored to do this by the use of certain weighting factors. In a preferred embodiment points from the register surface are sampled automatically according to a random distribution which 10 is weighted according to the weights dependent on the similarity of local surface neighbourhoods of candidate pairs of corresponding points eg in accordance with the second-mentioned aspect of the invention. The quality of selection will improve as the algorithm converges until in its final stages, only data from rigid regions of the surface are sampled.

In this and other embodiments in order to ensure that the selected region in the fregister' surface is indeed rigid, a 3D distance map (between the 2 surfaces) may be superimposed over one of the surfaces after the registration process is complete. The user may then check whether the points selected were indeed rigid by their projected 20 distance to the other surface. Any errors in selection can then be corrected by resampling only points which are particularly close to the reference surface, for example. This point resampling can be peformed automatically. This process can optionally be performed iteratively to optimise the resulting registration.

0 Preferred embodiments of the invention are described below by way of example only with reference to Figures I to 12 of the accompanying drawings, wherein:

Figure 1 is a diagrammatic representation of the movement of a triangulated register surface to superimpose a region of overlap on the corresponding region of a 30 tirangulated reference surface; Figure 2 is a diagrammatic representation of the nearly registered surfaces of Figure I showing the nearly minimised distances between pairs of corresponding points; Figure 3 is a series of diagrammatic tranverse views of two surfaces as they are registered by an ICP algorithm; Fi ure 4 is a flow diagram of the ICP algorithm; 9 6 Figure 5 is a sketch perspective view of two surfaces being registered by a method in accordance with first aspect of the present invention; Figure 6 shows the search space (in which points projecting to a peak or a ridge may lie) over a convex pyramidal region of the lower surface of Figure 5:

Figure 7 shows a modification in accordance with the first aspect of the invention of the flow diagram of Figure 4; Figure 8 is a diagrammatic transverse view of two surfaces being registered on the basis of a weighting factor based on the angle between normals, in accordance with the second aspect of the invention; Figure 9 is a diagrammatic transverse view of two surfaces being registered on the 15 basis of a weighting factor based on the distance between a pair of candidate corresponding points, in accordance with the second aspect of the invention; Figure 10 is a plot of the variation of the weighting factor of Figure 8 with 0, the angle between the normals from a pair of candidate corresponding points; Figure 11 is a plot of the variation of the weighting factor of Figure 9 with the distance d between a pair of candidate corresponding points, and Figure 12 is a diagrammatic representation of a computer-generated display of 25 registered semi-rigid facial representations of a human face in accordance with the third aspect of the invention.

Figures I and 2 have already been referred to.

Referring to Figure 3a), points Pa to pf are shown on a first surface and points qa to 30 qe are shown on a second surface. The first surface can be regarded as the reference surface and the second can be regarded as the register surface, ie the surface which is moved to register with the reference surface, or vice versa. Each set of points is shown connected (by hairlines) to the closest points of the other set.

These connections can be found by known methods such as those employing a k-d binary tree for example. The k-d binary tree (Henderson 1983, Friedman et al 1977) positions the data on a binary tree in which each point on the tree is connected to its nearest neighbours.

7 The aggregate of the squares of the distances along the above hairlines is minimised in accordance with expression (1) above.

At least one surface (in this case the register surface) is triangulated and, as shown in Figure 3b) after applying the rotation and translation determined from expression I to the set of points qa etc, new points qa', qb' qc' etc are generated by projecting normals (shown arrow-headed) to the facets. The distances along these normals are the distances from points Pa to Pf to the other surface and new points qa', qb' qc' are the closest points to P a to Pf.

The rotation and translation required to minimise the aggregate of the squares of the distances p - q' are found from expression (1) and applied to move the register surface containing points qa', qb' qc' etc towards the reference surface.

The process is repeated iteratively until the aggregate distance (or its mean) is below a predetermined threshold.

This process is summarised in Figure 4, which shows a first stage 100 involving finding the closest points q (on the register surface) to the points p (on the reference surface), a second stage 200 involving minimising the aggregate of the distance pq between the corresponding points to find R and T of expression (1) and a third stage 400 involving generating new candidate corresponding points q' (Figure 3b)) following rotation and translation of the register surface under R and T. When the mean of squared distances is below a predetermined threshold (intermediate step 300) the iteration terminates.

As noted above, we have found that the above iteration often does not converge satisfactorily or even becomes unstable and the reason for this appears to lie in the conventional methods of generating new candidate points by projecting normals to a triangulated surface.

This problem is illustrated in Figure 5. A triangulated register surface is composed of points q1 to q6 (in general, qi) and comprises a convex pyramid (ie pointing towards the other surface) defined by points q4, q5, q6, q7 and q8. A reference surface (not shown triangulated) comprises points pl. to p6.

It is assumed that the ICP algorithm has minimised the aggregate of distances plql, p2q2,.... piqi to find the rotation matrix R and translation vector T of expression (1).

8 In order to find an improved set of candidate closest points ql', q2', q3', q4', q5' q69 etc on the register surface the conventional procedure is to project normals (shown as arrow-headed lines) to a triangular facet surrounding the preceding closest point ql, q2 etc. This works satisfactorily for eg q3' which is generated by the inersection of a normal from p3 to facet qlq3q8. q1', q2' and q4' are also generated correctly by this procedure.

However p6 which lies over a ridge q5q6 cannot be projected normally onto a facet of the register surface and nor can p5 (which lies over the apex q5 of the pyramid referred to above). Accordingly the conventional ICP algorithm utilising point to surface matching would fail to generate q5' and q6' and this is believed to be the reason for the instability of the algorithm in many cases.

In accordance with a feature of the invention, q6' is generated correctly by projecting a perpendicular to the ridge q5q6 and q5' is generated at the same spot as q5 merely by asssociating p5 to the apex q5.

The geometry of the pyramidal portion bounded by q4, q6, q7 and q8 is illustrated in Figure 6. An irregular cross-shaped space defined by the sets of normals from the triangular faces of the pyramid diverges from ridges q4q5, q5q8, q5q7 and q5q6 and extends upwardly towards the reference surface. Any point on the reference surface lying within the central region RI of this space has q5 (the apex of the pyramid) as its closest point. Any point in one of the four regions R2 of this space has a closest point defined by the perpendicular from the underlying ridge to that point.

Accordingly the modified algorithm in accordance with the invention searches in spaces of the type shown in Figure 6 as shown in Figure 7. It should be noted that that points qi surrounded by convex regions made up of five or more triangles will have five or more branches in the regions of space corresponding to R2 in Figure 6 and that in general, sinfilar reasoning will apply to regions of convexity in surfaces having other types of simplex mesh, defining their facets - eg hexagonal or other regular or irregular polygonal facets. What generates the space of Figure 6 is the convexity of the underlying surface - if there is convexity in two orthogonal directions then a region similar to R1 is generated and if there is convexity in one direction only then a region similar to R2 is generated.

This aspect of the invention is not limited to surfaces defined by meshes, but is also applicable to surfaces defined by continuous curves or general surfaces containing discontinuities.

9 Returning to Figure 7, it will be seen that the modified algorithm searches for qi' in a space of type RI ie in a space defined by biconvexity (step 210) and if no such point is found searches in a space of type R2 ie in a space defined by convexity in one direction (step 220). If no such point is found then the algorithm proceeds to step 200' in which the aggregate of the distances pq is minimised as in the conventional algorithm shown in Figure 4. If step 210 results in a found point (which will be an apex) then the existing point q is identified as an apex and is not shifted in the subsequent iteration. If step 220 results in a found point (which will lie over a ridge) then the perpendicular from the ridge to the point p is constructed to locate the new point on the ridge.

A further modification to the ICP algorithm will now be described with reference to Figures 8 to 11 in conjunction with expression (2). A weighting factor wi for use in the expression (2) is generated as the following function of the angle 0 (in radians) between the normals Npi and Nqi of corresponding points pi and qi on the respective surfaces:

wi = exp(-02/a2x2) - exp(-1/4a2) (3) and this function is shown in Figure 10. As a result, during the minimisation of the aggregate of the distances between the sets of points pi, qi, those pairs of points which have a similarly oriented surrounding region (and are correspondingly likely to be true matching points) are driven together preferentially, resulting in a faster convergence of the algorithm. Preferably 0.1:S a:5 0.3.

Another weighting factor wi' can be multiplied or otherwise combined with with wi and is defined as follows with reference to Figure 9:

wi' = 1/(P + d) (4) where d is distance and P > 0 This weighting factor enhances the effect of the points on the surfaces which are already close together in driving the algorithm towards convergence.

Typical applications of the above aspects of the invention are:

a. The aligment of partially overlapping surfaces of the same object, captured from a 3D scanning device.

b. Fusion of multi-modal 3D surface data, e.g. MRI, C7, Surface scan data, etc.

Figure 12 shows a personal computer PC programmed to display a representation of a patient's face on a display D, both before and after facial surgery. In order to show the change in the patients facial appearance the computer is arranged to run a registration algorithm which is optimised to register the rigid (eg forehead regions) of pre- and post-operative 3D optical scans of the patient's face.

The apparatus is operated as follows:

i) a forehead region FH of eg the post-operative surface scan which appears to be rigid is selected ii) vertices from the selected region FH are sampled, and iii) the sampled vertices are fed together with the entire pre-operative surface scan (the latter acting as the reference surface) to an ICP based registration algorithm (preferably an algorithm in accordance with a previously-mentioned aspect of the invention).

The corresponding rigid points of the 'reference' surface are identified automatically as the algorithm converges and the altered areas of the patient's face are displayed selectively eg under the control of a mouse M or other pointing device.

In order to ensure that the selected region in the 'register' surface is indeed rigid, a 3D distance map (between the 2 surfaces) may be superimposed over one of the surfaces after the registration process is complete. The user may then check whether the points selected were indeed rigid by their projected distance to the other surface.

Any errors in selection can then be corrected by resampling only points which are, particularly close to the reference surface, for example. This point resampling can be peformed automatically. This process can optionally be performed iteratively to optimise the resulting registration.

Potential applications of this aspect of the invention include:

1. Aligning pre and post op data according to fixed baseline (as described above); 2. Registering surfaces with outliers (in which the outliers are discarded because they are not related by the common rigid transformation relating the majority of the surface points); 11 3. Extracting best possible plane of symmetry from a surface.

This would involve performing the following steps:

(a) fit a plane through the axis of approximate symmetry of the surface; (b) generate a new surface by reflecting all the vertices of the surface across the plane; (c) select regions on either surface which are reasonably symmetrical; (d) apply semi-rigid registration technique as described above to align the surface (from which symmetrical regions have been identified) with the other surface; (e) find the mid-point between the sets of corresponding vertices (i.e. original and 15 reflected) in the chosen symmetrical regions, and (f) perform a least-squares plane fit through the set of mid points to produce the final plane of symmetry.

The plane of symmetry cam be used as a basis for measuring asymmetry in the surface, particularlyin the case of the surface of a human face. One such application involves measuring the impaction of an orbit in a patient suffering from enophthalmus. The plane of symmetry is used to derive the saggital plane and hence the coronal plane (which is almost parallel with the front of the face). The centre of each orbit is projected onto the coronal plane to find the necessary distances from the centre of each orbit.

12

Claims

1. A method of registering two surfaces with corresponding surface regions by associating corresponding points of the surfaces wherein a first point on one surface facing a region of the other surface is associated with a corresponding point in that region by projecting the first point to a slope discontinuity (eg a ridge or peak) of 10 that region.

2. A method according to claim 1 wherein the slope discontinuity is a ridge.

3. A method according to claim 2 wherein an association between the first point and 15 a perpendicular extending from the ridge to the first point is established.

4. A method according to claim 1 wherein the slope discontinuity is a peak.

5. A method according to claim 4 wherein an association between the first point and 20 the peak is established.

6. A method of iteratively registering two spatial regions with corresponding parts by repeatedly:

7. A method according to claim 6 wherein the spatial regions are surfaces and the local spatial neighbourhoods are surface properties.

8. A method according to claim 7 wherein the weighting factor is inversely related to the angle between the surface normals of the candidate corresponding points.

9. A method according to claim 7 wherein the weighting factor is inversely related to the distance between the pairs of candidate corresponding points.

13

10. A method according to any of claims 7 to 9 wherein the weighting factor is dependent upon the relationship between respective surface curvatures or textures for example.

11. A method according to claim 7 wherein the following function is minimised:

I wi 11 qi' - (Rqi + T) 11 2 (2) where wi is a weighting factor which is dependent on the relationship between local surface neighbourhoods of candidate pairs of corresponding points qi' and qi.

12. A method as claimed in any preceding claim which comprises running a modified ICP algorithm.

13. A method of registering a partially rigid or semi-rigid first surface representation of a subject with a second surface representation of that subject wherein a relatively rigid portion of the first surface representation is selected and used as input to a registration algorithm so as to register with a corresponding portion of that surface representation.

14. A method as claimed in claim 13 as dependent upon any of claims I to 12.

15. A method as claimed in claim 13 or claim 14 which comprises:

i) selecting a region of a register surface which appears to be rigid, ii) sampling vertices from the selected region, and iii) feeding the sampled 'register' vertices together with the entire reference surface to an ICP based registration algorithm.

16 A method as claimed in claim 13 or claim 14 points from a register surface are sampled automatically according to a random distribution which is weighted according to the weights dependent on the similarity of local surface neighbourhoods of candidate pairs of corresponding points.

17. A method of surface registration substantially as described hereinabove with reference to any of Figures 5 to 11 or Figure 12 of the accompanying drawings.