WO2006092594A2

WO2006092594A2 - 3d ultrasound registration

Info

Publication number: WO2006092594A2
Application number: PCT/GB2006/000735
Authority: WO
Inventors: David Hawkes; Graeme Penney; Dean Barratt
Original assignee: Kings College London
Priority date: 2005-03-01
Filing date: 2006-03-01
Publication date: 2006-09-08
Also published as: WO2006092594A3

Abstract

Methods, apparatus and computer program code for registering a 3D ultrasound image with a 3D representation of a body part are described. The fit of the intensity values in the ultrasound image to the representation is optimising using at least one ultrasound scaling parameter as well as position and orientation parameters. The fit of the intensity can be optimised using at least two separate stages. The image and the representation can be converted into probabilistic images where the intensity of the pixels is proportional to the probability of that pixel containing a bone-to-soft tissue interface.

Description

3D ULTRASOUND REGISTRATION

The present invention relates to a computer implemented method of registering a 3D ultrasound image to another 3D image of a body part

Ultrasound imaging techniques can be used to generate 3D images of body parts. Known techniques involve tracking the location and orientation of an ultrasound transducer while a patient is being scanned. Information on the location and shape of the patient's body parts can be derived from the image data and the data from tracking the transducer. Tracking techniques which can be used in relation to the transducer are well known.

The 3D image of a body part (for example a bone) can be created by extracting points from the scanned images that lie on the surface of the body part and calculating their positions relative to a reference coordinate system. The bone surface is then represented in the coordinate system by a sparsely sampled collection of surface points that are derived from the ultrasound scan data. The surface points might be used in a preoperative model using a standard point to surface matching algorithm. The transformation which is calculated in this way provides an estimates of the physical to model transformation.

Localisation errors can arise from ultrasonic effects which can limit the accuracy of the registration of scan data to an image of patient tissue. Such errors can arise due to assumptions which are made in relation to ultrasound imaging techniques, for example concerning the speed with which sound travels through tissue, the path taken by the sound waves between the transmitter, the tissue and the receiver, and the effect of beam thickness. These errors can restrict the accuracy with which images can be created, and prevent the use of ultrasound techniques to image creation for applications such as preoperative planning in orthopaedic and other procedures.

The present invention provides a technique for registering an ultrasound image with a representation of a body part by optimising the fit of intensity values in the ultrasound image. Accordingly, in one aspect, the invention provides a computer implemented method for registering a 3D ultrasound image with a 3D representation of a body part, the method comprising optimising the fit of the intensity values in the ultrasound image to the representation using at least one ultrasound scaling parameter as well as position and orientation parameters.

Advantages arise from the method of the invention in that, by optimising the fit of intensity values, the need to perform a segmentation step on scan data is avoided. This can be a time consuming step in the analysis of scan data. The elimination of a segmentation step can allow scan data to be analysed in preparation for a procedure more quickly, and with less user intervention than has been required previously.

The method of the invention should preferably perform the registration between ultrasound image data and a representation of the volume of the body part. This can be contrasted with existing techniques in which a registration is performed with a representation of the surface of a body part. A further advantage arises from the method of the present invention when compared with existing techniques based on surface analysis techniques is that it can provide enhanced accuracy.

Preferably, the method of the invention involves altering the six rigid body parameters (location and orientation data in six degrees of freedom) and the or each relevant ultrasound scaling parameter so as to optimise the value of a similarity measure between the ultrasound image and the representation of the body part. Preferably, the image and the representation are then converted into probabilistic images where the intensity of the pixels is proportional to the probability of that pixel containing a bone-to-soft tissue interface. The conversion to probabilistic images can be carried out using probability density functions. These functions map two features extracted from the image and the represent- ation (such as the local image gradient) into a probability. The functions are calculated by a process which only needs to be carried out once. This process requires sets of training data relating to the ultrasound image and the representation, in which the bone-to-soft tissue interface has been manually segmented. The probability density function used for converting the ultrasound images can be based on at least one function intended to highlight bone edges in the ultrasound images. For example, the function can be intended to identify image features indicative bone to soft tissue boundary or a bone shadow region of low US intensity. More than one image feature can be used and preferably at least two image features are used. More than two image features can be used.

The probability density function used for converting the representation images can be based on at least one function intended to highlight edges in the representation. The probability density function used for converting the representation images can also be based on a function intended to identify different types of boundaries in the representation, preferably bone to soft tissue boundaries and skin to air boundaries. More than one image feature can be used and preferably at least two image features are used. More than two image features can be used.

Ultrasound scaling parameters can relate to, for example, the speed with which sound travels through tissue, the path taken by the sound waves between the transmitter, the tissue and the receiver, the effect of beam thickness, and ultrasound image scaling (for example in mm.pixel^"1) in the lateral and axial directions relative to the ultrasound beam.

Preferably, the ultrasound image data is subjected to an artefact removal step, to remove regions of the ultrasound image that are caused by artefacts such as shadowing behind boundaries between different body tissue components.

The similarity measure can be calculated as follows. An estimate of the rigid body parameters and the or each relevant ultrasound scaling parameter are used to reslice the probabilistic image derived from the representation in the plane of each ultrasound slice, the pixel values in these reformatted slices and in the probabilistic image derived from the ultrasound images are then compared using the normalised cross-correlation similarity measure. -A-

The algorithm proceeds by altering each of the rigid body parameters and the or each relevant ultrasound scaling parameter by a given step size, and calculating the value of the similarity measure. The algorithm then moves the current estimate of the solution in the direction which produces the greatest improvement in the similarity measure.

Preferably the fit of the intensity is optimised using at least two separate stages, and more preferably at least three separate stages. More than three stages can be used.

The or each optimisation stage can include altering a step size between parameters used during the optimisation stage. The step size can be reduced, and is preferably halved, between iterations of the optimisation during the or each optimisation stage.

The ultrasound scaling parameter or parameters can be used to optimise the fit in less than all of the stages. Preferably, the ultrasound scaling parameters are used to optimise the fit in only one stage, and more preferably in only the last stage.

The optimisation stages can have different resolutions. Preferably the first stage has a lower resolution than subsequent stages. The lower resolution stage or stages can include blurring the representation data.

The 3D representation of the body part can be derived using scanning techniques such as magnetic resonance, X-rays etc. The use of CT scanning techniques is particularly preferred, especially when the method of the invention is practised in relation to a bone part which comprises bone tissue.

Examples of surgical procedures to which the technique of the present invention is applicable include orthopaedic joint procedures (for example replacement of hip, knee, shoulder, ankle and elbow joints), peri-acetabular osteotomy, tibial osteotomy, distal radius osteotomy, anterior cruciate ligament reconstruction, osteoid osteoma excision, bone tumour resection, spinal procedures (for example in the placement of pedicle screws), and fracture surgery According to a further aspect of the invention, there is provided a data processing apparatus for registering a 3D ultrasound image with a 3D representation of a body part, the data processing apparatus including a data processor and a storage device storing computer program instructions which can configure the data processor to optimise the fit of the intensity values in the ultrasound image to the representation using at least one ultrasound scaling parameter as well as position and orientation parameters.

Preferred features of the method indicated above are also preferred counterpart features of the apparatus aspect of the invention.

Generally, embodiments of the present invention employ various processes involving data stored in or transferred through one or more computer systems. Embodiments of the present invention also relate to an apparatus for performing these operations. This apparatus may be specially constructed for the required purposes, or it may be a general- purpose computer selectively activated or reconfigured by a computer program and/or data structure stored in the computer. The processes presented herein are not inherently related to any particular computer or other apparatus. In particular, various general-purpose machines may be used with programs written in accordance with the teachings herein, or it may be more convenient to construct a more specialized apparatus to perform the required method steps. A particular structure for a variety of these machines will appear from the description given below.

According to a further aspects of the invention, there are provided computer program code for carrying out the method of the invention, and any preferred features, and a computer readable medium bearing such computer program code.

In addition, embodiments of the present invention relate to computer readable media or computer program products that include program instructions and/or data (including data structures) for performing various computer-implemented operations. Examples of computer-readable media include, but are not limited to, magnetic media such as hard disks, floppy disks, and magnetic tape; optical media such as CD-ROM disks; magneto- optical media; semiconductor memory devices, and hardware devices that are specially configured to store and perform program instructions, such as read-only memory devices (ROM) and random access memory (RAM). The data and program instructions of this invention may also be embodied on a carrier wave or other transport medium. Examples of program instructions include both machine code, such as produced by a compiler, and files containing higher level code that may be executed by the computer using an interpreter.

Embodiments of the invention will now be described in detail, by way of example only, and with reference to the accompanying drawings, in which:

Figure 1 shows a flow chart illustrating activities relating to a computer assisted surgical procedure in which the invention can be used;

Figure 2 shows an ultrasound probe which can be used to capture ultrasound images used in the invention;

Figure 3 shows a schematic illustration of a surgical setup which can be used to allow the invention to be carried out; Figure 4 shows a schematic block diagram of apparatus including the probe shown in Figure 2 and set up as shown in Figure 3 including data processing apparatus according to the invention;

Figure 5 shows a process flow chart illustrating a first embodiment of a registration method of the invention; Figure 6 shows a process flow chart illustrating a method for creating a probability density function for CT images used in the invention;

Figure 7 shows a process flow chart illustrating a method for creating a probability density function for US images used in the invention;

Figure 8 shows a process flow chart illustrating a second embodiment of the registration method of the invention; and

Figure 9 shows a schematic block diagram of a data processing apparatus according to the invention for carrying out the registration method illustrated in Figures 5 and 8.

Like items in different figures share common reference numerals unless indicated otherwise. With reference to Figure 1, there is shown a flow chart illustrating at a high level various activities related to a surgical procedure 100 in which the registration method of the present invention can be used. The registration method will be described within the context of a surgical procedure being carried out on a femur but the invention is not limited to that particular application.

In a first pre-operative step 101, probability density functions are created for ultra sound (US) and computer tomography (CT) images using US and CT images as training data and which are used subsequently during the registration method to create probabilistic images. In a probabilistic image, the intensity of any pixel, or voxel, is indicative of the probability of that pixel, or voxel, being at a bone to soft tissue interface. It will be appreciated that the probability density functions only need to be created once and can be reused subsequently. It will also be appreciated that the probability density functions can be created entirely independently of the subsequent activities illustrated in Figure 1 and also in a sequence different to that shown in Figure 1. All that is required is that the probability density functions are available for use in the registration step 106. Step 101 will be described in greater detail below with reference to Figures 6 and 7.

In a further pre-operative step 102, images of the patient's body part, in this instance the femur, are captured by carrying out a CT scan of the femur. The digitized CT scan data is then transferred, e.g. over a network, or otherwise made available to the computer of the computer assisted or aided surgery (CAS) system that will carry out the registration method of the invention.

After imaging step 102, a marker 4 trackable by a tracking system part 1 of the CAS system 10 is attached to the patient's body part to allow the position of the patients body part 6 in the reference frame of the tracking system to be dynamically determined. This trackable marker will be referred to herein as a Dynamic Reference Object or DRO. This intra-operative step is not illustrated in Figure 1.

Then at step 104 US images 5 of the patient's body part are taken intra-operatively and the US scan data is processed and made available to the computer part 12 of the CAS system 10. Then at step 106 the registration method of the invention is used to register the patient model 2 in the reference frame of the tracking system of the CAS system. Then, once the model 2 of the patient's body part has been registered the registered patient model can be used to navigate various instruments, implants and other entities during a navigated surgical procedure 108. For example, if a pre-operative surgical plan has been defined relative to the CT model 2, then registration of the CT model will register the surgical plan in the reference frame of the tracking system. The navigated surgical procedure may include various general operations such as planning the surgical procedure or image guided surgical steps with which a person of ordinary skill in the art will be familiar. Subsequent steps have not been described in greater detail so as not to obscure the present invention.

Having given an overview of how the US based registration method of the present invention fits into the overall CAS procedure, the surgical set up and associated apparatus will be described in greater detail with reference to Figures 2, 3 and 4.

Fig. 3 shows a schematic diagram illustrating of the set-up used for the registration method, together with the rigid-body transformations (T) that relate the 3D coordinate systems of the various elements indicated. The model of the patients body part in its frame of reference is illustrated by element 2 and Figure 3 shows the ultrasound probe 8 being used to capture an US image 5 through the patients femur 6.

The apparatus includes a standard clinical US scanner 14 (such as a Philips- ATL HDI- 5000 as provided by Philips Medical Systems, Bothell, WA) with a high-frequency linear- array scan-probe 8 (such as an L12-5 probe with a 5-12MHz broadband transducer) which are used to obtain the US images. The scan-probe can be tracked by a suitable trackable marker 3 such as an array of LEDs, as shown in Figs. 2 and 3. The 3D position of each LED can be measured by the optical localiser 1 part of the tracking system and these positions used to calculate the position and orientation of the marker 3 relative to the dynamic reference object marker 4. Suitable acquisition software, for example written in C++ using Visual Studio 6 (Microsoft Corp., USA), is used to synchronize US image capture with tracking measurements. US images are captured using an analogue-to-digital converter 16 (such as a Canopus ADVC-100, as provided by Canopus UK, Berkshire, UK) connected between the composite video output of the US scanner 14 and the computer 12 of the CAS system. Computer 12 includes or has access to a mass storage device 18 which can store various data required by the CAS system, including the US image data and the pre-operatively captured CT data or other data defining the model 2. Although elements are shown as separate in Figure 4 for the sake of clarity of explanation, it will be appreciated that the elements can be integrated into a single custom unit.

During US scanning, the scan-probe 8 can be covered with a plastic sheath containing a small quantity of US gel to maintain acoustic coupling. In order to minimize motion artefacts and ensure that captured images are closely synchronized with the tracking data, the scan-probe should be swept across the skin surface slowly during step 104.

The rigid-body calibration transformation for the 3D US system, indicated as T_PROBE<__US in Fig. 2, can be determined using a point-target phantom as is known in the art and as described, for example, in "Accuracy of an Electromagnetic Three-Dimensional Ultrasound System for Cartoid Artery Imaging", Ultrasound Med Biol, Vol. 27, pp. 1421- 1425, 2001 which is incorporated herein by reference for all purposes. The calibration procedure involves obtaining tracked images of a pinhead immersed in a water-glycerol solution from many different views. A plurality fo images captured from different directions are used to calculate the calibration transformation. The coordinates of the pinhead in each image are identified manually using custom- written software and the pixel scaling can be determined using electronic calliper measurements provided by the US scanner.

As illustrated by step 104, tracked US images of the surface of the femur are acquired from many different views of the anterior and the posterior femoral shaft, the greater trochanter, the femoral neck, and the epicondyles. US images can also be captured of the pelvis if the pelvis is to be registered instead of the femur. US images can be acquired in 2D compounding mode ('SonoCT™' with the level of compounding set to 'Survey') with a maximum penetration depth of 6cm and a single focal zone. Before images are acquired, the level of the focal zone can be adjusted to approximately correspond to the average depth of bone surface in order to maximise the image resolution at that depth. In compound mode, a number of images are formed by electronically steering the US beam through several preset angles. These images are then combined to produce a single compounded image. Compounding generally reduces speckle artefacts, but has also been found to improve the echo intensity from regions of a curved bone surface not perpendicular to the axial direction of the image.

The optical localiser 1 (such as an Optotrak 3020 as provided by Northern Digital Inc., Ontario, Canada) is used to track the US probe 8 and the dynamic reference object (DRO) 4 which is rigidly attached to the femur 6. The aim is to calculate the registration transformation, T_CT__DR0, which transforms physical positions in the coordinate system or reference frame of the DRO into voxel positions within the preoperative CT volume 2.

To calculate T_CT__DR0, 3D freehand US images are acquired of the bone surface at step 104. These images are then registered to the preoperative CT volume 2 to obtain transformation, T, which maps pixel positions, x_us, in the ith US image to voxel positions, x_cτ, within the CT volume 2. Transformation T is computed from four separate transformations: T = T_CT^DR0{T_QPT^DRQ) T_0PT__PR0BET_pR0BE^us

The calibration matrix, T_PROBE__US, which transforms positions from the US image to positions relative to the infrared-emitting diodes (IREDs) attached to the probe can be calculated as described above by scanning a pin-head in a water-glycerol solution. This calibration consists of eight parameters, six rigid-body parameters (three translations and three rotations) and two scaling parameters, the longitudinal or axial scaling parameter s (which is related to the speed of US) and the transverse or lateral scaling parameter. The transformation T'_OPT__DRO transforms positions relative to the IREDs on the dynamic reference object (DRO) to positions relative to the cameras on the Optotrak localiser. Similarly, the transformation T'_OPT_p_ROBE transforms positions relative to the coordinate system defined by IRED positions on the US probe to the Optotrak coordinate system.

An embodiment of the method of the invention is shown by way of example in the flow chart in Figure 5. An overview of an embodiment of the method of the invention will be described with reference to Figure 5, and a further embodiment of the invention will be described in greater detail with reference to Figures 6 to 8.

As illustrated in Figure 5, the method 50 of the invention calculates the transformation T by altering the six rigid-body parameters which define transformation T_CTVDR0 and the probe calibration scaling parameter in the vertical (axial or y) US image direction, s_y, and horizontal direction (lateral or x) S_x, in order to optimise the value of a similarity measure between the US and CT images. After optimisation, the rigid-body parameters (or the equivalent transformation i.e. T_CT>-DR0 ) ^can be used to determine the position and orientation of any tracked and calibrated object (such as a surgical drill) in the CT scan, and so relate the intraoperative positions of such instruments to a preoperative plan described with reference to the CT image 2.

The method 50 is an intensity based method and so does not require segmentation of the US images or CT volume. Further, the method includes optimising at least one US scaling parameter in order to improve the registration accuracy.

At step 52, the US and CT images are converted into probabilistic images using probability density functions that have previously been created using US and CT training data 54 which was manually segmented. In the probabilistic images created at step 52, the pixel or voxel intensity is proportional to the probability of the pixel or voxel containing a bone to soft tissue interface. At step 54, a current estimate of the six rigid body parameters for T_clv_DRO and the two ultrasound scaling parameters are used to reslice the CT probabilistic image in the plane of each US slice.

Then at step 58, the similarity of the pixel values in the reformatted CT slices and the US slices are compared. The similarity measure for the method is calculated as follows. The current estimate of the rigid-body parameters, s_y and S_x are used to reslice the CT probability image in the plane of each US slice. The pixel values in these reformatted slices and in the US probability images are then compared using a normalised cross- correlation similarity measure (also known as linear correlation, or "Pearson's r"). At step 60, the rigid body parameters and US parameters are incremented or decremented and the similarity measure is re-calculated. The current estimate of the rigid body parameters and the US parameters is then changed at step 62 to values that improve the similarity measure. It is then determined, as illustrated by return loop 64, whether further changing the parameters will increase the similarity measure, i.e. improve the registration. If so, then processing returns to step 58 in which the similarity measure is again calculated with the CT slices reformatted using the current estimate of the rigid body and US scaling parameters.

When it is determined that the similarity measure is not increasing, then the best fit has been found and loop 64 stops and the method ends and outputs data 66 indicating the values for the six rigid body transformation parameters and the values for the two US scaling parameters corresponding to the maximum similarity measure. Hence, registration has been achieved.

The invention will now be described in greater detail with reference to Figures 6 to 8.

Figure 6 shows a flow chart illustrating a method 70 for creating a probability density function from training data in the form of a set of CT images, and which corresponds generally to a part of step 101 of Figure 1. The method 70 allows a function to be created which can convert the pre-operatively captured CT images of the patient into probabilistic images. At step 72 a set of CT images is captured as training data. Then an initial CT image is selected at step 74 and at step 76 the bone is segmented. This can be carried out semi-automatically, for example, by delineating the bone boundary in each transverse slice using the software package Analyze 6.0 as provided by the Mayo Foundation, Rochester, MN. Then at step 78 a binary image of the bone is formed by setting voxels outside the delineated boundary to zero and voxels within the boundary to one. Then at step 80, all those voxels lying on a bone to soft tissue boundary, S_edgeCT, are identified. This can be done by applying vertical and horizontal Sobel gradient operators to each slice of the binary volume. AU the voxles in the training data, S_cτ, are also identified. Then a next CT image from the training data is processed in the same way as indicated by loop 84 until all the training data CT images have been processed. Then at step 86 the probability density function p_cτ(a,b) is determined using:

(_a b) = ^C^^{XG S}edgeCT>fcTl (*) ^{= a} JcT2 (*) ⁼ fil)

where the function C calculates the cardinality (i.e. the number of elements) of the set.

In p_cτ, two image features are used. Both CT features are calculated from two dimensional image slices by applying operators to each CT slice in turn. The first CT feature f_CTi(x) is the intensity of a gradient image G_2D(x). The gradient image is set to be the magnitude of the vertical and horizontal gradient images calculated by convolving each 2D CT slice I_2D(x) with 3 x 3 vertical and horizontal Sobel operators. The second CT feature f_CT2(^x) *^s set equal to the maximum value under a 3 x 3 mask centred on pixel x. The effect of the first feature is to highlight edges within the CT volume, while the second feature helps distinguish between bone to soft tissue edges and skin to air edges.

Figure 7 shows a flow chart illustrating a method 90 for creating a probability density function from training data in the form of a set of US images, and which corresponds generally to a part of step 101 of Figure 1. The method 90 allows a function to be created which can convert the intra-operatively obtained US images of the patient into probabilistic images.

Initially US image data is collected 92 to provide a training set. The captured images are then processed to remove artifacts. Artifact removal uses the fact that most strong reflections cause a loss of signal intensity along the direction of the beam. Artifact removal 94 begins by starting at the bottom of each column of pixels in the US image and moving upward toward the transducer face. The pixels of the image are labelled as artifact until a threshold value of pixel intensity T_art is reached. For example a threshold value of 40, in arbitrary intensity units, can be used. Pixels in a small region immediately adjacent to the transducer face are also labelled as artifacts. For example those pixels in a region within approximately 3mm of the transducer face. This helps to remove artifacts caused by imperfect coupling at ths skin surface. Image regions labelled as artifacts are not used in subsequent image processing. This step effectively produces a mask so that only pixels in the relevant region of the image are further processed.

Then a first US image is selected at step 96 and a number of points are manually picked along the centre of the bone edge at step 98. Then the points are connected by straight lines at step 100. Then at step 102, the closets pixel in each column of pixels to the line defining the bone edge is identified. These pixels define the set S_edgeUS being the pixels defining the bone to soft tissue boundary in the US image. The pixel set S_us is also determined, being all of the non-artifact pixels. Then a next US image from the training set is selected at step 104 and processing loops as indicated by line 106, until S_edgeUS and S_us have been determined for all the training data.

The probability density function for US images, p_us, is then calculated at step 108, using:

_(a b₎ = ^C(fr 1 ^{X 6 S}edgeUS>fuSl (*) ^{= g} *fθS2 (*) ⁼ frP

Paste ^~ _c([φ_{eScτ}fusi (x) = a}j_{us2(χ) = b])}

where again C calculates the cardinality of the set. The US probability density function p_us is also based on two image features. The first US feature, f_usl(x), is the intensity of the US image I_us(x). The second US feature, f_US2(x) is the number of pixels between x and the bottom of the image which were not labelled as an artifact. Both of these features help to highlight bone edges due to the large change in acoustic impedance between soft tissue and bone, which results in the high US intensity values at the boundary (detected by fusi(x)) and a bone shadow region of very low US intensities (detected by f_US2(x))-

Figure 8 shows a process flow chart illustrating a further embodiment of the invention corresponding generally to step 106 of Figure 1. The registration method 110 is generally similar to that shown in Figure 5 above. The computer implemented registration method 110 begins and at step 112 artifacts are removed from each US image slice through the patient in the same manner as described above in connection with creation of the US probability density function. Then at step 114, the non-artifact part of each US image slice is converted from an intensity image I(x) into a corresponding probabilistic image P(x)_us, where x represents the image position, and P(x)us represents the probability of the pixel containing a bone to soft tissue interface. The probability images are calculated using the US probability density function p_us, which, given a set of image features F(x)_us returns an estimate of the probability that position x is a bone edge, i.e. P(X)_Us=P_Us(F(X)_Us)- The probability density function for the US image slices, p_us, is based on two image features F(x)us ⁼ F(fusi(^x)> W(X)) ^as described above.

Then at step 116 the CT image is also converted into a probabilistic image in a similar manner. The CT images are converted from an intensity image I(x) into a corresponding probabilistic image P(x), where x represents the image position, and P(x)_cτ represents the probability of the voxel containing a bone to soft tissue interface. The probability images are calculated using the CT probability density function p_cτ, which, given a set of image features F(x)_cτ returns an estimate of the probability that position x is a bone edge, i.e. P(X)_CT=P_CT(F(X)_CT)- The probability density function for the CT image slices, p_cτ, is based on two image features F(x)_cτ = F(f_CT1(x), f_CT2(x)) as described above.

Then at step 118, a current estimate of the six rigid body parameters and s_y are used to reslice the CT probability image in the plane of each US slice. The method includes three different sets of optimisation parameters which are used in three iterations of the general method. The US scaling parameter s_y is only actually used in the final iteration of the method. After reslicing, the similarity of the resliced CT data and US data is compared at step 120. In particular the pixel values in the reformatted CT slices and in the US probability images are compared using the normalised cross-correlation similarity measure, also known as linear correlation, or Pearson's r.

As indicated above the method uses a three stage optimisation approach. A low resolution stage, a high resolution stage and then a further high resolution stage including optimisation of s_y. The following table summarises the optimisation parameters used for each stage.

Except where specified, all measurements are in mm. The blurring and subsampling apply to the CT volume and σ denotes the standard deviation of a Gaussian blurring kernel. The step sizes indicate the maximum and minimum step sizes, and % change for s_y, used in the method at each stage of the optimisation of the registration. The US image resolution is not altered at each stage, but it is reduced by a factor of four prior to registration by Gaussian blurring and subsampling to give pixel sizes of 0.486 x 0.486mm².

At step 122 it is determined whether the similarity has been maximised, that is whether the fit has been optimised. Initially the fit will not be optimised and so processing proceeds to step 124 at which at least one or some of the parameters are varied and processing returns to steps 118 and 120 and the similarity value is calculated again. The method then updates the current optimum value of the registration parameters by moving in the direction which produced the greatest improvement in the similarity measure. Processing continues to loop until no improvement in the similarity measure for the optimum parameters is found.

Processing then proceeds to step 126 at which, for the current parameter set, the parameter step size is reduced by a factor of two. Processing then returns to step 118 and the proceeding steps are repeated as above. When no improvement in similarity occurs then again at step 126, the parameter step size is reduced by a factor of two and processing returns to step 118. Processing proceeds in this way until the minimum step size has been reached for the parameter set, in the first stage, the low res parameters in the table.

Then when it is determined at step 130 that the minimum step size has been reached, then processing proceeds to step 132 at which the next set of optimisation parameters are used in a second optimisation stage which has a higher resolution than the first optimisation. Steps 118 to 132 are repeated as indicated above until finally a third optimisation is carried out using a third set of optimisation parameters in which the value of s_y is also optimised. After the optimisation has been completed, the output of the method is the six rigid body parameters and US scaling parameter used to register the CT data in the reference frame of the tracking system. Figure 9 illustrates a typical computer system that, when appropriately configured or designed, can serve as the data processing apparatus or computer of the CAS system according to the invention. The data processing apparatus or computer 500 includes any number of processors 502 (also referred to as central processing units, or CPUs) that are coupled to storage devices including primary storage 506 (typically a random access memory, or RAM), primary storage 504 (typically a read only memory, or ROM). CPU . 502 may be of various types including microcontrollers and microprocessors such as programmable devices (e.g., CPLDs and FPGAs) and unprogrammable devices such as gate array ASICs or general purpose microprocessors. As is well known in the art, primary storage 504 acts to transfer data and instructions uni-directionally to the CPU and primary storage 506 is used typically to transfer data and instructions in a bi-directional manner. Both of these primary storage devices may include any suitable computer-readable media such as those described above. A mass storage device 508 is also coupled bi-directionally to CPU 502 and provides additional data storage capacity and may include any of the computer-readable media described above. Mass storage device 508 may be used to store programs, data and the like and is typically a secondary storage medium such as a hard disk. It will be appreciated that the information retained within the mass storage device 508, may, in appropriate cases, be incorporated in standard fashion as part of primary storage 506 as virtual memory. A specific mass storage device such as a CD-ROM 514 may also pass data uni-directionally to the CPU.

CPU 502 is also coupled to an interface 510 that connects to one or more input/output devices such as such as video monitors, track balls, mice, keyboards, microphones, touch- sensitive displays, transducer card readers, magnetic or paper tape readers, tablets, styluses, voice or handwriting recognizers, or other well-known input devices such as, of course, other computers. Finally, CPU 502 optionally may be coupled to an external device such as a database or a computer or telecommunications network using an external connection as shown generally at 512. With such a connection, it is contemplated that the CPU might receive information from the network, or might output information to the network in the course of performing the method steps described herein. Although the above has generally described the present invention according to specific processes and apparatus, the present invention has a much broader range of applicability. In particular, aspects of the present invention is not limited to any particular kind of orthopaedic procedure and can be applied to virtually any method in which US images are sufficient to allow the anatomical structure of interest to be distinguished. Thus, in some embodiments, the techniques of the present invention could be used to register non-rigid anatomical structures, i.e. non-bony structures, such as organs. One of ordinary skill in the art would recognize other variants, modifications and alternatives in light of the foregoing discussion.

Claims

CLAIMS:

1. A computer implemented method for registering a 3D ultrasound image with a 3D representation of a body part, the method comprising optimising the fit of the intensity values in the ultrasound image to the representation using at least one ultrasound scaling parameter as well as position and orientation parameters.

2. A method as claimed in claim 1 , which includes a step of altering the six rigid body parameters and the or each relevant ultrasound scaling parameter so as to optimise the value of a similarity measure between the ultrasound image and the representation of the body part.

3. A method as claimed in claim 2, which includes the subsequent step of converting the image and the representation into probabilistic images where the intensity of the pixels is proportional to the probability of that pixel containing a bone-to-soft tissue interface.

4. A method as claimed in claim 3, in which the conversion step is carried out using probability density functions.

5. A method as claimed in claim 4, wherein the probability density function used for converting the ultrasound images is based on at least one function intended to highlight bone edges in the ultrasound images.

6. A method as claimed in claim 4 or 5, wherein the probability density function used for converting the representation images is based on at least one function intended to highlight edges in the representation.

7. A method as claimed in claim 6, wherein the probability density function used for converting the representation images is also based on a function intended to identify different types of boundaries in the representation.

8. A method as claimed in claim 2, in which the similarity measure is calculated by a method which comprises estimating the rigid body parameters and the or each relevant ultrasound scaling parameter, reslicing the probabilistic image derived from the representation in the plane of each ultrasound slice, and comparing the pixel values in these reformatted slices and in the probabilistic image derived from the ultrasound images using the normalised cross-correlation similarity measure.

9. A method as claimed in claim 1 , in which the ultrasound scaling parameters relate to at least one of the speed with which sound travels through tissue, the path taken by the sound waves between the transmitter, the tissue and the receiver, the effect of beam thickness, and ultrasound image scaling in the lateral and axial directions relative to the ultrasound beam.

10. A method as claimed in claim 1, which includes a step of removing regions of the ultrasound image that are caused by artefacts.

11. A method as claimed in claim 1, wherein the fit of the intensity is optimised using at least two separate stages.

12. A method as claimed in claim 11, wherein the ultrasound scaling parameter is used to optimise the fit in only one of the stages.

13. A method as claimed in claim 11 or 12, wherein the optimisation stages have different resolutions.

14. Computer program code executable by a data processing device to provide the method of claim 1.

15. A computer program product comprising a computer readable medium bearing computer program code as claimed in claim 14.

16. A data processing apparatus for registering a 3D ultrasound image with a 3D representation of a body part, the data processing apparatus including a data processor and a storage device storing computer program instructions which can configure the data processor to optimise the fit of the intensity values in the ultrasound image to the representation using at least one ultrasound scaling parameter as well as position and orientation parameters.