3D long bone reconstruction based on level sets

Computerized Medical Imaging and Graphics 28 (2004) 377–390 www.elsevier.com/locate/compmedimag

3D long bone reconstruction based on level sets S. Morigi*, F. Sgallari Department of Mathematics, University of Bologna, P.zza di Porta San Donato 5, 40127 Bologna, Italy Received 4 January 2004; revised 21 June 2004

Abstract In medical imaging a three-dimensional (3D) object must often be reconstructed from serial cross-sections to aid in the comprehension of the object’s structure as well as to facilitate its automatic manipulation and analysis. The most popular interpolation scheme for a sequence of image slices is the shape-based method, where object information extracted from a given 3D volume image is used in guiding the interpolation process. The paper presents a level set reformulation of the well-known shape-based method as well as a new automatic level set method, which offers better performance. In particular, we focus on X-ray examinations of long bones, which also requires us to deal with the problem of an optimal slice positioning. To this aim, a 2D version of the proposed algorithm will be used to localize a subset of slices from the entire volume image. A number of experiments were performed on computed tomographic real images to evaluate the proposed approach. The experimental results show a substantial improvement of visual effects (qualitative evaluation) using the proposed method in comparison to both the conventional gray-level interpolation scheme and the shape-based method. Compared with the shape-based interpolation scheme the proposed method has much lower computational cost. q 2004 Elsevier Ltd. All rights reserved. Keywords: 3D image interpolation; PDE; Image processing; Level set

1. Introduction Image interpolation is of great importance in biomedical visualization and analysis. Most three-dimesional (3D) biomedical volume images are sampled anisotropically, with the distance between consecutive slices significantly greater than the in-plane pixel size. Image interpolation creates a number of new slices between known slices in order to obtain an isotropic volume image. From this volume image the 3D object of interest is then reconstructed. A major application of multidimensional image interpolation is the 3D reconstruction of bone geometry and tissue density features starting from non-invasive measures, such as computerized tomography. In computerized tomography a bone, or other 3D objects of interest, can be reconstructed from a series of two-dimensional (2D) parallel slices. Situations where these reconstructions are needed include: prosthesis design, surgery planning. It is also used in the 3D * Corresponding author. E-mail addresses: [email protected] (S. Morigi), sgallari@dm. unibo.it (F. Sgallari). 0895-6111/$ - see front matter q 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.compmedimag.2004.07.002

visualization of anatomy for purposes of investigating the structural and functional behaviors of human bones, this being accomplished by finite-element modelling. During in vivo computed tomographic (CT) data acquisition, the number of scans must be limited in order to protect patients from the risks related to X-ray absorption. Since, the 3D reconstruction accuracy depends on which set of image slices are used in the reconstruction process the first problem to cope with is to define the location of a set of image slices which lead to optimal 3D reconstruction. The main goal in this process is then maximize the density of information, minimizing the X-ray absorption. In clinical practice the position of the CT scans is decided by the radiologists on the basis of an X-ray anterior–posterior projection of the bone segment under examination (the socalled scout image). In order to solve this first problem, this work will present a heuristic method to detect a subset of slices that leads to a pseudo-optimal 3D reconstruction of the bone. The second problem is the choice of an image interpolation method to finalize the reconstruction of the missing slices.

378

S. Morigi, F. Sgallari / Computerized Medical Imaging and Graphics 28 (2004) 377–390

In this work, we present a new approach for the image interpolation problem which, while retaining important features of the existing methods, overcomes most of their limitations. Our method is based on the level set ideas developed by Osher and Sethian in Ref. [1], which has been proved powerful in front tracking, and it has recently been adopted in many image processing applications [2,3]. Partial differential equations (PDE) models have recently emerged in image processing as an alternative to traditional approaches [4–6]. The preliminary results we have obtained investigating the use of a PDE approach to the image interpolation and reconstruction problem was promising [7] and this has ‘encouraged’ us to finalize a suitable PDE model for this interesting problem. We adopt the level set technique for representing feature shapes in cross-sections of the object that has to be reconstructed, then we follow the propagation of these shapes between given slices. In general, we do not have any a priori knowledge about the topology of the structures of the shape to be reconstructed. Among other advantages, the level set approach handles naturally the topological changes in the evolving level set curves. Thus, the method can be successfully applied to data set containing tree anatomical structures. A similar philosophy, in a discrete setting, has been followed by the shape-based method, which is one of the most popular among the methods for image interpolation [8]. In this work, we reformulate the shape-based method as an interface propagation problem. This leads us to an intuitive comprehension of the differences with our proposal. The rest of the paper is organized as follows. Section 2 presents an overview of the interpolation methods in medical image processing, while in Section 3 we briefly describe the underlying PDE perspective on moving interfaces. The evolution of the front can be represented both as a boundary value PDE problem, and as a time dependent initial value problem. Then, we give a new formulation of the well known shape-based method for 3D reconstruction in terms of boundary value PDE problem. In Section 4, we propose a new level set evolution approach which efficiently solves the 3D reconstruction problem preserving and enhancing structures. Section 5 describes numerical implementation details of the proposed method and discusses the computational advantages of this approach. In Section 6 we will see how a simplified version of the proposed method can be used as a strategy to improve the selection of the CT scans for the 3D reconstruction. In Sections 7 and 8 we investigate results from a comparison of the performance of the methods on the basis of qualitative and quantitative analysis, according to the common figures of performance merit [9]. The proposed interpolation algorithm has been tested successfully to CT images for a long bone reconstruction case study, although it could be easily applied to different kind of medical images, such as Positron Emission Tomography (PET) and Magnetic Resonance (MR) images.

2. Review of 3D image interpolation methods Traditional interpolation techniques in medical image processing, can be divided into two groups: scene-based and object-based. Scene-based methods determine pixel values for interpolated cross-sections directly from the density values of given slices. In object-based methods, some object information extracted from the slices, such as shape, are used in guiding the interpolation process. The most commonly used scene-based methods are the straightforward nearest-neighbor and linear interpolation of the gray intensity values [10,11]. Higher order techniques such as spline-based image interpolation have been used to fill the inter-slice spaces [12], obtaining slightly superior results. Our experience confirmed the results in Ref. [13]. In fact, we did not find great difference between the performance of gray-level linear and gray-level cubic spline. The latter has been implemented as a non-negative function interpolation problem. Regarding a qualitative comparison, we noticed insignificant differences for small slice spacing, which becomes more pronounced as the slice spacing increases. Of course, the gray-level cubic spline seems to produce smoother boundaries than the gray-level linear approach. However, if the cross-section morphology of individual structures changes or is displaced significantly from slice to slice, image intensities from one tissue type may be interpolated with those of another, producing artifacts in the images and in the resulting 3D reconstruction. Hence, image interpolation methods of this class often do not provide an appropriate correspondence between structures of the same tissue type in adjoining slices, especially near the object boundaries. This problem is accentuated as the distance between adjacent slices increases. Summarizing, this class presents three major disadvantages: low frequency information is introduced, artifacts appear, structures and anatomy of the object tend to be lost. Thus, when objects are extracted and displayed from the interpolated image, these objects often exhibit a ‘staircase’ and generally unsatisfactory appearance. In Section 8 we will show results which better present these limits. Object-based interpolation methods offer an improvement over scene-based interpolation [8,14–17]. Raya and Udupa in Ref. [13] have introduced the shape-based interpolation technique that has become the most relevant of the object-based interpolation class. For segmented binary images, the shape-based interpolation process consists in, first, converting the binary slice into a gray level image, in which the gray value approximates the distance of the pixel to the nearest point on the boundary of the object. This procedure is known as distance transform. The intermediate binary images are estimated by interpolating the distance maps (gray-value images) and thresholding the results at zero. This strategy has been then extended to gray level input data set in Ref. [16], where

S. Morigi, F. Sgallari / Computerized Medical Imaging and Graphics 28 (2004) 377–390

the authors propose a method that is applicable to gray scale images with k gray-level values. The shape-based method interpolates the images incorporating the correspondence (in shape and position) between objects on adjacent slices. Recently, shape-based interpolation methods have been improved using weighted distance transform and morphing [18] with the goal of providing a better approximation of the Euclidian distance for the distance transform step. Methods belonging to the object-based interpolation class preserve the anatomy and the structures of the objects and solve the artifact problems. However, they are extremely time consuming, they produce poor reconstructions of treelike structures, and, finally, they are very sensitive to noise, since a noisy point affects its whole neighborhood, causing images to appear blotchy in noisy areas. Our method can be classified as an object-based method that works directly on the gray intensity values of the objects, and it offers an improvement to the shape-based method. In Section 3 we will present the common PDE framework that allows us to unify our method and the shape-based one, and to outline a clear comparison between the two approaches. For all the experiments (see Section 8), we have considered comparisons with the most representative algorithm of each class of interpolation methods. For the scene-based class we have chosen the gray-level linear interpolation method, while for the object-based class we have considered the shape-based algorithm, using the weighted distance transform.

3. A boundary value PDE formulation for the shape-based method In this section, we present a new formulation for the shape-based method by a boundary value PDE model based on the level set framework. This allows us to make an intuitive comparison between the new proposed method discussed in Section 4 and the shape-based method—arguably the most promising method in literature—both in terms of computational costs, and of numerical solutions. Let us briefly introduce the level set techniques as a mean to numerically approximate the equations of motion for a propagating front. More precisely, given an initial position for an interface G, where G is a closed curve in R2 ; and a speed function F, which gives the speed of G in its normal direction, the level set methods view the moving interface G as the zero level set of a function f(x,tZ0) from R2 to R. Let f(x,tZ0)Z Gd, where d is the distance from xZ(x,y) to G, defined by dðxÞ Z distðx; GÞ dminp2G kx K pk; x Z ðx; yÞ 2R2

(1)

379

and the plus (minus) sign is chosen if the point (x,y) is outside (inside) the initial curve G. The evolution of this level set function f is then linked to the propagation of the front itself through a time dependent initial value problem: IVP :

vf C FjVfj Z 0; vt

fðx; t Z 0Þ given:

(2)

This is the level set equation given by Osher and Sethian [1]; for certain forms of the speed function F one obtains a standard Hamilton–Jacobi equation. Consider now the case where the speed F is always either positive or negative. Let the single valued function T(x,y) be the time at which the curve (the front) crosses each point (x,y), thus the front motion is characterized as the solution to a boundary value problem: BVP : jVTjF Z 1;

Tðx; yÞ Z 0 on G:

(3)

If the speed F depends only on position, then the equation reduces to what is known as the ‘Eikonal’ equation. It can be shown that the distance function (1) is the viscosity solution to Eq. (3). Osher et al. in Ref. [19] provided a link to evolution equation (2) by proving that the t-level set of T(x,y) is the zero level set of the viscosity solution of the evolution equation (2) at time t. The position of the front at time t is given either by the zero level set f(x,y,t)Z0 of the evolving function f or by the level set T(x,y)Zt of the boundary value solution. This set need not be a single curve, and it can break and merge as t advances. In fact, an advantage associated with both these two perspectives on propagating interfaces is that topological changes in the evolving front G are handled naturally. In recent years, level set methods have been used in a variety of settings for problems in computer vision and image processing. The idea behind these approaches is to represent an image as a family of level lines, from which the image can be reconstructed. Let an image I be an intensity map I(x,y) given at each point of a 2D domain. The range of the function I depends on the type of the image; for grayscale images it is a function mapped between 0 and, for example, 255. We define level lines as boundaries of upper level sets, which are defined at each gray-level k, by fðx; yÞ;

s:t: Iðx; yÞR kg:

In Fig. 1 (center), the graph of a test function Iðx; yÞZ 2 2 x eðx Cy Þ ; is shown. Since I(x,y) has negative and positive values, it can be rescaled in RC in order to represent it as a gray-level image (Fig. 1, left). Some of its associated level lines are shown in Fig. 1(right). We should note that this level lines representation is invariant with respect to any increasing contrast change. The basic idea in image interpolation/recovery is to extract particular shapes from that image, where to extract means to produce a mathematical description of these shapes which can be used as initial fronts in the reconstruction process. The iso-intensity contours or level lines representation of an image provides a natural way to obtain this.

380

S. Morigi, F. Sgallari / Computerized Medical Imaging and Graphics 28 (2004) 377–390

Fig. 1. Evolution of the level lines of a test image for different time steps: (a) tZ0, (b) tZ1 and (c) tZ2. The image (left), its representation as graph of a function u(x,y,$) (center), and 10 of its level lines (right).

In the following a new formulation for the shape-based method based on the level set framework using the boundary value PDE model (3) is presented. Suppose we wish to reconstruct a sequence of 2D images (slices). A sub-set of N of these slices is given as input data set, which we refer to as reference slices, while in-between slices have to be reconstructed by the interpolation method. In what follows, we focus on the reconstruction of n slices between two consecutive reference slices U(x,y) and V(x,y), assuming that there is smooth change between the two slices. A small point to note is that the entire data set will be given by applying the interpolation process NK1 times. From a level set point of view, the classical steps of the shape-based method (see Ref. [14]) can be revisited as follows. Step 1 Lifting: the two gray-level images U(x,y) and V(x,y) are represented as a collection of, respectively, p and q level sets.

Step 2 Distance transform: according to Eq. (3), we need to solve the following BVP’s: jVTUðiÞ j Z 1

TUðiÞ ðx; yÞ Z 0

jVTVðjÞ j

TVðjÞ ðx; yÞ

Z1

Z0

on GðiÞ U;

i Z 1; .; p

GðjÞ V ;

j Z 1; .; q

on

(4)

ðjÞ where GðiÞ U and GV are the level lines. Step 3 Interpolation: a linear interpolation process between TUðiÞ and TVðjÞ is then applied to any level set with a common associated gray-value in order to get the n slices between U(x,y) and V(x,y). Each of the gray-level resulting images is then thresholded to the level set associated value, in order to represent it as collections of level lines. Step 4 Collapsing: convert back from the level line representation to the gray scale intensity representation. The most computationally expensive work is obviously done by Step 2, where the distance transform has to be

S. Morigi, F. Sgallari / Computerized Medical Imaging and Graphics 28 (2004) 377–390

computed for each of the pCq images. One of the most popular level set algorithms which solves the problem is the fast marching method (FMM), which requires total O(S log2 S) operations, where S is the number of grid points. Recently, an improvement to O(S) complexity has been introduced [14] by the so-called fast sweeping algorithm.

4. A new level set evolution approach We want to construct an evolution process which moves simultaneously both U(x,y) and V(x,y) with respect to the time t, toward each other, describing by u(x,y,t) and v(x,y,t), a path that goes from u(x,y,0)ZU(x,y) and v(x,y,0)ZV(x,y) toward a common transition image A resulting interpolated sequence of images is produced by the solution of u(x,y,t) forward in time, and the solution of v(x,y,t) backward in time. As we want to provide a shape-based approach to image interpolation that allows us to take into account for inter-pixel relationships, the strategy for the definition of the evolution process should be to move the structures contained in U rather than pixel by pixel. To this aim we use the level set representation of an image and we can see the image interpolation as an evolution of the level set shapes that allows us to account for inter-pixel relationships. Fig. 1 shows u(x,y,t) in three consecutive time steps t of the evolution of the level sets which characterize the image u(x,y,0) in Fig. 1, left. In Fig. 1 each image u(x,y,$) is illustrated together with its representation as a bivariate function. We can now formulate the interpolation process as the evolution of each level set of an image u to become more like the corresponding level sets of the target v, and vice versa. In order to get our final model let us first represent the kth level set of a single image implicitly, by embedding it in a higher dimension function u(x,y,t) which evolves in time according to Eq. (2) vu Z FjVuj vt

uðx; y; 0Þ Z Uðx; yÞ;

381

vðx; y; 0Þ Z Vðx; yÞ:

The term sgn($) defines if a particular point on a given level set is inside or outside the corresponding level set of the target so that the level line can move forward (growing) and backward (shrinking). Since the motion speed of the level sets is not influenced by their structure, in Eq. (6) we take FZsgn($)2 {K1,0,1}, where FZ0 represents the perfect match uZv. We suppose that, under a certain resolution, the shape of the objects represented by their level sets, should change smoothly, and the brightness of the data set does not change significantly, that is the case of CT images here considered. Thus, big variations can be only due to errors in the choice of the scanning sequence. In Ref. [20], model (6) is applied to image blending, and the work presented here takes a similar philosophy. The solutions u and v of Eq. (6) progress toward each other and meet at a single transition image. That is, the differential equation (6) must share a single steady-state solution which forms the midpoint of our image interpolation. For a detailed discussion on this analysis we refer the reader to Ref. [20]. Model (6) can be derived by variational formulation where the distance metric considered is based on the differences between the areas contained within the kth level sets of the two images u and v. This leads to the minimization, over the image domain D, of the area of the kth level set in u and v that does not overlap, that is the sum of the two areas represented by the shaded regions ðin Fig. 2: 3k ðtÞ Z Hðuðx; y; tÞ K kÞ½1 K Hðvðx; y; tÞ K kÞdxdy D

ð C Hðvðx; y; tÞ K kÞ½1 K Hðuðx; y; tÞ K kÞdxdy D

where H(uKk) is the Heaviside function on the values of u that is 1 inside the area (u(x,y,t)Rk) and 0 outside. This reflects the object-based nature of our approach.

(5)

where F is the speed term in the outward normal direction. Eq. (5) describes the time evolution of the level set function u in such a way that the propagating level line is always identified with the kth level set of this evolving function [2,5]. Therefore, the model that establishes the entire evolution process is represented by a coupled pair of PDE given by 8 vu > < Z sgnðv K uÞjVuj vt (6) vv > : Z sgnðu K vÞjVvj vt with Neumann boundary conditions and initial conditions given by

Fig. 2. Areas contained within the kth level set of two images when the shapes have overlap in space. Four subregions are formed by the intersection of level lines.

382

S. Morigi, F. Sgallari / Computerized Medical Imaging and Graphics 28 (2004) 377–390

In a different way, we could follow the most common way to minimize the squared difference (L2 norm) between the gray values in the two images u and v, that is ð 1 3ðtÞ Z ðuðx; y; tÞ K vðx; y; tÞÞ2 dxdy 2 D

1 Z kuðx; y; tÞ K vðx; y; tÞÞk22 2 ( jVujki z

(7)

where Duki and Dvki are discrete approximations in space to the left-hand sides of Eq. (6). The right-hand side of Eq. (6) stems from the hyperbolic dilation/ erosion equation vtuZGjVuj. Consequently, it is advantageous to approximate the gradient by an upwind scheme [1], where the direction of motion of the level set in Eq. (6) is given by the sign of ðuki K vki Þ: Thus, we get 1

jVKujki ¼ ½maxðDKx uki ; 0Þ2 þ minðDþx uki ; 0Þ2 þ maxðDKy uki ; 0Þ2 þ minðDþy uki ; 0Þ2  2 ; 1

jVþ ujki ¼ ½minðDKx uki ; 0Þ2 þ maxðDþx uki ; 0Þ2 þ minðDKy uki ; 0Þ2 þ maxðDþy uki ; 0Þ2  2 ;

if F% 0; if FO 0;

This measure would lead to a scene-based method, with the limitation that we have previously mentioned. The solutions of Eq. (6) are a couple of functions u($,t), v($,t), with t2[0,N[, where

where DCx, DKx, DCy and DKy denote forward, respectively, backward approximations of the spatial derivatives. In the numerical implementation the term sgn(x) is obtain by

 lim uð$; tÞ Z lim vð$; tÞ Z u:

sgnðuki K vki Þ Z

t/N

t/N

Thus, we need to select n slices between U and V, from those produced by the algorithm. Suppose, for the sake of simplicity, that nZ2mC1. The strategy we chose for sampling n interpolated slices considers the following difference metric

for a given s/0. The updates on u and v are therefore given by Duki Z

f ðtÞ dkuðx; y; tÞ K vðx; y; tÞk2 ; with f(0)ZkUKVk2. Since f ðtÞ/t/N 0; we can take the stationary solution u as slice in the middle of the finite sequence, while for the other m, left and right, interpolated slices we choose m values for t such that f ðtÞ Z f ð0Þ K iDf ;

Df Z

f ð0Þ ; m

i Z 1; .; m:

5. Numerical implementation In order to derive a numerical algorithm one has to consider discretizations of space and time. We employ discrete times tkdkt, where k2N, and t denotes the time step size. Moreover, each slice (image) is divided by a uniform mesh of spacing hZ1 into grid nodes. To simplify the notation in the following, a discrete image u(xa,yb), a,bZ1,.S is represented further on by a vector u 2RS!S ; whose components ui,i2{1,.,S!S}, contain the pixel values. For instance, this vector might be constructed by concatenating the rows of an image. Thus uki denotes the approximation of ui at iteration k. Let vki be a similar approximation to v. Hence, using a finite forward difference scheme in time, Eq. (6) are discretized as ( kC1 ui Z uki C tDuki (8) vkC1 Z vki C tDvki i

uki K vki ðs2 C ðuki K vki Þ2 Þ1=2

Dvki ¼

ðs2

vki K uki jVKujki C ðvki K uki Þ2 Þ 1 2

uki K vki jVKvjki ðs2 þ ðuki K vki Þ2 Þ1=2

for F%0. The use of a higher order scheme is not taken into account due to the fact that a propagating level line can develop corners and discontinuities as it advances, since preserving corners may be an important issue, and higher order schemes can smooth them out. Following Refs. [1,20] in order to satisfy stability conditions and avoid the situation where the functions u and v could pass through one another, t in Eq. (8) is required to be t%1/4. However, from our experimentation this constraint is not severe.

6. 2D-CT scanning plan In this section we investigate the problem of defining, in vivo, the optimal scanning sequence for 3D object reconstruction, given a value N representing the finite number of CT slices that can be acquired and a 2D pre-scanning image (e.g. Anterior–Posterior (AP) scout image of the femur bone segment). These N selected slices will be the reference slices for the interpolation method. The detected scanning plan can be defined as optimal if the set of slices automatically chosen by the strategy applied to the 2D pre-scanning image, minimizes the overall error of the 3D reconstruction for a given interpolation method. To this aim the image interpolation technique adopted both

S. Morigi, F. Sgallari / Computerized Medical Imaging and Graphics 28 (2004) 377–390

for the 2D image and for the 3D data set is the level set interpolation method discussed in Section 4 applied along the longitudinal axis. The CT scanning plan strategy we propose provides the selection of a set of N scan locations corresponding to the N slices which minimize the reconstruction error of a given 3D data set made up by M[N slices. According to Ref. [21], to address the bone reconstruction problem we can consider MZ450 which represents a femur longitudinal dimension of 450 mm and NZ50 slices as a reasonable length for the scanning plan. A first step in the scanning plan procedure consists of a 2D reconstruction of the scout using a scanning step given by M/N. Then the global root-mean-square error (RMSE) between the interpolated scout image and the original one is computed. We refer to the RMSE as a function of the scan location x. If we denote by [0,L] the interval along the longitudinal axis where the M slices are located, then for any segment [a,b]3[0,L], the number of slices that we can locate in that segment is given by n[a,b]N, where Ðb RMSEðxÞdx n½a;b Z ÐaL : 0 RMSEðxÞdx For a validation, we compare the CT scanning procedure with the standard radiological protocol and with a uniform distribution of the reference slices along the longitudinal axis. According to standard radiological CT protocol for custom made hip joint prosthesis [21], the position of the CT scans is decided by the radiologist which divides the anterior–posterior scout image into three regions: from the head of the femur to lesser trochanter, the diaphysis, and the distal ephiphysis. The space in-between two consecutive acquired slices in the second region is double with respect to the first and third regions. This protocol is standard for all patients. In Table 1 we report, for three different human femur data sets (F1,F2,F3), the 2D RMSE of the 2D reconstruction with respect to the original 2D data sets using the standard radiological protocol, the uniform scanning plan and the proposed CT scanning procedure. The number of slices selected is NZ50 which is considered, in clinical practice, a safe amount of radiation absorbed by the patient. The methodology based on the 2D scout image was introduced in Ref. [21] where the authors provided an heuristic procedure for detecting a sub-optimal CT scanning plan adopting gray-level linear interpolation. The algorithm turns out to be more computationally expensive than Table 1 2D reconstruction error (RMSE2D) for different CT scanning plans for NZ50 selected slices

F1 F2 F3

Uniform

Radiological

Our

6.78 6.92 5.55

7.43 7.57 6.67

6.13 6.26 3.53

383

Fig. 3. The 2D root mean square error as function of an increasing number of reference slices.

the proposed one, and it reduces the 3D reconstruction error of about 20%. In Fig. 3 the RMSE of the 2D reconstruction of the scout image is drawn for increasing values of the number of reference slices. Fig. 4 shown as the 2D-CT optimal scanning plan improves significantly the 3D reconstruction error of the composite femur, compared with the standard radiological protocol. In fact, the 2D-CT scanning plan procedure reduces the 3D reconstruction error by about 37% with respect to the standard protocol. In Table 5 in Section 8 we will better appreciate the reduction of the error obtained considering our 2D-CT scan strategy with respect to the uniform distribution. In fact, Table 5 in Section 8 reports a measure of the improvement in the reconstruction of the entire femur bone due to the mentioned 2D-CT scan strategy. An explanation of the quantitative measures used in Table 5 will be given in Section 8.

7. Methods of evaluations There are several ways one can objectively compare the proposed method with the methods commonly used, described in Section 2. We first want to discuss the criteria used and the methods of comparison among interpolations. We want to take into account both qualitative as well as quantitative evaluations for assessing interpolations. We considered 3D data sets derived from tomographic data. We intentionally sub-sampled the data set by removing slices from the original isotropic data set and then estimated interpolated slices exactly at the locations corresponding to the removed slices using gray-level linear, shape-based and level set interpolation methods. For the quantitative evaluation we consider the following figures of merit as measures of similarity between

384

S. Morigi, F. Sgallari / Computerized Medical Imaging and Graphics 28 (2004) 377–390

the global root mean square error in intensity over the whole volume (FOM1); the number of sites of disagreement for a given tolerance toll (FOM2); and the largest difference in term of gray-level between the reconstructed and the original data set (FOM3). All the figures of merit are computed on the entire data set. Objective criteria to compare the quality of 3D medical images generated by different techniques have not yet evolved in 3D imaging. Thus, in the qualitative evaluation that follows, the comparison is by subjective observation of the results. For a qualitative evaluation we compared each estimated slice to the original slice, as well as the 3D reconstructed part of the data set to the original one, to appreciate the reconstructed silhouette.

8. Experiments and discussion

Fig. 4. The 3D root mean square error projected on the (x,z) plane by computing the RMSE along the y axis using, respectively: (up) the radiological protocol (down) the 2D-CT scan strategy. The number of reference slices used is 50.

the original volume image I($,$,k) and the interpolated volume image u($,$,k): FOM1 "

IX max JX max Kmax X 1 ðuði;j;kÞKIði;j;kÞÞ2 Z Imax!Jmax!Kmax iZ1 jZ1 kZ1

#1=2

FOM2Zjfði;j;kÞ; s:t:juði;j;kÞKIði;j;kÞj!tollgj FOM3Zmaxi;j;k juði;j;kÞKIði;j;kÞj where Imax!Jmax!Kmax is the total number of points in the original data set, u(i, j, k) is the value of the point (i, j, k) in the estimated data set, and I(i, j, k) is the value of the point (i, j, k) in the original data set. Since it is desirable for the interpolated volume image and the original volume image to have similar intensity distributions, we have considered three different measures of the variation in the histograms of the two volume images:

For quantitative error evaluation, three in vitro CT sequences of a cadaver femur were used, selecting a subset of N reference slices upon the set of M available slices. This would simulate the in vivo investigations where the number of CT scans has to be limited because it determines the quantity of radiation absorbed by the patient [21]. The femurs were positioned with the longitudinal axis parallel to the CT scanning direction and were scanned with a 1 mm step on a GE Sytec 3000 scanner. Each slice is a 512!512 array where each cell has a density value (HU) which is transformed in a correspondent gray value so obtaining a gray-level image. Two of the data sets have slice with pixels of gray-level values in the range [0, 4095], while the third one with pixels of gray-level values in the range [0, 255]. The AP scout images and the corresponding 3D data sets were kindly provided by Prof. Angelo Cappello (D.E.I.S, University of Bologna, Italy) and Ing. Marco Viceconti (Istituti Ortopedici Rizzoli, Bologna, Italy). According to clinical experience we considered in our experiments three different data sets representative, respectively, of three parts of the bone: the proximal epiphysis (data set A), the diaphysis (data set B), and the distal epiphysis (data set C). In the discussion of the experiment results, we use the following abbreviations for these three methods: GL, graylevel linear, SB, shape-based, LS, level set. We express the density of a uniform sub-sampling by a number called interpolation factor, which we define as the slice spacing between two consecutive reference slices in the original data set. 8.1. Example 1 In this first example we apply the interpolation methods GL, SB, and LS with an interpolation factor 8 on different data sets composed of 32 slices from different regions of a femur (regions A,B,C) for a dataset with gray-level values in the range [0, 255]. The results are summarized in Table 2 according to the three figures of merit introduced in

S. Morigi, F. Sgallari / Computerized Medical Imaging and Graphics 28 (2004) 377–390 Table 2 Comparison of different interpolation methods for several data sets from regions A, B and C of the femur bone GL

SB

LS

A

FOM1 FOM2 FOM3

9.14 347,048 150

10.13 286,937 169

8.56 282,300 158

B

FOM1 FOM2 FOM3

1.88 61,716 65

2.16 61,158 77

1.88 59,792 65

FOM1 FOM2 FOM3

9.85 608,887 112

11.29 506,768 127

9.71 504,515 122

C

Section 7. Part B of the bone present no particular change either in topology or in gray-values, or in the structures, thus a simple GL method, at least quantitatively, could be sufficient to reconstruct that part of the bone, but as we will see in the following qualitative examples, the visual effect using GL is not very satisfactory. Then we tested the interpolation methods for increasing values of the interpolation factor (IF) according to the three quantitative FOMs. Table 3 reports the results obtained applying the GL, SB and LS interpolation methods to a dataset of 50 slices in region A and B for a femur dataset with gray-level values in the range [0, 4095]. We observe that LS method produces overall a better estimate of the missing slices than the other two methods. It scores over in all FOM’s except in FOM3. The latter indicates that the largest difference between gray-level values is bigger in the LS method, but this difference or less is present in a fewer number of pixels (see FOM2). A very convenient feature of the LS algorithm presented here is that its complexity is a function only of the size of the data set, as with the scene-based methods. Thus, the computational effort required does not depend either on the number of different gray levels in each slice, or on the complexity of the object shape. Differently, shape-based methods are computationally more expensive because of the Table 3 Comparison of different interpolation methods for several interpolation factor values IF

GL

SB

LS

2

FOM1 FOM2 FOM3

19.98 1,323,373 952

19.52 1,499,186 950

19.60 1,421,491 952

4

FOM1 FOM2 FOM3

46.74 2,947,226 1593

42.32 29,117,048 1616

42.94 2,882,647 1488

8

FOM1 FOM2 FOM3

75.93 3,908,244 2104

69.38 3,833,130 2259

68.92 3,804,959 2064

10

FOM1 FOM2 FOM3

97.64 4,340,122 2372

88.81 4,247,981 2450

91.38 4,179,332 2284

385

Table 4 Computational results (in minute) of the different interpolation methods IF

GL

SB

LS

2 4 8 10

0.072 0.105 0.116 0.124

141 120 115 83

19 14 9 7

calculation of the distance transforms for each gray level in each image. Table 4 is a representative timing table where we reported the execution times (in minutes) of the application of the three interpolation methods to a dataset of 512!512!50 CT slices (16 bit per pixel) on a PC Pentium IV 1700 MHz 512 MB RAM. For the quantitative results of the experiments reported in Table 4 we refer to Table 3. Analogous timing results were obtained with the other two datasets. We can observe that although LS takes longer time than the simple GL interpolation, it is much faster than the SB method. Moreover, the execution times increase for increasing interpolation factors when we apply the GL method because its computational cost depends only on the number of slices that have to be reconstructed. While the computational costs of the SB and LS methods decrease for increasing IF values, that is for a decreasing number of reference slices. In fact, in these cases, the total computational cost depends on how many times we need to solve the PDE problems for SB or LS, respectively, described in Sections 3 and 4, that is on how many couple of reference slices we have to process. Each single PDE problem that reconstructs the images between two given reference slices has a cost which does not depend on the number of slices that have to be produced. 8.2. Example 2 In our second experiment we recover the structure of the entire femur bone. We considered both a uniform sub-sampling with an interpolation factor 9, and a non-uniform sub-sampling using 2D-CT scanning strategy proposed in Section 6. The nonuniform strategy considers 50 reference slices in order to compare the results with the uniform reconstruction with interpolation factor 9, starting from a data set of 420 total slices. The quantitative results are reported in Table 5, while the reconstructed femurs are displayed in Fig. 5, where an isosurface visualization of the femur bone is shown. The original Table 5 Complete data set reconstruction by the different interpolation methods

FOM1 FOM2 FOM3

GL

SB

LS (unif.)

LS (nonunif.)

7.320 3,135,053 191

7.130 2,675,444 190

7.000 3,153,265 254

5.376 2,882,478 253

386

S. Morigi, F. Sgallari / Computerized Medical Imaging and Graphics 28 (2004) 377–390

Fig. 5. 3D reconstruction of the femur bone: (a) original data set; (b) LS; (c) GL; (d) SB.

data set in Fig. 5(a) can be compared with the reconstructed data sets using LS (Fig. 5(b)), GL (Fig. 5(c)), SB (Fig. 5(d)). It can be seen from this example that the SB interpolation gives a ‘cone-like’ connection between boundaries of adjacent reference slices. The problem is that the straight connection is not always the best estimation for the structures to be interpolated. On the other hand GL interpolation, given by Fig. 5(c), resulted in bad staircase effects. This is because GL purely interpolates missing data along the coordinate axes without any structure or object surface information. The results produced by the LS method shown in Fig. 5(b) demonstrate the improved quality of reconstruction.

8.3. Example 3 In the third example, we present the results of an experiment set up to analyse the qualitative accuracy of the interpolation methods in case of topological changes. We selected a sequence of four slices which contain a bifurcation between the first and the last slice in the sequence. Fig. 6 shows the slices of the sequence while the bifurcation is visible in Fig. 7 where the 3D reconstruction is shown. We then compared the iso-surface image of the reconstructed part of the bone, with the iso-surface representation of the original data set. It is clear that both LS

S. Morigi, F. Sgallari / Computerized Medical Imaging and Graphics 28 (2004) 377–390

387

Fig. 6. Slice sequence from data set C: (a) the original slices in anticlockwise order, where (1) and (4) are reference slices; (b), (c), and (d) are the slices (2) and (3) interpolated, respectively, by LS, GL, and SB. Fig. 7. 3D reconstruction of part of data set A: (a) original data set; (b) LS; (c) GL; (d) SB.

and SB interpolation schemes are able to effectively preserve the shape of the object, even if more accurate details can be observed using the LS scheme. In fact, the undesirable effect of chopping part of the bone that the SB method produces is due to the fact that it does not propagate at all every gray-level value which is present in one reference slice, but not in the adjacent reference slice.

While the LS method produces a more natural blend between two adjacent reference slices since it takes into account the shrinking/growing of the corresponding level sets. In case a level set is not present in the next reference slice it smoothly shrinks that level set. Let us investigate this issue with another experiment.

388

S. Morigi, F. Sgallari / Computerized Medical Imaging and Graphics 28 (2004) 377–390

comparison of Fig. 8(b) and (d), shows that, while the LS algorithm gives an accurate reconstruction of the missed part of the bone, the SB produces a poor interpolation of the slices in-between, with the lost of several gray-level details. Finally, GL interpolation (Fig. 8(c)) introduces clear artifacts not present in the original slices (Fig. 8(a)). 8.5. Example 5 In Fig. 9 results are shown from the reconstruction of a part of the femur bone using 10 reference slices sub-sampled from 40 slices in data set C. The reconstructions are displayed by their iso-surfaces obtained by LS interpolation (Fig. 9(b)), by GL interpolation (Fig. 9(c)), and by SB algorithm (Fig. 9(d)). Significantly smoother and less

Fig. 8. Slice sequence from data set C: (a) the original slices in anticlockwise order, where (1) and (4) are reference slices; (b), (c), and (d) are the slices (2) and (3) interpolated, respectively, by LS, GL, and SB.

8.4. Example 4 We now consider an intrinsic critical aspect of the overall interpolation methods that is the presence of an object in a reference slice that does not appear in the next reference slice. We examined the different interpolation methods in such a situation presented in a sequence of four slices belonging to data set C. Fig. 8(a) illustrates the original slices where the slices (1), (4) are reference slices. Visual

Fig. 9. 3D reconstruction of data set C and part of B: (a) original data set; (b) LS; (c) GL; (d) SB.

S. Morigi, F. Sgallari / Computerized Medical Imaging and Graphics 28 (2004) 377–390

‘staircasing’ representations of the femur bone are produced using LS models than those created using GL interpolation. In the 3D reconstructions it is not hard to observe that the SB method leads to unnatural reconstructions (see for example the ‘cone’ effect in Fig. 9(d)), while the LS method improves the silhouette reconstruction. We can make the following observations on the basis of these experiments and a variety of other similar experiments we have carried out: (1) the new LS method produces, in general, an overall better estimate of the missing slices than GL and SB interpolation methods, although there are tests where the differences in the FOMs are subtle and GL or SB result slightly better; (2) comparing the efficiency of the methods, a significant improvement can be appreciated using LS interpolation with respect to the SB method; (3) the qualitative results illustrate how the LS method overcomes problems such as artifacts and staircase effects typical of the GL interpolation, as well as loss of detail and ‘cone effects’ of the SB interpolation; (4) on the other hand, the LS algorithm will not produce desirable results in cases where corresponding objects in the two reference images do not overlap, because it follows a local interpolation. This is also a limit of GL and SB methods. In this case, a preliminary segmentation phase is required and the method should determine a sort of correspondence between the two segmented reference images. We should note that all the CT images we analysed for the experimental tests present a very low noise level as we can observe from the images reported later in this work which did not require a pre-filtered step. Obviously, for processing noisy images, such as ultrasound images, a denoising prefiltering step is essential in order to avoid numerical problems arising from derivative calculations. However, since the PDE level set approaches, as it is well known, are well suited to denoise medical images, a more general model could be developed in order to integrate a denoising operator in the reconstruction model we propose.

9. Conclusions and future plans In this paper, we investigated the problem of 3D reconstruction from cross-sectional images. The approach we presented is based on a PDE model and consists of a pair of coupled, non-linear, first order, PDE’s that model multidimensional level set interpolation. The qualitative and quantitative comparisons with most of the representative standard methods have shown effective improvements for 3D long bone reconstructions starting from 3D CT data set. Thus, we can conclude that using the LS approach we preserve and improve the accuracy of the object-based methods while maintaining a good performance which above all improves the dramatic costs of the shape-based method. Further developments will address the application of the algorithm to medical images with low resolution and with noise. In order to cope with noisy images we want to

389

introduce a regularization term in the model and we will consider the idea to extend the model presented in Ref. [7] for a system of equations. Acknowledgements The authors thank Prof. Angelo Cappello for his support from a bioengineering point of view which is very helpful and important for our current and further research. Special thanks go to Luigi Mattei for assistance with the implementation. References [1] Osher S, Sethian J. Front propagating with curvature dependent speed: algorithms based on the Hamilton–Jacobi formulation. J Comput Phys 1988;79:12–49. [2] Sethian JA. Level set methods and fast marching methods. Evolving interfaces in geometry, fluid mechanics, computer vision, and material science Cambridge monographs on applied and computational mathematics. Cambridge: Cambridge University Press; 1999. [3] Sethian JA. Numerical algorithm for propagating interface: Hamilton– Jacobi equations and conservation laws. J Diff Geom 1990;31:131–61. [4] Weickert J. Anisotropic diffusion in computer vision. Stuttgart: Teubner; 1998. [5] Sapiro G. Geometric partial differential equations and image analysis. Cambridge: Cambridge University Press; 2001. [6] Special issue on partial differential equations and geometry-driven diffusion in image processing and analysis. IEEE Trans Image Process 1998;7(3). [7] Morigi S, Sgallari F, Zannoni C, Cappello A. An anisotropic nonlinear diffusion approach to image interpolation. Proc Conf MTNS 2000, France 2000. [8] Barret WA, Stringham RR. Shape-based interpolation of gray-scale serial slice images. SPIE Proc, 1898 1993;105–15. [9] Grevera GJ, Udupa JK. An objective comparison of interpolation methods. IEEE Trans Med Imaging 1998;17:642–52. [10] Higgins WE, Ledell BE. A non linear filtering approach to the grayscale interpolation of 3D medical images. SPIE Proc Image Proc 1994;2167:284–95. [11] Goshtasby A, Turner DA, Ackerman LV. Matching of tomographic slices for interpolation. IEEE Trans Med Imaging 1992;11:507–16. [12] Herman GT, Rowland SW, Yau MM. A comparative study of the use of linear and modified cubic spline interpolation for image reconstruction. IEEE Trans Nucl Sci 1979;NS-26(2):2879–94. [13] Raya SP, Udupa JK. Shape-based interpolation of multidimensional object. IEEE Trans Med Imaging 1990;9(1):32–42. [14] Herman GT, Zheng J, Bucholtz CA. Shape based interpolation. IEEE Comput Graph Appl 1992;12(3):69–79. [15] Lotufo RA, Herman GT, Udupa JK. Combining shape-based and gray-level interpolations. Visualization in biomedical computing. SPIE Proc Vis Biomed Comput 1992;1808:289–98. [16] Grevera GJ, Udupa JK. Shape-based interpolation of multidimensional gray-level images. IEEE Trans Med Imaging 1996;15(6):881–92. [17] Grevera GJ, Udupa JK. Shape-based interpolation of nD gray scenes. SPIE Proc 2707 Med Imaging 1996;106–16. [18] Bin L, Hancok ER. Slice interpolation using the distance transform and morphing Proceedings of the 13th International Conference on Digital Signal Processing, vol. 2 1997 p.1083–86. [19] Tsai R, Cheng L, Osher S, Zhao H. Fast sweeping algorithms for a class of Hamilton–Jacobi equations. SIAM J Numer Anal 2003;41(2): 673–94. [20] Whitaker RT. A level-set approach to image blending. IEEE Trans Image Process 2000;9(11):1849–61.

390

S. Morigi, F. Sgallari / Computerized Medical Imaging and Graphics 28 (2004) 377–390

[21] Zannoni C, Cappello A, Viceconti M. Optimal CT scanning for long bone 3D reconstruction. IEEE Trans Med Imaging 1998;17:663–6.

Serena Morigi graduated magna cum laude in 1992, with the Laurea in Computer Science at University of Bologna, Italy. In 1997, she received a PhD Degree in Applied Mathematics from the University of Padova, Italy. In 1997, she spent a post-doc fellowship at the Vanderbilt University, Tennessee, USA. She is now Associate Professor of Numerical Analysis at the II Faculty of Engineering, University of Bologna. Her research interests include geometric modelling, surface reconstruction, image processing and computer graphics. She is also interested in discrete ill-posed problems. She serves as reviewer for international journals in numerical analysis and computer graphics. She took part to several national and international research projects.

Fiorella Sgallari graduated magna cum laude in 1976, with the Laurea in Mathematics at University of Bologna, Italy. She is now full professor of Numerical Analysis at the Faculty of Engineering, University of Bologna. Her current research interests are in the field of Numerical Partial Differential Equations in image processing, with particular attention to medical imaging, linear and nonlinear systems of large dimensions and discrete ill-posed problems. She is author of more than 80 papers in national and international journals, and conference proceedings and she currently serves as reviewer for international journals in numerical analysis. She coordinated several national and international research projects and she has been responsible for an Associate Research Center of CNR in Bologna. She was co-organiser of national and international conference and at the moment of ScaleSpace 2005 and ALGORITMY 2005. She is a member of IEEE, SIAM, UMI and GNCS-INDAM.