Please verify that (1) all pages are present, (2) all figures are acceptable, (3) all fonts and special characters are correct, and (4) all text and figures fit within the margin lines shown on this review document. Return to your MySPIE ToDo list and approve or disapprove this submission.

On the Use of Lossless Integer Wavelet Transforms in Medical Image Segmentation Nithin Nagaraj, Yogish Mallya Imaging Technologies Lab GE Global Research John F. Welch Technology Center Bangalore 560066, India. ABSTRACT Recent trends in medical image processing involve computationally intensive processing techniques on large data sets, especially for 3D applications such as segmentation, registration, volume rendering etc. Multi-resolution image processing techniques have been used in order to speed-up these methods. However, all well-known techniques currently used in multi-resolution medical image processing rely on using Gaussain-based or other equivalent floating point representations that are lossy and irreversible. In this paper, we study the use of Integer Wavelet Transforms (IWT) to address the issue of lossless representation and reversible reconstruction for such medical image processing applications while still retaining all the benefits which floating-point transforms offer such as high speed and efficient memory usage. In particular, we consider three low-complexity reversible wavelet transforms namely the - Lazywavelet, the Haar wavelet or (1,1) and the S+P transform as against the Gaussian filter for multi-resolution speed-up of an automatic bone removal algorithm for abdomen CT Angiography. Perfect-reconstruction integer wavelet filters have the ability to perfectly recover the original data set at any step in the application. An additional advantage with the reversible wavelet representation is that it is suitable for lossless compression for purposes of storage, archiving and fast retrieval. Given the fact that even a slight loss of information in medical image processing can be detrimental to diagnostic accuracy, IWTs seem to be the ideal choice for multi-resolution based medical image segmentation algorithms. These could also be useful for other medical image processing methods. Keywords: Multi-resolution Analysis, Medical Image Segmentation, Wavelet Transforms, Lossless, Segmentation, Integer Wavelet Transforms, Image Processing.

1. INTRODUCTION The advent of latest technologies for generating high resolution images have brought new opportunities for imaging but has also created challenges for data handling, post-processing, networking and archiving. The enormous data owing to ever reducing slice thickness of the subject image volume may slowdown computer based image processing techniques that are part of the clinical workflow2,11. Modern imaging modalities are capable of generating thinner and more numerous slices with a richer description of the internal anatomy but also invariably resulting in large amounts of image data to be processed and stored. A thoracic CT scan formerly produced approximately 30 sections with 10-mm collimation that was the standard for many years12. The same type of scan, with the 1.25-mm collimation available on state-of-the-art multidetector scanners, now generates about 240 section images for radiologists to interpret. A CTA runoff case with 0.25mm slice spacing covering a longitudinal range of 150cm would generate 6000 images, constituting total memory of 3.0GB on a 512x512 pixel matrix. The increasing volume of image data being acquired owing to these advancements in imaging modalities have motivated researchers to employ multi-resolution (or multiscale) image processing techniques to speed-up the otherwise computationally intensive image processing algorithms. Multi-resolution representation for image processing was first made popular by Burt and Adelson3,4 where they employed the Laplacian and Gaussian pyramids to generate a hierarchical representation of the image at different scales

5747-221 V. 4 (p.1 of 12) / Color: No / Format: A4 / Date: 2005-01-17 03:46:59 SPIE USE: ____ DB Check, ____ Prod Check, Notes:

Please verify that (1) all pages are present, (2) all figures are acceptable, (3) all fonts and special characters are correct, and (4) all text and figures fit within the margin lines shown on this review document. Return to your MySPIE ToDo list and approve or disapprove this submission.

mirroring the way in which processing is done by the Human Visual System. Wavelet based hierarchical representation of data also known as Multi-Resolution Analysis (MRA) is extensively used in most of the present state-of-the-art image compression algorithms. JPEG20005 - the latest international standard for image compression incorporates the Daubechies 9/7 biorthogonal wavelet for lossy compression and it’s integer equivalent for lossless compression. The wavelet representation allows approximate representation of the image at dyadic resolutions with most of the energy being concentrated in the low-resolution approximation images while the high-resolutions images contain the details (edge information). Image compression algorithms exploit this clustering of energy by allocating more bits to the lowresolution coarse images than the high-resolution details. Additionally, it is possible to obtain a resolution embedded bitstream which essentially means that the low-resolution images are reconstructed first and subsequent information enables the reconstruction of a higher resolution image. Such a framework is heavily used in tele-radiology applications6 where such a method is efficient across networks with low-bandwidths since the successive dyadic lowresolution images reduce the waiting time at the decoder, thereby enhancing the productivity of the remote radiologist. Furthermore, it was found that the wavelet representation offers higher compression ratios for the same quality than traditional methods that employ a fixed transform (Discrete Cosine Transform in JPEG). The fact that wavelets can be designed to best suit the data and the application in hand is an added advantage over other fixed transforms (transforms whose kernel does not change with data and/or application). Motivated by the immense success of wavelets for lossless and lossy image compression, other image processing algorithms began to exploit the multi-resolution framework to gain in terms of speed and memory usage. Image registration was one of them, where the mutli-resolution framework gained tremendous popularity. Image registration of two large volumes is essentially an optimization problem where the moving volume is being continuously transformed such that a certain similarity metric between the moving and the fixed volume is minimized (or maximized depending on the objective function). Mutual Information (MI) based image registration is the most popular among these owing to its high accuracy. However, computation of MI and convergence of the metric is a computationally intensive task. A multi-resolution framework was the automatic choice where convergence is first achieved on low-resolution data sets of the volume and the transformation parameters were used as initial conditions so that convergence is achieved faster in subsequent higher resolutions. It is also well known that multi-resolution image registration methods are more robust in terms of convergence than ordinary methods. It may be interesting to note that multi-resolution methods are not as popular in medical image segmentation as they are in registration. Of recent, there are methods that employ multi-resolution methods for image segmentation7 mainly for speed-up. One reason for their relatively less popularity could be that the quality of segmentation affected by such imperfect processing is non-negligible. There is a need for a more “perfect” representation of the image. In this paper, we investigate the use of Integer Wavelet Transforms8 [IWT] for multi-resolution based medical image segmentation. IWTs map integers to integers while obtaining a multi-resolution representation of the image data. The most notable advantage of IWTs over floating-point filters is their ability to perfectly reconstruct the data and their low computational complexity (since they involve only integer arithmetic and low-level operations such as bit-shifts and adds which make them inherently faster). Section 2 talks in detail on these aspects. In Section 3, we consider the example of an Automatic Bone Removal algorithm for Abdomen CTA for removing occluding 3D structure, primarily bone, when trying to analyze vascular regions while using volume rendering methods like Maximum intensity projection (MIP) to demonstrate the utility of these integer transforms. We describe the method adopted for testing these transforms for this application. In particular we consider 3 elementary IWTs namely – Lazy wavelet, the Haar or (1,1) filter and the S+P transform9. In Section 4 we follow up by experiments on two abdomen CTA data-sets and a brief discussion regarding accuracy, speed and implementation and workflow related issues. We summarize our findings and conclude in Section 5.

2. MOTIVATION As was previously mentioned, one of the growing needs of medical image processing is a smart and efficient representation of voluminous medical data to enable complex image processing techniques to perform their tasks in less time and at the same time not compromising on the accuracy, which is paramount for diagnosis. Multi-resolution representation of the image data has been very popular to this end. The idea is to represent the image volume (any dimension) as a pyramid as depicted in Figure 1.0 with the top regions of the pyramid comprising of lower resolution

5747-221 V. 4 (p.2 of 12) / Color: No / Format: A4 / Date: 2005-01-17 03:46:59 SPIE USE: ____ DB Check, ____ Prod Check, Notes:

Please verify that (1) all pages are present, (2) all figures are acceptable, (3) all fonts and special characters are correct, and (4) all text and figures fit within the margin lines shown on this review document. Return to your MySPIE ToDo list and approve or disapprove this submission.

approximations of the original data and the bottom regions of the pyramid represent higher resolution approximations. The base of the pyramid is the original data at its acquired resolution. Image processing algorithms could then operate at any height of the pyramid. It is quite obvious that the volume of the data to be processed decreases exponentially as we traverse from the bottom to the top of the pyramid. By operating at the higher altitude of the pyramid, image processing algorithms could gain significantly in speed at the cost of accuracy of the final output. It should be noted that the computation time (one time) and storage of the entire pyramid would be quite expensive. Therefore in any multiresolution representation, the following features are most desirable: 1.

In-place computation of only a particular resolution (approximation) of the pyramid without computing (and storing) the entire pyramid itself.

2.

Ability to obtain higher or lower resolutions of the data if required from the resolution that was originally chosen.

3.

Ability to accurately obtain the entire image data or smaller Regions of Interest (RoIs) at the highest (original) resolution if required by the image processing algorithm. It is highly desirable that this is obtained without storing the entire pyramid.

4.

Low computational time to obtain the approximate image (low resolution) from the original image and also for reconstructing the original image (or RoI) from the approximate (low resolution) image.

5.

Low memory requirements. Quantitatively, it is reasonable to desire that the storage requirements should not exceed 50% beyond the original volume.

6.

Ability to recovery the original data in a lossless manner at any step as desired by the image processing algorithm.

The above mentioned features would be quite beneficial for any image processing algorithm in general and computationally expensive ones in particular, both from the point of time (speed) and memory. For medical image processing, feature (6) i.e. lossless reconstruction of the data at any step in the image processing algorithm is a very important requirement. Any loss in the data representation and processing could significantly affect the diagnostic accuracy (we do not study the effects of data representation on diagnostic accuracy in this paper).

Figure 1.0. The multi-resolution pyramid. The base of the pyramid represents the original data (highest resolution = N x N). For any height ‘k’ from the base of the pyramid, the resolution is M x M where M = N/k. Typically, ‘k = 2n’ where ‘n’ is a positive integer, in which case the resolutions decrease dyadically.

5747-221 V. 4 (p.3 of 12) / Color: No / Format: A4 / Date: 2005-01-17 03:46:59 SPIE USE: ____ DB Check, ____ Prod Check, Notes:

Please verify that (1) all pages are present, (2) all figures are acceptable, (3) all fonts and special characters are correct, and (4) all text and figures fit within the margin lines shown on this review document. Return to your MySPIE ToDo list and approve or disapprove this submission.

INTEGER WAVELET TRANSFORMS (IWT) We wish to bring to light a transform that can provide all of the aforementioned desirable features and more. The class of transforms we allude to are known as Integer Wavelet Transforms8 (IWT). Since they map integers to integers, they are named thus. They are also known as lossless wavelet transforms, reversible wavelets and perfect reconstruction wavelets in the signal processing/compression literature. IWTs are quite popular in the field of image compression, especially for lossless compression8 applications such as medical image compression. Since they map integers to integers, they do not require any quantization and hence do not suffer from round-off errors unlike their floating-point counterparts. Hence, they are ideally suited for the lossless compression of medical data where FDA regulations severely prohibit any loss for medical data for diagnosis. In this paper, we study three integer wavelet filters. Although there are several IWTs in literature8, our choice was motivated mainly by factors such as ease of implementation and intuitive understanding. We would also like to point that our choice should be looked as a demonstration of the utility of these exciting transforms for medical image processing and not made with any particular engineering trade-off in mind. A prudent selection of a particular IWT should be made depending on the type of image processing application and the nature of the image data. It is also possible to design your own wavelets10 but it is outside the purview of this paper. A. LAZY WAVELET This is the simplest of the IWTs and as the name suggests, it is ‘Lazy’ because it does ‘nothing’. The image data itself is separated into even and odd indices (if the data-length is not even, then it is made even by zero padding). This may be performed in all or any of the dimensions and comprises of one ‘level’ of decomposition. The image represented by the even indices is the low-resolution approximation and the remaining data (the odd indices) is the high-resolution detail. After one level of decomposition, the image size reduces by ½ along that particular dimension. For more levels of decomposition, the same procedure is applied on the low-resolution approximation image (the even indices) that was obtained in the previous level. The action described above of separating the data into odd and even indices is performed for all the dyadic-IWTs and hence we say that the ‘Lazy’ wavelet does nothing beyond this. Other IWTs perform more than just this kind of re-ordering of the data. The equations that describe the Lazy wavelet are as follows: S0[n] = x[2n] D0 [n] = x[2n+1] where x[.] represents the input image samples. S[.] stands for the low-resolution approximation samples and D[.] stands for the details. The subscript for ‘S’ and ‘D’ represent the level of decomposition. In this particular case, S1[n]= S0[n] and D1[n]= D0[n]. B. HAAR (1, 1) The Haar8 IWT, also known as (1, 1) where the number in parenthesis represent the number of vanishing moments for the low pass and high pass filters (the filter which generates the low-resolution approximation is known as the low-pass filter and the filter that generates the details is known as high-pass filter). The Haar IWT is the next simplest filter to the Lazy wavelet. As in the case of the Lazy wavelet, the data is divided into odd and even indices and the equations describing the Haar (1, 1) IWT are as follows: S1[n] =  ½ (S0[n] + D0[n] ) D1[n] = D0[n] - S0[n]

5747-221 V. 4 (p.4 of 12) / Color: No / Format: A4 / Date: 2005-01-17 03:46:59 SPIE USE: ____ DB Check, ____ Prod Check, Notes:

Please verify that (1) all pages are present, (2) all figures are acceptable, (3) all fonts and special characters are correct, and (4) all text and figures fit within the margin lines shown on this review document. Return to your MySPIE ToDo list and approve or disapprove this submission.

C. S+P TRANSFORM As can be guessed, the S+P transform9 is the next improvement over the Haar. S+P has the following equations: D1 (1)[n] = D0[n] - S0[n] S1[n] = S0[n] +  ½ (D1 (1)[n] ) D1[n] = D1 (1)[n] +  α-1(S1[n-2] - S1[n-1]) + α0 (S1[n-1] - S1[n]) + α1 (S1[n] - S1[n+1]) - β1(D1 (1)[n+1] ) where α-1 = 0 , α0 = 2/8, α1 = 3/8 and β1 = -2/8 as suggested by Said and Pearlman9. For the above IWTs, the inverse transformation equations can be trivially deduced from the forward transform equations, and thus are not given here. There are many other IWTs namely (2, 2), (4, 2), (2+2, 2), (6, 2) etc. The reader is referred to the literature8 to study these more complex IWTs. It may be useful to note that one can design one’s own IWT by a method known as ‘Lifting’10. Since our aim is to demonstrate the utility of these transforms for medical image segmentation, we concern ourselves with the simplest of these transforms (Lazy, Haar and S+P). At every level of decomposition of the IWT, the image is divided into what are called as ‘subbands’. The subband structure for a 2D image after one level of decomposition of an IWT is shown in Figure 2.0.a. After two levels of decomposition, six detail subbands and one low-resolution approximation subband is created (Figure 2.0.b). It is important to note that, the calculation of the detail and the approximation subbands are ‘in-place’ which means that the total number of samples is always equal to the acquired resolution. Unlike the traditional multi-resolution pyramid, there is no need to store low-resolution approximations at every level because we could always generate lower resolutions by applying more levels of decomposition and if we wish a higher resolution, we can apply the inverseIWT. All arithmetic is integer-based and it is possible to reconstruct the original image perfectly (lossless). Hence, these are known as lossless, invertible, reversible or perfect reconstruction transforms.

(a)

(b)

Figure 2.0. (a) Left: After one level of decomposition of an IWT. (b) Right: After two levels of decomposition of an IWT. The entire pyramid need not be stored. The details at each level are required for the reconstruction of the low-resolution image at the next higher level. The original image can be ultimately reconstructed in a lossless fashion. Both the details and the approximations are integers.

5747-221 V. 4 (p.5 of 12) / Color: No / Format: A4 / Date: 2005-01-17 03:46:59 SPIE USE: ____ DB Check, ____ Prod Check, Notes:

Please verify that (1) all pages are present, (2) all figures are acceptable, (3) all fonts and special characters are correct, and (4) all text and figures fit within the margin lines shown on this review document. Return to your MySPIE ToDo list and approve or disapprove this submission.

3. INTEGER WAVELET TRANSFORMS FOR MEDICAL IMAGE SEGMENTATION In this section, we study the utility of the three IWTs that were described in the previous section in the context of an Automatic Bone Removal algorithm (AutoBone) for CT-Angiography data (CTA)1 . The Automatic Bone Removal algorithm automatically identifies and removes obstructing anatomical structures like bone from CT-Angiography for detailed 3D visualization of the underlying structures such as vasculature. The approach leverages various observed characteristics of the anatomical region, underlying imaging modality physics, pathology information, and known specification limits of the underlying 3D structure to segment voxels corresponding to obstructing structures. For example, one of the characteristics of the anatomical region is the standard deviation of intensity values for each region. Since the vessel lumen visible under contrast enhancement has relatively uniform density, it is expected to have a relatively lower standard deviation as compared to bone with high intensity on the periphery (Cortical) and low in the interior (Traebecular). Similarly, various other key parameters such as mean intensity, covariance, circularity, size of the region etc. are used in a sequential rule-based classifier to distinguish the bone regions from the vessel regions. The final step in the algorithm is 3D vessel and bone connectivity. The 3D connectivity starts by extracting entire vascular and circulatory regions by considering regions identified as vessel for a constrained region growing to select all vascular and circulatory regions in the image (kidneys, urethra, bladder, etc). The bone mask is obtained by considering all regions classified as bone and 3D region growing into all non-vessel regions in the volume dataset. As indicated in the system block diagram in Figure 3.0, one level of decomposition in the X and Y directions were applied on the input CTA volume (digital format). The low-resolution images were then segmented using the bone-removal algorithm described above. The Maximum Intensity Projection (MIP) of the segmented output is displayed on the screen. The MIP thus obtained is half the original resolution and can be resampled to the original dimensions by 2D interpolation (which is fast). For a fair evaluation of the method, the bone removal algorithm was also applied on the original data (full resolution) followed by the application of one level of IWT to obtain lowresolution images from these segmented high resolution data sets. This formed our gold standard for visual comparison (Figure 4.0).

Figure 3.0. System block diagram describing the method. The IWT is applied along X and Y dimensions (one level) and the Automatic Bone Removal algorithm is applied on the low-resolution image. The MIP of the output is displayed.

Figure 4.0. System block diagram for obtaining the Gold Standard. The Automatic Bone Removal algorithm is applied directly on the CTA image volume at its original resolution and the IWT is applied along X and Y dimensions of the resulting output volume. The MIP of the output is displayed. We show this to compare against the method in Figure 3.0.

5747-221 V. 4 (p.6 of 12) / Color: No / Format: A4 / Date: 2005-01-17 03:46:59 SPIE USE: ____ DB Check, ____ Prod Check, Notes:

Please verify that (1) all pages are present, (2) all figures are acceptable, (3) all fonts and special characters are correct, and (4) all text and figures fit within the margin lines shown on this review document. Return to your MySPIE ToDo list and approve or disapprove this submission.

4. EXPERIMENTAL RESULTS AND DISCUSSION Two data-sets were used in the study – Case1 comprising of 331 slices of 512 x 512 resolution (spacings – 0.703mm x 0.703mm x 1.25mm) and Case2 comprising of 145 slices of 512 x 512 resolution (0.703mm x 0.703mm x 1.25mm). The method was compared with the Gaussian floating point transform. The radius and standard deviation of the Gaussian kernel was chosen to be 1 pixel and 1 grayscale respectively. The Gaussian kernel was applied along the X and Y-axis of the data (followed by downsampling by a factor of two) to facilitate a fair comparison with the IWT. The same bone-removal algorithm was applied on the resulting low-resolution images. Since the Gaussian is a floating-point filter, the filtered output values were truncated to unsigned short data type. The timing results of these experiments are shown in Table. 1.0. Maximum Intensity Projection (MIP) images of the final processed volumes were used for comparison. These are shown in Figure 5.0-5.1. The Gold Standard generated for qualitative comparison is shown in Figure 6.0. Table 1.0. Timing performance of IWT and Gaussian-based multi-resolution methods compared with the Automatic Bone Removal algorithm (AutoBone) on original resolution

Methods AutoBone (Original resolution) AutoBone (Low-resolution) Gaussian + AutoBone (Low-resolution) Lazy IWT+ AutoBone (Low-resolution) Haar IWT+ AutoBone (Low-resolution) S+P IWT+ AutoBone (Low-resolution)

Timing in seconds on a 1.4 GHz, Intel Xeon Processor, 2 GBytes RAM Case1 (512 x 512 x 331) Case2 (512 x 512 x 145) Total Time (secs) secs/slice Total Time (secs) secs/slice 86.00

0.26

35.00

0.24

21.00

0.06

8.00

0.05

55.67

0.17

24.43

0.17

39.43

0.12

16.19

0.11

39.28

0.12

15.86

0.11

49.57

0.15

20.42

0.14

4.1. TIMING PERFORMANCE As it is evident from Table 1.0, the AutoBone algorithm applied on the low-resolution output of the three IWTs perform faster than the algorithm on original resolution (~50% gain) as well as on the Gaussian low-resolution output (23%-30% gain). This is not surprising as the three IWTs are of low computational complexity. The Lazy IWT consists of only re-ordering of data and no arithmetic or logical operations; the Haar IWT involves 2 adds and 1 bit-shift operations (total = 3 operations) and the S+P IWT involves 7 adds, 4 shifts and 1 multiplier (total = 12 operations). This clearly demonstrates the utility of IWT as against floating point transforms like the Gaussian for speed-up of image segmentation algorithms.

5747-221 V. 4 (p.7 of 12) / Color: No / Format: A4 / Date: 2005-01-17 03:46:59 SPIE USE: ____ DB Check, ____ Prod Check, Notes:

Please verify that (1) all pages are present, (2) all figures are acceptable, (3) all fonts and special characters are correct, and (4) all text and figures fit within the margin lines shown on this review document. Return to your MySPIE ToDo list and approve or disapprove this submission.

Original (High Resolution)

(a) Gaussian

(c) Haar (1, 1)

(b) Lazy

(d) S+P

Figure 5.0. Comparison of methods: MIPs for Case1 (512 x 512 x 331).

5747-221 V. 4 (p.8 of 12) / Color: No / Format: A4 / Date: 2005-01-17 03:46:59 SPIE USE: ____ DB Check, ____ Prod Check, Notes:

Please verify that (1) all pages are present, (2) all figures are acceptable, (3) all fonts and special characters are correct, and (4) all text and figures fit within the margin lines shown on this review document. Return to your MySPIE ToDo list and approve or disapprove this submission.

Original (Highest Resolution)

(a) Gaussian

(b) Lazy

(c) Haar (1, 1)

(d) S+P

Figure 5.1. Comparison of methods: MIPs for Case2 (512 x 512 x 145).

5747-221 V. 4 (p.9 of 12) / Color: No / Format: A4 / Date: 2005-01-17 03:46:59 SPIE USE: ____ DB Check, ____ Prod Check, Notes:

Please verify that (1) all pages are present, (2) all figures are acceptable, (3) all fonts and special characters are correct, and (4) all text and figures fit within the margin lines shown on this review document. Return to your MySPIE ToDo list and approve or disapprove this submission.

(a) Gaussian

(c) Haar (1, 1)

(b) Lazy

(d) S+P

Figure 6.0. IWTs on output of Automatic Bone Removal algorithm applied on original highest resolution Case1 (obtained as shown in Figure 4.0). This is used as a Gold standard to compare against the method proposed in this paper. A similar gold standard was generated for Case2 for visual comparison with outputs in Figure 5.0-5.1.

4.2. ACCURACY If we closely observe Figure 5.0 and 5.1, we find that there is a break in the vessels in the output of the Lazy wavelet based method (Figure 5.0.b). The output of the other two IWTs and the Gaussian based methods are comparable in quality to the original AutoBone output on full-resolution image. However, as we noted earlier, the IWTs are faster than the Gaussian based method. Moreover, for the Gaussian, the radius (R) and standard deviation (SD) parameters were chosen after several trials to yield the best output (R = 1 pixel, SD = 1 grayscale). If the standard deviation is a bit high (2 grayscales) the Gaussian based output is of poor subjective quality as can be seen in Figure 7.0. For smaller radius, there appears a cut in one of the main vessels (Figure 7.0.a).

5747-221 V. 4 (p.10 of 12) / Color: No / Format: A4 / Date: 2005-01-17 03:46:59 SPIE USE: ____ DB Check, ____ Prod Check, Notes:

Please verify that (1) all pages are present, (2) all figures are acceptable, (3) all fonts and special characters are correct, and (4) all text and figures fit within the margin lines shown on this review document. Return to your MySPIE ToDo list and approve or disapprove this submission.

(a) Gaussian (R = 0.1, SD = 2.0)

(b) Gaussian (R=2.0, SD = 2.0)

Figure 7.0. Effect of varying the radius (R) of the Gaussian on the AutoBone output with SD = 2.0. For lower radius, there is a break in one of the main vessels (as shown by black arrow in (a) above). For higher radius small vessels are missing (black arrow in (b) above). For a higher standard-deviation some amount of bone is left behind. This shows the sensitivity of parameters for Gaussian.

For a larger SD, the Gaussian output is more homogeneous (smoother) and this leads to a failure of the classifier used in the automatic bone removal algorithm. Since the vessel lumen visible under contrast enhancement has relatively uniform density, it is expected to have a relatively lower SD as compared to bone with high intensity on the periphery (Cortical) and low in the interior (Trabecular). The smoothing operations of the Gaussian changes the characteristics of the regions under consideration and the classifier loses the ability to distinguish the vascular regioins from the bone regions at higher values of SD and radius factor of the Gaussian filter. This is the reason that we see some bone debris scattered in the output (Figure 7.0.b). This sensitivity of parameters for the Gaussian is a problem for segmentation. The IWTs come with no such parameters and this makes them quite robust. We have used only one level of decomposition for multi-resolution based segmentation using IWTs. More levels would imply further speed-up but at the cost of accuracy.

4.2. IMPLEMENTATION, STORAGE AND OTHER ISSUES The IWTs were implemented in C. The AutoBone is implemented as a filter in VTK. Regarding storage requirements, the IWTs are in-place lossless transforms and no additional memory is required for storage during the execution of the algorithm. This is very useful for large medical data-sets such as runoffs which could easily involve 1500 slices. In such cases, creating a complete multi-resolution pyramid is infeasible and in-place transforms are the only solution. With respect to the clinical workflow, the AutoBone is a one-click protocol. With the new proposed IWT based multiresolution segmentation, no change in workflow is needed. This is because, unlike the Gaussian, there are no sensitive parameters to be adjusted for the IWT. The proposed method was found to be consistent and robust with other data-sets as well. It is important to note that IWTs are lossless and if desired, the original data-set at its acquired resolution could be reconstructed in a lossless manner at any step in the algorithm. This makes it ideal for other computationally intensive image processing methods as well. Since IWTs involve integer arithmetic, they are easy to implement in hardware.

5. CONCLUSIONS We summarize our findings as follows: 1. There is a growing need for multi-resolution image processing methods owing to increasing data volumes in medical imaging. Particularly, it is desirable to have a data-representation frame-work which allows for fast

5747-221 V. 4 (p.11 of 12) / Color: No / Format: A4 / Date: 2005-01-17 03:46:59 SPIE USE: ____ DB Check, ____ Prod Check, Notes:

Please verify that (1) all pages are present, (2) all figures are acceptable, (3) all fonts and special characters are correct, and (4) all text and figures fit within the margin lines shown on this review document. Return to your MySPIE ToDo list and approve or disapprove this submission.

2. 3. 4. 5. 6. 7. 8. 9.

execution of otherwise computationally intensive image processing algorithms. Accuracy can’t be compromised and lossless reconstruction is highly desirable for medical image processing. IWTs map integers to integers and hence lossless recovery of the original data is possible both for the entire image and only for a particular RoI (if desired). Their implementation is of low computational complexity, as they comprise mostly of adds and bit-shifts. They could be easily implemented in hardware. We demonstrated their utility in the context of an automatic bone removal algorithm for CTA abdomen datasets. Particularly we investigated the Lazy wavelet, Haar and the S+P transform. A 50% gain in speed was obtained with respect to the AutoBone on original resolution. IWTs were found to be more accurate and faster (23-30%) than the floating-point Gaussian filter. The subjective quality of the segmented output on the low-resolution images of the IWTs was comparable (visually) to those of the Gold Standard. There are several IWTs available and new ones can be designed using the lifting method. The choice of the IWT would depend on the application, algorithms and the type of medical images at hand. The IWT segmented volumes could be losslessly compressed and archived for future retrieval studies.

A comprehensive clinical evaluation on a database of medical images (of different modalities, patients, anatomies, scanners etc.) would be the next appropriate research activity to completely validate the claims made in this paper. However, we hope that this study would motivate the use of IWTs for medical image processing algorithms (not just limited to segmentation) for efficient representation, speed-up, low memory requirements, compact storage and handling; and for obtaining high accuracy by exploiting the power of lossless reconstruction of IWTs which floating point transforms fail to offer.

REFERENCES 1.

R. Mullick, R. Avila, J. Platt, Y. Mallya, J. Knoplioch, F. Poupon, B. Senzig, Automatic Bone Removal for Abdomen CTA: A Clinical Review, Vol 225, pp. 646, RSNA 2002. 2. S. Napel, Dealing with Voluminous MDCT Datasets: What Will It Take? 6th Annual International Symposium Multidetector-Row CT, pp 200-202, June 2004. 3. P. J. Burt, E. H. Adelson, The Laplacian Pyramid as a Compact Image Code, IEEE Transactions on Communications, pp. 532-540, April 1983. 4. E. H. Adelson, C. H. Anderson, J. R. Bergen, P. J. Burt, J. M. Ogden, Pyramid Methods in Image Processing, RCA Engineer, Vol. 29, pp. 33-41, 1984. 5. D. Taubman, M. Marcellin, JPEG2000: Image Compression Fundamentals, Standards, and Practice, Kluwer International Series in Engineering and Computer Science, Kluwer Academic Publishers, Boston, November 2001. 6. N. Nagaraj, S. Mukhopadhyay, F. W. Wheeler, R. S. Avila, Region of Interest and Windowing-based Progressive Medical Image Delivery using JPEG2000, Proceedings of SPIE, Vol. 5033, pp. 382-391, May 2003. 7. B. Kim, J. Shim, D. Park, Fast Image Segmentation Based on Multi-resolution and Wavelets, Pattern Recongnition Letters, Vol. 24, No. 16, pp.2995-3006, December 2003. 8. A. R. Calderbank, I. Daubechies, W. Sweldens, B. L. Yeo, Wavelet Transforms that Map Integers to Integers Journal of Applied and Computational Harmonic Analysis, Vol. 5, No. 3, pp. 332-369, July 1998. 9. A. Said, W. A. Pearlman, An Image Multiresolution Representation for Lossless and Lossy Image Compression, IEEE Transactions on Image Processing, Vol. 5, No. 9, pp. 1303-1310, 1996. 10. W. Sweldens, The Lifting Scheme: A Custom-Design Construction of Biorthogonal Wavelets, Journal of Applied and Computational Harmonic Analysis, Vol. 3, No. 2, pp. 186-200, July 1996. 11. Alexander Lembcke, Mukta Joshi, Dushyant Sahani, Neeraj Rastogi, Sanjeeva Kalva and Sanjay Saini, Workflow Challenges in Multislice Helical CT with Isotropic, High Resolution Imaging, and High-Volume Coverage in the Daily Clinical Routine, 0858PP-e,RSNA 2004. 12. Samuel G Armato III, CAD dissects growing volume of data from lung CT exams, www.diagnosticimaging.com/advancedct2003/

5747-221 V. 4 (p.12 of 12) / Color: No / Format: A4 / Date: 2005-01-17 03:46:59 SPIE USE: ____ DB Check, ____ Prod Check, Notes: