WO2001031581A1

WO2001031581A1 - Image processing method and apparatus

Info

Publication number: WO2001031581A1
Application number: PCT/GB2000/004113
Authority: WO
Inventors: Ivan Daniel Meir; Norman Ronald Smith
Original assignee: Tct International Plc
Priority date: 1999-10-26
Filing date: 2000-10-25
Publication date: 2001-05-03
Also published as: GB9925343D0; GB2355789A; AU1042001A

Abstract

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 function (2) where wi is a weighting factor which is dependent on the relationship between local surface neighbourhoods of candidate pairs of corresponding points pi and qi. 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

Image processing method and apparatus

Field of the invention

The present invention relates to a method and apparatus for processing 3D datasets, particularly but not exclusively 3D surface representations.

Background

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 1 and 2.

Figure 1 shows a register surface S-^_Q which is shown in the form of a triangulated mesh of points ql, q2,...ql3....qm (in general, q^ 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 be aligned optimally with a reference surface S_REF similarly comprising a simplex mesh of points pi, p2, p3 pl3 (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 S_REG and their corresponding closest points on S_REF. Accordingly the expression:

(1)

is minimised in the ICP algorithm (N being the number of points in the region of overlap and p_; being the closest point on S_REG to q_; ) by assuming that the closest points correspond and then iteratively solving equation (1). Thus the sum of the squares 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 S_REG, the resulting aggregate of distances d_π, d₂₂, d₃₃ .... d_i; being minimised. In practice the sum of 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 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 Nol 15 No 6, December 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 each of 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 approximately 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 Surfaces International Journal of Computer Vision Nol.13, Νo.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.

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 implementation of the ICP algorithm, using point to triangle matching only, can exhibit an unstable behaviour, inaccuracy, or lack of convergence.

Summary of the invention

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 accordance with the invention, a more reliable convergence can be achieved without being perturbed by the occurrence of slope discontinuities such as ridges or peaks. 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: S^b -C^ + r)!² ι=l

(2) where w_; is a weighting factor which is dependent on the relationship between local surface neighbourhoods of candidate pairs of corresponding points p_; and q_;.

It will be noted that the weighting factor utilised by Maurer et al supra was a non- relative weighting factor, ie dependent only on local surface neighbourhood of a point, not 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 post-operatively 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 surf ce representation. For example 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 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 appears to be rigid,

ii) sampling vertices (possibly randomly) 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 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 '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 performed automatically. This process can optionally be performed iteratively to optimise the resulting registration.

Brief description of the drawings

Embodiments of the invention are described below by way of example only with reference to 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 triangulated reference surface; Figure 2 is a diagrammatic representation of the nearly registered surfaces of Figure

1 showing the nearly minimised distances between pairs of corresponding points;

Figure 3 is a series of diagrammatic transverse views of two surfaces as they are registered by an ICP algorithm;

Figure 4 is a flow diagram of the ICP algorithm; 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 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 θ, 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 registered semi-rigid facial representations of a human face in accordance with the third aspect of the invention. Detailed description

Figures 1 and 2 have already been referred to.

Referring to Figure 3a, points p_a to p_f are shown on a first surface and points q_a to q_e 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.

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 1 to the set of points q_a etc, new points q_a', q_b' q_c' etc are generated by projecting normals (shown arrow-headed) to the facets. The distances along these normals are the distances from points p_a to p_f to the other surface and new points q_a', q_b' q_c' are the closest points to p_a to p_f.

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 q_a', q_b' q_c' 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 ql to q6 (in general, q;) 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 pi to p6.

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

(I)-

In order to find an improved set of candidate closest points ql', q2', q3', q4', q5' q6' 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 intersection of a normal from p3 to facet qlq3q8. ql, 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 associating 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 register surface lying within the central region Rl 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 points q- 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, similar 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 Rl 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.

Returning to Figure 7, it will be seen that the modified algorithm searches for qi in a space of type Rl 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 he 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 θ (in radians) between the normals Npi and Nqi of corresponding points pi and qi on the respective surfaces:

wi = exp(-θ²/αV) - exp(-l/4 ²) (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 < a <0.3.

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

wi = l /(β+ d) (4)

where d is distance and β > 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 alignment 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, CT, Surface scan data, etc.

Figure 12 shows a personal computer PC programmed to display a representation of a patients 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 patients 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 improved weighting factor 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 patients 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 performed 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);

3. Extracting best possible plane of symmetry from a surface.

This involves 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 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 can be used as a basis for measuring asymmetry in the surface, particularly in 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.

From the foregoing it will be appreciated that the methods according to the invention can be embodied in a computer program to be run on a computer, for example the computer shown in Figure 12, for manipulating data corresponding to the reference and register surfaces.

Many modifications and variations fall within the scope of the invention. For example whilst in Figure 3, the register surface is shown to be triangulated and the reference surface is moved to be aligned with it, instead the reference surface may be triangulated and the method may be used to move the register surface towards alignment with the reference surface. An advantage of the invention is that when the iteration described with reference to Figures 4 and 7 is carried out. N pairs of corresponding points p, q are available for each step in the iteration. Because each pair can be reliably formed, even when the normal to the triangle does not form the pair, in the event of a slope discontinuity. The consistent provision of N pairs for successive iterations results in a reliable convergence of the iteration process.

Furthermore, whilst the surfaces are described hereinbefore as triangulated they may be segmented in any other appropriate way, with polygons or other shapes.

Whilst the description and claims make reference to a surface and surfaces, it will be understood that the or each surface may be represented by data corresponding to surface points thereon, which may be manipulated by data processing apparatus such as the computer shown in Figure 12, and the following claims are to be interpreted accordingly.

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 of that region.

2. A method according to claim 1 including attempting to find said corresponding point by projecting a line from said first point on said one surface to pass through the other surface at a predetermined angle to the other surface thereby to define said corresponding point thereon, and in the event that such a line cannot be projected from said first point to the other surface, then defining said corresponding point by projecting said first point to said slope discontinuity .

3. A method according to claim 1 or 2 including simulating movement of the surfaces relative to one another such as to minimise the aggregate distance between a plurality of said first points on one of the surfaces and a corresponding associated points on the other of the surfaces.

4. A method according to claim 3 including iteratively performing said simulated movement of the surface until the aggregate distance converges to less than a predetermined threshold.

5. A method according to claim 2 wherein the other surface is defined by a plurality of planar elements and said projecting of the line comprises projecting a normal to one of the planar element that make up the surface.

6. A method according to any preceding claim wherein the slope discontinuity is a ridge.

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

8. A method according to any one of claims 1 to 5 wherein the slope discontinuity is a peak.

9. A method according to claim 8 wherein an association between the first point and the peak is established.

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

i) finding corresponding sets of points in the respective spatial regions,

iii) finding new sets of corresponding points,

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

12. 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.

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

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

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

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

17. 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.

18. A method as claimed in claim 17 as dependent upon any of claims 1 to 16.

19. A method as claimed in claim 17 or claim 18 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.

20. A method as claimed in claim 17 or claim 18 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.

21. A computer program to be run on data processing apparatus to perform a method as claimed in any preceding claim.