BACKGROUND OF THE INVENTION In neurosurgery or endovascular surgery, vessel system understanding of the patient is of fundamental importance. Further, accurate curvilinear modelling of anatomical tubular structures may become a crucial step in computer-assisted surgery planning. Curvilinear models of a tubular entity is an abstraction in object representation that retains only those geometric properties that are essential to length, angle and tortuosity measurements while excluding all potentially interfering details such as widths, width changes (e. g. aneurysms) and branching. Previous methods for extracting centrelines (or medial axes) have in many situations not been able to construct an accurate curvilinear model.

A classical method for curvilinear modelling is the skeletonisation technique. This is often achieved by topology-preserving thinning using an algorithm that deletes simple cells after region growing (see N. Nikolaidis and 1. Pitas. 3-D Image Processing Algorithms. John Wiley, 2001).

Alternatively voxels in the middle of a curvilinear structure can be located by smoothing the image using a multi-scale response. A model-based scale-space filtering approach has been also proposed by several authors. The centreline is located by

detecting optimal responses derived from Hessian matrices at a range of scales of interest. This class of algorithms show good results but some local problems occur at junctions or tangent vessels.

Direct tracking from a seed point has also been used to extract the optimal path.

From a given point, the most likely local orientation of the vessel can be searched or predicted. This involves an iterative search, prediction, centering, and image re- sampling. However, the fidelity of the model to the"true"centerline thus extracted depends heavily on the centering algorithm.

Graph-search principles have also been used to calculate the optimal path. To achieve this task, Dijkstra's algorithm has been the most efficient and widely used one for this purpose.

A recent approach proposed by T. Deschamps and L. D. Cohen"Fast extraction of minimal path in 3D images and applications to virtual endoscopy", Medical Image Analysis, 5: 281-299,2001 incorporates multiple improvements over the standard fast marching algorithm, including a technique of using multiple passes of fast marching (from different initial positions) to centre the track.

Unfortunately, using existing centreline extraction algorithms, a high level of accuracy cannot be guaranteed. A common drawback of the existing techniques is they are not sufficiently resistant to neighbourhood interference. Examples of the types of neighbouring interference are shown schematically in Fig. 1 with an initial Vessel 1 having a centreline 2. The vessel and centreline is then shown distorted by a branch 4, another anatomical structure 5, vessel compression 6 and by an aneurysm 7.

Consequently, existing algorithms may introduce spurious tortuosity at places of asymmetricity, e. g. vessel branching, aneurysms, neighbouring objects touching the

object of interest. This may significantly impact the accuracy of any length and tortuosity measurements based on the model.

SUMMARY OF THE INVENTION It is an object of the present invention to provide for an improved method of curvilinear modelling of tubular type structures.

In accordance with a first aspect of the present invention, there is provided a method of constructing a curvilinear model of a tubular structure from a 3-dimensional data set, the method comprising the steps of : (a) forming a substantially minimal cost path for the tubular structure by the steps of : (al) constructing a distance field; (a2) utilising the distance field to construct a propagation map; (a3) utilising the propagation map to construct the substantially minimal cost path; (b) for each point on the substantially minimal cost path, providing the steps of : (bl) sampling the 3-dimensional data set in a substantially locally perpendicular axis to the minimal cost path; (b2) constructing a deformable tube around the substantially minimal cost path determined by alterations in pixel values around the minimal cost path; (b3) extracting an axial measure of the deformable tube as the curvilinear model of the tubular structure. The method can further include the step (b4) of transforming the curvilinear model back into the original data set.

The step (b2) further can comprise the steps of constructing a resample data around the substantially minimal cost path and then constructing the deformable tube within the resample data. The distance field can be constructed utilising a Chamfer algorithm for a distance measure. The distance field can be thresholded. The propagation map can be constructed utilising Fast Marching Method techniques.

The tubular structure can comprise a tubular structure within the human body and the 3-dimensional data set can comprise a scan of a portion of the human body.

BRIEF DESCRIPTION OF THE DRAWINGS Preferred embodiments of the present invention will now be described with reference to the accompanying drawings in which: Fig. 1 illustrates schematically various forms of possible interference present inside the human body; Fig. 2 is a flow chart of the process of extraction of a minimum cost path of the preferred embodiment; Fig. 3 is a flow chart of the steps of formation of a deformable tube of the preferred embodiment; and Fig. 4 illustrates one form of computer system suitable for implementing the preferred embodiment.

DESCRIPTION OF PREFERRED AND OTHER EMBODIMENTS The method of the preferred embodiment for accurate curvilinear modelling substantially reduces the problems described earlier with the existing approaches by using shape knowledge in addition to wave propagation. This is provided by a modified tubular deformable model to implement the constraints.

The preferred embodiment is described with reference to the person skilled in the art of interactive computer graphics design and medical image processing. In particular, it is assumed that the person skilled in the art will be readily familiar with the techniques discussed in papers such as T. Deschamps and L. D. Cohen,"Fast extraction of minimal

path in 3D images and applications to virtual endoscopy"Medical Image Analysis, 5: 281-299,2001.

Fig. 2 illustrates a flow chart in the initial steps of the preferred embodiment in extraction of a suitable minimal cost path. A distance field is first constructed 11 and then uniformly thresholded 12 to obtain suitable propagation channels. A propagation map is then constructed 13 from a start point using Fast Marching Methods (FMM), and then backtracked when the end is reached to generate the minimal cost path (MCP) 14.

Fig. 3 illustrates the subsequent processing steps 20. At each point on the MCP, the original data is resampled 21 in planes perpendicular to the local orientation of the MCP. A deformable tubular model 22 is constructed inside the resampled tube and allowed to deform in order to recover the object of interest. Finally, the curvilinear model is obtained by extracting the spine of the tubular model and transforming it back into the original data. This process ensures symmetric rapid changes, such as a sharp turns, are not lessened while many asymmetric changes, such as a branching event or an aneurysm, can be removed to a large extent. In addition, deforming the model in the transformed data domain vastly simplifies the computation, making it efficient.

Instead of being used to extract the centreline through topology-preserving thinning or ridge enhancement as provided previously, distance transforms are used in the method of the preferred embodiment as a pre-processing step for larger structures to reduce the widths of propagation channels. Doing so achieves two purposes. Firstly, this helps the tubular model in the next stage to maintain faithfulness to the structure of interest, especially at high curvature parts of the structure, and promotes faster convergence of the deformation of the model. This is because better initialization for the tubular model can be achieved with narrower propagation channels relative to the axial curvature. Secondly, higher efficiency can be achieved by reducing the size of

propagation channels, primarily because this reduces the extent of branching and neighbourhood touching-two major causes of propagation leakage-in much the same way as erosion.

This pre-processing is performed by first thresholding the image, then computing a distance map using a Chamfer algorithm, and finally thresholding the distance map.

Using the Chamfer algorithm eliminates the need to administer multiple passes of FMM, as proposed by Deschamps et al and the additional complexity that it brings (T. Deschamps and L. D. Cohen,"Fast extraction of minimal path in 3D images and applications to virtual endoscopy"Medical Image Analysis, 5 : 281-299,2001.

Although the resultant distances are not exactly Euclidean, the approximation was found to be adequate as only the order of distances, rather than the distance values per se, are of interest to the method of the preferred embodiment, and using any of the common Chamfer matrices (e. g. G. Borgefors. Distance transformations in arbitrary dimensions.

Comput. Vision Graphics Image Process, 27: 321-345,1984) that order is preserved on the size of grids for CT or MRI scans is suitable.

Finding a Minimal cost path The minimal cost approach is both effective and robust for extracting curvilinear structures such as the vasculature and the colon. It is robust to many conditions such as vessel stenosis, partial volume effects and noise. Unfortunately, minimizing the path integral of a cost function is an intractable NP-complete problem.

Simulated wave propagation has recently emerged as a desirable approach to finding the minimal cost path. Not only is it efficient, but it can also give sub-grid accuracy. The latter is a distinct advantage over the more traditional graph search algorithms as it overcomes an inherent ambiguity associated with a discrete grid.

The most commonly used method for efficiently tracking a propagating front is the Fast Marching Method (FMM) proposed by Sethian (J. A. Sethian. A fast marching level set method for monotonically advancing fronts. Proc. of the National Academy of Sciences of the USA, 93 (4): 1591-1595). Provided that the speed of the propagation is monotonic, interface evolution is governed by the Eikonal equation |VT|F = I where T is the arrival time of the front, and F its speed. The front speed can be defined as: G * I max (G<, *I) where I is the image, GXJ is a Gaussian kernel and a is a small constant. A small global directional bias term can be optionally added to speed up the fast marching process, but a larger value of impacts the accuracy. One of the purposes of normalising of the intensity is to minimise the need to adjust parameters such as a and fl across different data sets.

Efficient entropy-satisfying numerical methods have been developed by Sethian and his associates based on the Hyperbolic Conservation Laws. After a forward marching phase, backtracking from anywhere in the field in the direction of steepest decent will reach the origin of the propagation. Thus, from a pre-determined end point, a minimal cost path can be found via wave propagation simulated using FMM.

The safest way to avoid spurious tortuosity being introduced, while ensuring that any real tortuosity is not underrepresented by the model, is to segment out any bumps, branches or foreign objects in the neighbourhood. To be able to recognise those, however, one needs to apply an a priori knowledge about the desired and irrelevant structures.

A tubular shape model : Active-surface-based tubular deformable model can be utilised for the recognition. Suitable models can be found in D. Terzopoulos, A. Witkin, and M. Kass. Constraints on deformable models: Recovering 3D shape and nonrigid motion. Artificial Intelligence, 36 (1) : 91-123,1988. The axial constraints of this model means it will reject many sudden changes in the axial direction. Consequently, the irrelevant structures that are filtered out will involve those that include sudden changes in diameter.

Local orientation of the vessel track One commonality among vessels, ducts, bronchi and the colon are that they have circular or elliptical cross sections and smoothly varying radii. To exploit this fact, a tubular shape model is used as a vehicle to carry knowledge about the desired structure for the purpose of filtering out any irrelevant bumps, branches or foreign objects in the neighbourhood. The knowledge is embedded not only in the mesh structure, but also in the internal forces of the model. The combination of the model's intra-ring forces (explained below) and the inflation force alone favours such an outcome that each cross-sectional ring of the resultant tubular model has a nearly constant area. Although the inter-ring forces help enforce this, in the present application they only play a secondary role. The internal forces can be designed so that the constant area favoured by the model approximately corresponds to that of the structure being modelled. The use of a deflation force (see below) near the end of the deformation provides further assurance. The image force localizes the model.

A feature that makes the model simple and computationally efficient is that the model's deformation takes place in a transformed image. As mentioned above, the original data are resampled in planes perpendicular to the local orientation of the MCP.

The resampled data are then stacked up to form a new, transformed data volume. An initial thin tube is constructed in the middle of the new data and allowed to evolve to

minimize the following"energy"functional E of the model in a space of permissible deformation: where v (s, r) = (x (s, r), y (s, r), z (s, r) ) is a parametric surface on a parameter region H, s and r are the parameterisation in the cross-sectional-tangential and axial directions respectively, P (v) is the potential associated with external forces and can be defined as - (V) I where Iis the image. wlo, wol, wl l, w20, w02 control surface properties of tension, rigidity and resistance to twist (they are not necessarily constants). Used in our model, the coefficients wlo and woo encode the strengths of the intra-ring forces mentioned above, while wo, and wo2 represent those of the inter-ring forces. wn is associated with a combination of intra and interring forces. An inflation force is also used throughout the process.

A minimum of E can be reached by solving the associated Euler-Lagrange equation where s and r are the parameterisation in the cross-sectional-tangential and axial directions respectively. Similar methods can be found in"Finite-Element Methods for Active Contour Models and Balloons for 2-D and 3-D Images", by L. Cohen and I.

Cohen, IEEE Transactions on Pattern Analysis and Machine Intelligence, 15 (11): 1131- 1147, Nov 1993.

Using image gradients as the potential, local minima of the functional are usually encountered when the surface passes through some gradients. A technique that can be employed to speed up the model's convergence is that, since the value at any position in the distance field from a binary image corresponds to the smallest distance from the point to the background, it is possible to use the distance value at the points on the MCP, as the initial radii in the shape model described above. This technique promotes a quicker convergence.

After the model converges at such an energy minimum, the medial axis is transformed back to the coordinate system of the original data. This transformed axis is the curvilinear model that the present method produces.

This method provides a high degree of accuracy in length, angle and tortuosity measurements for certain surgical planning (e. g. endovascular) purposes providing all possible measurements concerning a tubular entity, potentially interfering irrelevancies can be filtered out before extracting a model suitable for the measurements. The method provides an implementation of the curvilinear approach and demonstrates that it is resistant to introducing spurious curvatures, those that are not due to any genuine change in the local orientation of the object of interest, while faithfully reproducing"real"high curvatures, those that are actually part of the object in question.

Further improvements can be made. In the implementation presented above, a deformable tubular model and a front propagation approach are the two main elements used in the filtering and the extraction. However, it is possible to modify or replace either of these. For example, incorporating adaptivity in wo, in equation 1, in a fashion such as that expressed in the equation below can help improve the robustness of this approach.

wo, (i, i+1)+α[max(ai,ai+1)-min(ai,ai+1)]+ß, (3) where ai and ai, l are areas of two successive rings of the model. However, experience has suggested that this adaptivity is not strictly necessary.

An a Priori estimation of the vessel should be used to give such a scale of smoothing that it facilitates real edge finding, such as = X where dce is the estimated width of the propagating channel. An inflation force adaptive to the sign and magnitude of the gradient can also be used.

The preferred embodiment is preferably implemented on a standard workstation type computer arrangement utilising standard computer graphics programming languages. A suitable hardware arrangement for programming can be as shown in Fig. 4 with central processing unit 41, memory 42, Disk store 43 and 10 Device Controller 44 arranged around a central bus 45. The arrangement 40 can be suitably programmed in the usual manner to carry out the operations of the preferred embdodiment.

The foregoing describes only preferred embodiments of the present invention.

Modifications, obvious to those skilled in the art, can be made thereto without departing from the scope of the invention.