1 ISSN LABORATORY OF MATHEMATICS IN IMAGING Fast Sub Voxel Reinitialization of the Distance Map for Level Set Methods Karl Krissian Carl Fredrik Westin N 0001 November 2003 Technical Report Laboratory of Mathematics in Imaging Harvard Medical School Brigham and Women s Hospital department of Radiology 75 Francis Street 02115 Boston MA USA Fast Sub Voxel Reinitialization of the Distance Map for Level Set Methods Karl Krissian Carl Fredrik Westin 1 February 2 2004 1 Harvard Medical School Brigham and Women s Hospital Dep of Radiology Thorn 323 Boston MA 02115 USA Fax 1 617 264 6887 Emails karl westin bwh harvard edu Abstract Redistancing the implicit surface is currently the most time consuming stage in Level Set Methods usually accomplished by applying the Fast Marching algorithm to a binarized image We propose to apply a faster and linear approximation of the Euclidian Distance while maintaining the sub voxel accuracy of the interface Keywords Level Set Methods Euclidean Distance Map Segmentation Contents 1 2 3 4 5 Introduction Sub Voxel reinitialization of the Narrow Band 2 1 Sub Voxel reinitialization 2 2 Comparison with a binary reinitialization Narrow Banded Fast Distance Transform 3 1 Chamfer Distance Transform 3 2 Algorithm Interpretation Experiments 5 1 Accuracy experiments 5 2 Speed experiments 5 3 Segmentation of the white matter 2 2 3 4 5 7 7 9 11 11 12 14 1 Introduction 1 2 Introduction Although the Level Set Method OS88 Set99b has only relatively recently been introduced as an image processing tool it has already prooved useful in delineating contours and in segmenting objects Because it relies on an implicit representation of the evolving contour as the zero crossings of an image it deals naturally with changes in topology However since this method represents contours as the zero crossings of an image the dimensionaly of the problem is increased by one and processing time may increase considerably To offset this problem we typically create a narrow band of the image points lying within a given distance to the contour and then process the Partial Derivative Equation of the contour evolution inside this narrow band The use of the narrow band and numerical considerations require computing the distance to the evolving contour every few iterations Several techniques have been proposed for computing this distance SF99 Set99a In SSO94 the authors propose to solve the following equation in order to compute the signed distance to the interface x 0 0 t sign 0 1 k k 1 In particular in SF99 their formulate a volume preserving constraint designed to prevent the straying of the zero level set from the initial position even after many iteration their demonstrate that they only need to solve the equation up to time t L where L is the thickness of the narrow band about the interface in which a distance function is needed and they use high order methods in time as well as space for solving the equation However this method which has good properties is computationally expensive especially noticeable when dealing with high resolution three dimensional images In RS00 the authors identify the problem of preserving the exact position of the interface and propose to solve a new Partial Derivative Equation for this purpose In Set99a a nlog n algorithm is proposed to compute the distance by propagation until a given distance to the contour is reached In this paper we propose a novel fast algorithm which both preserves the sub voxel accuracy of the contour position and computes its signed narrow banded distance with a linear complexity Add Ross Whitaker Tsikilis Strain Osher 2 Sub Voxel reinitialization of the Narrow Band The Level Set Method used for segmentation is intrinsically a sub voxel estimation of the contours of the object because it represents the object as the zero crossings of a gray scale image The surface is generally extracted at the end of the evolution using the Marching Cubes Algorithm LC87 which hypothesises tri linear interpolation between 2 Sub Voxel reinitialization of the Narrow Band 3 the voxels It is common to reinitialize the Distance Transform from a binary image every few iterations however doing this compromises the sub voxel accuracy 2 1 Sub Voxel reinitialization We will limit our discussion to the isosurface of value zero Hence we designate a voxel neighbor to the isosurface if this voxel is either of intensity zero or if the isosurface crosses one the six segments or four in 2D images with its direct neighbors To calculate this distance for each point designated as neighbor to the surface we could compute the real distance to the interpolated isosurface as the minimum of the distance to the close triangles However such a formula would require extracting the triangles and computing for each triangle its distance to several close voxels a process far too time consuming for our application Alternatively we could estimate the distance to the surface while retaining as much as possible of the same interpolated surface To achieve this we could maintain the relative position of the points of the surface with their neighboring voxels We denote I x y z the intensity of the level set image at the point x y z N3 A point of the surface M when the surface is interpolated by trilinear interpolation is always a linear combination of two voxels V1 x1 y1 z1 and V2 x2 y2 z2 M 1 V1 2 V2 2 2 0 1 and 2 1 1 0 1 where I1 I V1 I x1 y1 z1 and with 1 I2I I 1 I2 I V2 V1 d1 I V2 M 2 d2 1 Figure 1 Distance from the neighbor points to the interpolated surface We denote I the new image with estimates of the distance to the isosurface for voxels neighbor to the isosurface Thus the values of I at V1 and V2 are estimated by projecting 2 Sub Voxel reinitialization of the Narrow Band 4 the vector V1 V2 in the gradient direction I 1 d1 I1 and I V I V2 d2 I2 with I 1 I2 I V 1 V2 k Ik 1 3 Moreover the vector V1 V2 is always oriented in one of the axes of the image co I ordinates and the scalar product V1 V2 k Ik is thereby reduced to taking one of the components of the unit vector in the gradient direction Remarks Owing to the characteristically discrete grid it is not possible to retain the same interpolated surface after the distance estimation in the neighborhood In the simple case shown in fig 1 however the position of the interpolated surface is preserved Because the same voxel can share several edges that intersect with the isosurface we need only capture the minimal distance value estimated