University of Southern California Dissertations and Theses

Random access to compressed volumetric data

(USC Thesis Other)

Random access to compressed volumetric data

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content RANDOM ACCESS TO COMPRESSED VOLUMETRIC DATA by Ayberk Ozturk ______________________________________________________ A Thesis Presented to the FACULTY OF THE VITERBI SCHOOL OF ENGINEERING UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree MASTER OF SCIENCE (ELECTRICAL ENGINEERING) August 2007 Copyright 2007 Ayberk Ozturk ii DEDICATION I would like to dedicate this work to my entire family my father Nurettin Ozturk, my mother Nesrin Ozturk, my brother Burak Ozturk and my girlfriend Martha Northcutt. Without their continuous support and motivation, this work would have not been possible. iii ACKNOWLEDGEMENTS I would like to express my appreciation and thankfulness to Prof. Antonio Ortega for giving me an opportunity to work in this project. His guidance and support helped me to reach my goals and broadened my knowledge of digital image compression. iv TABLE OF CONTENTS DEDICATION………………………………………………………………………....ii ACKNOWLEDGEMENTS……………………………………..…….........................iii LIST OF TABLES………………………………………………..……………….......vi LIST OF FIGURES..…………………………………………………………………vii ABSTRACT…..……………………………………………………………………..viii CHAPTER 1: INTRODUCTION……………………………………...........................1 1.1 Thesis Outline………………………………………………………………….6 CHAPTER 2: WAVELET THEORY........................................................................... 8 2.1 Introduction ……………………………………………………………………8 2.2 Continuous Wavelet Transform (CWT)............................................................ 8 2.3 Discrete Wavelet Transform ............................................................................. 9 2.3.1 Classification of Wavelets................................................................... 10 2.3.2 Wavelet Families................................................................................. 11 2.3.3 Perfect Reconstruction ........................................................................ 12 2.4 The Lifting Scheme Model ............................................................................. 13 2.4.1 Split ..................................................................................................... 14 2.4.2 Predict ................................................................................................. 14 2.4.3 Update ................................................................................................. 15 2.4.4 Inverse Lifting Scheme Wavelet Transform....................................... 16 2.5 Multi-Resolution Analysis .............................................................................. 16 2.6 2-D Separable Wavelet Transform.................................................................. 18 2.7 3-D Separable Wavelet Transform.................................................................. 20 CHAPTER 3: IMAGE CODER: PROGRES.............................................................. 24 3.1 Introduction 24 3.2 Hierarchical Coefficient Dynamic Ranges...................................................... 25 3.3 Coding Algorithm ........................................................................................... 28 3.4 An Example of PROGRES ENCODER ......................................................... 30 CHAPTER 4: DETERMINING WAVELET COEFFICIENTS FOR VISUALIZING A 2-D PLANE………………………………………………………………………..36 4.1 Introduction…………………………………………………………………..36 4.2 Region of Interest............................................................................................ 38 4.3 Discovering the Voxels on the Plane of Region of Interest............................ 39 4.3.1 Algorithm Description ........................................................................ 39 v 4.4 Determining the Wavelet Coefficients............................................................ 42 4.4.1 Extracting Coefficients in Two Dimensional Wavelet Domain ......... 43 4.4.2 Extracting Coefficients in Three Dimensional Wavelet Domain ....... 47 CHAPTER 5: VOLUMETRIC IMAGE CODER ...................................................... 50 5.1 Introduction…………………………………………………………………..50 5.2 3-D PROGRES ............................................................................................... 52 5.3 Block-based Coding........................................................................................ 54 5.3.1 Encoding ............................................................................................. 55 5.4 Designing a Significance Map ........................................................................ 56 CHAPTER 6: 3-D RANDOM ACCESS DECODER ................................................ 60 6.1 Introduction…………………………………………………………………..60 6.2 Decoding wavelet coefficients ........................................................................ 62 6.3 3-D Inverse Haar Wavelet Transform............................................................. 64 6.4 Experiments..................................................................................................... 64 CHAPTER 7: CONCLUSION.................................................................................... 72 7.1 Thesis Summary.......................................................................................... 72 7.2 Future Work ......................................................................................... ...…75 REFERENCES…………………………………………………………………….... 76 APPENDIX……………………….…………………………………………….…....78 vi LIST OF TABLES Table 3-1: PROGRES encoder algorithm.................................................................... 29 Table 3-2: Encoding of wavelet coefficients by PROGRES-1 .................................... 32 Table 3-3: Encoding of wavelet coefficients by PROGRES-2 .................................... 33 Table 3-4: Encoding of wavelet coefficients by PROGRES-3 .................................... 34 Table 3-5: Encoding of wavelet coefficients by PROGRES-4 .................................... 35 Table 4-1: Steps of determining the voxels on a given 2-D plane............................... 40 Table 4-2: Steps of locating coefficients in two dimensional wavelet domain............ 45 Table 4-3: From 4-2a through 4-2e, wavelet coefficients for reconstructing each block of pixels are shown. 4-2f shows combined wavelet coefficients................................. 46 Table 4-4: The wavelet coefficients needed for reconstruction are shown in Table 4-5b through 4-5e. ................................................................................................................ 49 Table 6-1: Lossless Compression Performances (Bit/Voxel) ...................................... 65 Table 6-2: Random access decoding of diagonal cuts with 45° and 135° from resolution 256x164 to 224x164.................................................................................... 68 Table 6-3: Random access decoding of diagonal cuts with 45° and 135° from resolution 192x164 to 160x164.................................................................................... 69 Table 6-4: Random access decoding of diagonal cuts with 45° and 135° from resolution 128x164 to 96x164...................................................................................... 70 Table 6-5: Random access decoding of diagonal cuts with 45° and 135° with resolution 64x164......................................................................................................... 71 vii LIST OF FIGURES Figure 2-1: 2 Level Lifting Scheme Diagram................................................................... 13 Figure 2-2: 1-D Wavelet Transform ................................................................................. 17 Figure 2-3: The filtering steps of 2-D separable wavelet transform ................................. 18 Figure 2-4: 2-D Forward and Inverse Wavelet Transform ............................................... 19 Figure 2-5: 3-D Wavelet Subbands................................................................................... 21 Figure 2-6: 3-D Forward and Inverse Wavelet Transform ............................................... 23 Figure 3-1: Coefficients of trees ....................................................................................... 26 Figure 3-2: Scanning order of coefficients........................................................................ 28 Figure 3-3: Wavelet Coefficients of a 16x16 block.......................................................... 30 Figure 4-1: A 2-D diagonal plane cuts the volumetric data.............................................. 40 Figure 4-2: Blocks of pixels which are reconstructed for visualizing the line.................. 45 Figure 4-3: Generating a 3-D volume data from combining multiple slices .................... 47 Figure 5-1: 3-D Spatial Orientation Tree.......................................................................... 54 Figure 5-2: Wavelet coefficients of a first slice of a 16x16x16 block.............................. 56 Figure 5-3: A significance map representing 16x16x2 coefficients of a unit block ......... 57 Figure 6-1: A significance map representing 16x16x2 coefficients of a unit block and blue colored dynamic range value includes the seeking coefficient at (9,4,6) ................. 63 Figure 6-2: The 2x2x2 subblock of blue colored dynamic range value in Figure 6-1...... 64 Figure 6-3: Diagonal cuts through the volume data with an angle of 0°, 45°, 90°, 135°. 66 Figure 6-4: Lossless bit rate (bit/voxel) for 2-D diagonal planes with different angles ... 67 Figure 6-5: Lossless bit rate (bit/voxel) and number of decoded coefficients.................. 68 viii ABSTRACT Volumetric medical data is a three dimensional view of the object being captured, examples include imaging modalities such as magnetic resonance (MR), computer tomography (CT) and positron emission tomography (PET). For compression of these data, a new 3-D wavelet-based lossless coding scheme, which is an improved version of 2-D PROGRES, is developed. In this compression method, a significance map, which consists of dynamic range values of coefficients of a slice, is designed for reconstructing individual voxels. Significance maps allow any region of interest to be decoded by retrieving only necessary wavelet coefficients. This reduces the number of decoded voxels and the decoding time very significantly. A key novelty of the proposed compression scheme is that it achieves random access decoding of diagonal cuts through the compressed data. This allows flexible display, which can be useful for diagnostic, while not requiring manipulation of very large uncompressed data sets. 1 CHAPTER 1 INTRODUCTION Random access decoding is an important feature provided in volumetric image coding schemes such as 3-D SPIHT 2 , 3-D SBHP 7 , and PROGRES 3 , which all allow regions of interest in a slice to be decoded without decompression of whole volume data. In this research, a diagonal cut through the whole volume data set, representing a 2-D plane, can be visualized by decoding only necessary wavelet coefficients from the compressed volume data. Thus, the new compression algorithm enables rapid random access decoding of regions of interest, including diagonal cuts through the compressed volume data. Volumetric data is used in areas such as medical imaging, fluid dynamics and scientific simulations. Volumetric medical image data is generated by magnetic resonance (MR), computer tomography (CT) and positron emission tomography (PET), which produce a three dimensional view of the object being captured. These three dimensional data are represented as hundreds of two dimensional slices which require very large storage space. Due to the volume data size, using such data sets in uncompressed form can only be possible on computers with extensive amounts of memory. Compression methods reduce the storage space needed for such data. One approach for using volumetric data efficiently is to allow users to access these 2 volumes remotely in compressed form over the network and request individual subvolumes of interest for visualization. Compressing very large volume data at one time is not achievable because the functioning memory space of a typical encoder is not sufficient. 32 bit microprocessors hold 4 GB (2 32 ) of physical memory space. Overhead can happen if the data sets are larger than the maximum memory space. Therefore, very large volume data requires partitioning data sets into blocks and then encoding each block separately. Another benefit of block-based coding is fast compression and decompression. A requested subvolume can be extracted from the compressed data, transmitted in compressed form over the network and then decoded on the user side. One of the most important requirements in developing a compression algorithm is that it must allow quick random access to an individual voxel of compressed data. By enabling this, diagonal cuts through the volume data can be decoded for visualization without requiring complete decoding of all intersecting slices. None of the previous volumetric data compression algorithms provide efficient decoding of diagonal cuts. Some algorithms, such as 3-D SBHP 7 and PROGRES 3 , achieve random access to regions of interest on a slice by slice basis. Other algorithms such as those proposed by Wu 18 and Menegaz 9 provide random access to independent slices. Achieving random access to diagonal cuts from compressed data is potentially very significant since diagonal cuts can provide important information to help understand the whole volume data. Enabling this may allow radiologists to perform better medical examinations. A diagonal cut is visualized as a 2-D plane. Each row of the 2-D image plane represents 3 all the voxels from the cut of each slice. Thus, the height of the 2-D plane is based on the number of slices which are cut. Volumetric data sets are viewed as a collection of slices. The slices represent cross sections of the subject at different positions along a third axis. It is possible to compress the slices independently. However, neighboring slices convey similar spatial details and correlation exists along the third axis. Applying 2-D wavelet transform to volume data will not take advantage of all the redundancy to improve coding efficiency. Thus, a 3-D separable wavelet transform is required to remove inter-slice redundancy. Medical image applications have very strict quality requirements on image compression algorithms. All digital master files from medical imaging devices should be archived in lossless format. For lossless compression of volumetric images, 3-D integer wavelets are preferable since invertible mapping is required from an integer input to an integer wavelet representation. The lifting scheme model presented by Sweldens 4,15 is applied for construction of integer to integer wavelet transform. Many encoders are not well suited to random access of individual data sets because some of them require decoding the whole data to visualize regions of interest in a slice. For instance, run length coding has many benefits given that most of the wavelet coefficients are zero. However, for access to particular wavelet coefficients at the decoder, a large percentage of bits in the bit stream may have to be decoded. Any entropy coding such as arithmetic or Huffman coders will also be limited in terms of allowing extraction of particular piece information quickly from the bit stream. Also, 4 some encoders such as EZW 13 , SPIHT 12 and PROGRES 3 require decoding extra wavelet coefficients for visualization. For instance, if coefficients are encoded in tiles and the region of interest is smaller than the unit tile, there is going to be redundancy. Ihm and Park 5 considered random access ability. In their approach, volume data was partitioned into 16x16x16 unit blocks, and then the 3-D Haar wavelet transform was applied. Wavelet coefficients below specific block thresholds were set to zero after transformation. The transformed coefficients of each unit block were quantized and coded by using cell tag table and cell information, where each cell represented a 4x4x4 subregion of a unit block. The encoding resulted in two bit streams: a code stream and an index stream which provided random access to individual wavelet coefficients. Rodler 11 considered volumetric image encoding as a video coding utilizing bidirectional temporal prediction of frames. The 2-D Haar wavelet transform is applied for the 2-D slices. After transformation, coefficients below slice specific thresholds were set to zero. For fast decoding, two bit stream were designed, one encoding the significance map, by exploiting the hierarchical structure of blocks in a transformed slice, and the other one encoding wavelet coefficients. This method considerably improved both rate and distortion results over Ihm and Park’s algorithm 5 . The EZW 13 algorithm was used for compression of volumetric medical images. The EZW approach was improved by Said and Pearlman with the SPIHT 12 algorithm. This approach has been extended to a 3-D SPIHT 2 algorithm using an integer wavelet packet transform for volumetric image data. EZW and SPIHT are based on bit-plane 5 coding. While bit-plane coding facilitates rate-control, it causes overhead in computation because encoding only one bit of every coefficient is actually necessary for its significance decision. EZW and SPIHT greatly improved the performance of a wavelet encoder but at a much slower encoding and decoding time. PROGRES 3 (progressive resolution decompression) was introduced to reduce the time complexity of set-partitioning based image coding by not implementing the bit-plane model. Instead, a two stage prediction of the dynamic range of coefficients magnitudes is performed. The dynamic range of coefficients magnitudes in a spatial orientation tree is hierarchically represented. This hierarchical dynamic range coding enables resolution scalability. Zero-tree representation is used to improve the compression ratios. In this research, the PROGRES algorithm is improved to enable fast random access to individual voxels in the compressed data, so that the minimum number of wavelet coefficients needs to be decompressed. Fast reconstruction of an individual voxel is a central theme of this study. A new compression scheme that can be used effectively in visualizing very large volume data is introduced. This compression scheme, which is designed for users who have computers with limited memories, can visualize any particular area of a very large volume data rapidly. Most of the previous developed compression techniques trade off random access ability for higher compression ratios. However, the data storage and communication bandwidth costs per megabyte keep dropping. Therefore, in designing the proposed method, random access ability becomes the most important feature of the new compression scheme. 6 Aiming to improve the productivity of radiologists, we strive to achieve the fastest run-time random access ability including diagonal cuts from compressed volume data. Random access decoding is achieved without decompression of the minimum number of wavelet coefficients. 1.1 Thesis Outline Chapter 2 introduces the underlying wavelet theory used throughout the thesis. The discrete wavelet transform and wavelet filter banks were initially used for designing separable wavelet transform. Perfect reconstruction is required for lossless compression methods so the lifting scheme model is applied as integer to integer wavelet transform. 2-D separable wavelet transform is extended to 3-D format for use on volumetric data sets. Chapter 3 presents the 2-D PROGRES 3 algorithm. PROGRES is a very fast, low complexity algorithm and enables resolution scalability and random access ability. The algorithm is based on set-partitioning and does not use bit-plane coding. The coding algorithm is explained extensively and an example showing the encoding steps of PROGRES is provided. Chapter 4 describes the steps which enable random access to diagonal cuts through the volume data. Only the voxels which are cut are examined in order to determine the necessary wavelet coefficients for reconstructing those voxels. Extracting coefficients in two dimensional wavelet domain is introduced. Then, this method is extended for three dimensional data sets. 7 Chapter 5 introduces the 3-D encoder which is an improved version of 2-D PROGRES. Entropy coding is not used in order to preserve random access ability and a significance map is designed for fast reconstruction of individual voxel values. This is a lossless compression method where the compression ratio achieved is 2:1. Chapter 6 explains the 3-D decoder, where any individual voxel values can be reconstructed very rapidly. Only the necessary wavelet coefficients are retrieved for the inverse wavelet transform. Finally, Chapter 7 concludes this thesis with overall conclusion and discussion of further work. 8 CHAPTER 2 WAVELET THEORY 2.1 Introduction The wavelet transform provides a time-frequency representation of the image. The discrete wavelet transform uses the multi-resolution technique by which different frequencies are analyzed with different resolutions. Wavelets are localized basis functions that have their energy concentrated in time or space and are suited for analysis of signals. For the wavelet analysis, first the signal to be analyzed is multiplied with a wavelet function and then the transform is computed for each segment generated. The width of the wavelet function changes with each spectral component. Therefore, at high frequencies the wavelet transform provides good time resolution and poor frequency resolution, while at low frequencies poor time resolution and good frequency resolution are obtained. 2.2 Continuous Wavelet Transform (CWT) The continuous wavelet transform equation is provided by Equation 2-1 where x(t) is the signal to be analyzed and ψ (t) is the basis function or mother wavelet. The scale parameter, s is defined as |1/frequency| and corresponds to frequency information. Scaling either compresses or expands the signal. Low scales (high frequencies) compress the signal and provide the most information of the signal, while 9 high scales (low frequencies) expand the signal and provide the detailed information of the signal. The translation factor is τ. γ(s,τ) = ∫ x(t)ψ* s,τ (t)dt (2-1) The wavelets are generated from a mother wavelet ψ (t) by scaling and translation: ψ s,τ (t) = (1 / √s) ψ((t- τ) / s) (2-2) The wavelet transform is occurring by constantly shifting a continuously scalable function over a signal and calculating the correlation between signals 14 . However, these scale functions do not form an orthogonal basis 6 and the obtained wavelet coefficients will be highly redundant. Thus, CWT is not practical due to its redundancy. Another problem with CWT is that an infinite number of wavelets are required. To overcome these problems, discrete wavelet transform and filter banks are introduced. Discrete wavelets are scaled and translated in discrete steps. 2.3 Discrete Wavelet Transform For discrete wavelet transform, the scaling function takes the values 2 k and the translation component, τ is chosen as 2 k m where k and m are both integers. In this case, Equation 2-2 can be rewritten as: ψ k,m (t) = 2 −k/2 ψ(2 −k t- m) (2-3) These wavelets are used to decompose the image data into a superposition of dilated and shifted wavelets. Thus, x(t) can be expressed as: x(t) = ∑ α k,m ψ k,m (t) (2-4) 10 α k,m is the wavelet coefficient of x(t) and can be calculated as: α k,m = 2 −k/2 ∫ x(t) ψ(2 −k t- m)dt (2-5) The discrete wavelet transform can be implemented using filter banks. A time- scale representation of the digital image is obtained by using digital filtering techniques. The image to be analyzed is passed through filters with different cutoff frequencies at different scales. Discrete wavelets can be made orthogonal by special choices of the mother wavelets. 2.3.1 Classification of Wavelets Wavelets can be classified into two major classes: orthogonal and biorthogonal. Wavelets with symmetry and asymmetry properties are important. Daubechies proved that except for the Haar wavelet, there are no symmetric real orthogonal compactly supported wavelets. Orthogonal Wavelet Filter Banks The coefficients of orthogonal filters are real numbers. The filters are the same length and are not symmetric. The low pass filter L and high pass filter H are related. The two filters are mirror images of each other. Perfect reconstruction is possible with alternating flip. Also for perfect reconstruction, the forward (analysis) filters are identical to the inverse (synthesis) filters except for a time reversal. Wavelet families that form orthogonal bases are preferred because squared error metrics (as may be used to evaluate encoding efficiency) are the same in the time domain and in the 11 transfer domain. Finite impulse response (FIR) filter can be used for both orthogonal and biorthogonal decompositions, but symmetric filters are only available with biorthogonal representations. Biorthogonal Wavelet Filter Banks The coefficients of the biorthogonal filters are either real numbers or integers. The low pass and high pass filters do not have the same length. The low pass filter is always symmetric while high pass filter could be either symmetric or asymmetric. For perfect reconstruction, the analysis filters can be symmetric with odd length filters or one symmetric and the other asymmetric with even length filters. Also, the two sets of analysis and synthesis filters must be dual. The linear phase biorthogonal filters are the most popular filters for data compression applications. 2.3.2 Wavelet Families There are a number of basis functions that can be used as the mother wavelet for wavelet transformation. The mother wavelet determines the characteristics of the wavelet transform because it produces all wavelet functions used in the transformation through scaling and translation. Therefore, the appropriate mother wavelet should be chosen in order to use the wavelet transform effectively. Haar wavelet is the simplest wavelet and capable of performing perfect reconstruction. The wavelets are chosen based on their shape and their ability to analyze the signal in a particular application. 12 2.3.3 Perfect Reconstruction For perfect reconstruction, the following conditions must be achieved; (L a (-z) × L s (z)) + (H a (-z) × H s (z)) = 0 (2-6) (L a (z) × L s (z)) + (H a (z) × H s (z)) = 2 x z -d (2-7) L a and H a correspond to the low pass and high pass analysis filters and L s and H s are low pass and high pass synthesis filters respectively. The first condition implies that the reconstruction is aliasing free and the second condition implies that the amplitude distortion has constant frequency response. The time reversal of the analysis filters is necessary to balance for the delays in the filters. Without the time reversion, it would be impossible to arrive at a non-delayed, perfectly reconstructed signal. If the conditions are satisfied, then all the aliasing caused by subsampling will be cancelled. Perfect reconstruction by classical wavelet transform can be achieved but it is computationally very expensive. In general, while the input images have integer values, the filtered output does not consist of integers due to floating point arithmetic. Loss occurs while rounding off the coefficients. Thus, quantization is applied. For lossless coding, an invertible mapping from an integer input to an integer wavelet representation must be achieved. The lifting scheme is used for constructing integer to integer wavelet transforms 4 . Biorthogonal wavelets constructed by the lifting scheme provide good results for lossless image compression. The image can be reconstructed losslessly because all the coefficients are integers so they can be stored without rounding-off error. Quantization is not required. Another benefit of using the lifting 13 scheme is its computational complexity. The lifting scheme reduces the computational complexity by a factor of 2 as compared to classic wavelet transform. The computation is faster due to integer arithmetic. Since it is an integer wavelet transform, the inverse wavelet transform can be immediately calculated by reversing the operations of the forward transform. 2.4 The Lifting Scheme Model Consider an image s n with 2 n samples. The output of one level of the lifting scheme becomes the input for the next level since the lifting scheme model is recursive. Each step operates on 2 n-k samples where k = 1,2,3….n-1. The 2-level lifting scheme is shown in Figure 2-1. The lifting scheme consists of three steps: split, predict and update. Figure 2-1: 2 Level Lifting Scheme Diagram 14 2.4.1 Split In this stage, the input is split into two disjoint sets of samples. One set consists of the even indexed samples, even i,j , and the other set consists of the odd indexed samples, odd i,j . Each set includes half of samples of the input image. This splitting into even and odd sets is called the Lazy wavelet transform. 2.4.2 Predict The even and odd sets are combined. If these two sets are highly correlated, then it is possible to predict one of the sets when given the other set. In this case, the even set is used for prediction. The forward transform equation is: odd i+1,j = odd i,j – P(even i,j ) (2-8) The odd set odd i,j is subtracted from the predicted even function, P(even i,j ), to obtain a detail image. The predicted transform is calculated without using any auxiliary memory locations by overwriting the locations that hold odd coefficients with the values of detail image. In the Haar case, the predicted even function equals: P(even i,j ) = even i,j (Haar bases) (2-9) If the predicted function in Equation 2-8 is replaced with the value in Equation 2-9, the result is: d = odd i,j – even i,j (2-10) 15 The detail image is d. In the Haar lifting scheme transform, the predicted function replaces an odd coefficient with the difference between odd and even coefficients (these are consecutive coefficients). 2.4.3 Update Odd set is used for the update. The forward transform equation is: even i+1,j = even i,j – U(odd i,j ) (2-11) The even set even i,j is subtracted from the updated odd function, U(odd i,j ), to obtain a coarser image. The locations that hold even coefficients are overwritten by the coarser image. In the Haar case, the updated odd function equals: U(odd i,j ) = (even i,j + odd i,j ) / 2 (Haar bases) (2-12) If the updated function in Equation 2-11 is replaced with the value in Equation 2-12, the result is: a = even i,j – (even i,j + odd i,j ) / 2 a = (even i,j + odd i,j ) / 2 (2-13) The coarser image is represented by a. In the Haar lifting scheme wavelet transform, the update function replaces an even coefficient with the average of the even and odd coefficients (these coefficients are a pair). 16 2.4.4 Inverse Lifting Scheme Wavelet Transform The odd and even set values are replaced by their average, a, and difference d: a = (even i,j + odd i,j ) / 2 d = odd i,j – even i,j The even and odd coefficients are required to be reconstructed. To do that, these two equations are added and solved for odd and even coefficients: even i,j = a – d/2 and odd i,j = a + d/2 (2-14) The inverse transform can be immediately found by undoing the operations of the forward transform. 2.5 Multi-Resolution Analysis Wavelets can be defined by iterations of filters with rescaling. The resolution of the image is determined by the filtering operations and the scale is determined by downsampling and upsampling operations. Lowpass and highpass filtering of the discrete time-domain signal are applied. This connects the continuous-time resolution to discrete-time filters. Multi-resolution analysis was introduced by Mallat 8 . In Figure 2-2, the image is denoted by the sequence x[n] where n is the number of pixels. In this case, n is supposed to be a power of 2. The low pass filter is denoted by L while the high pass filter is denoted by H. At each level, the low pass filter produces coarse approximations of the image, a i [n], and the high pass filter produces detail information of the image, d i [n]. The subscript i denotes the number of 17 decomposition levels. This is a 1-D wavelet transform. The low pass filtered coefficients contain most of the information of the signal so the resulting output coefficients are further decomposed into a coarser approximation and a detail signal. However, the decomposition level is limited by the amount of data. As mentioned above, data can be represented by 2 m and m is the maximum number of decomposition levels. The DWT of the original image is reconstructed by retrieving all the coefficients, a i [n] and d i [n], starting from the last decomposition level. a i [n] and d i [n] are upsampled at every level, passed through low pass filter L` and high pass filter H` and then they are combined. The result is transferred to the next decomposition level, a i-1 [n]. This process is applied through all the decomposition levels and the last result will give the reconstructed image. Figure 2-2: 1-D Wavelet Transform 18 2.6 2-D Separable Wavelet Transform For image compression, the one dimensional wavelet transform can be extended to two dimensions. By choosing separable filters, a two dimensional subband coding scheme can be implemented. This is performed by applying the same filters along the two dimensions iteratively. For an MxN image, first each row of N pixels is filtered and then downsampled by 2. The results are represented in two M x N/2 subbands, which correspond with the high and low pass filtered coefficients of the rows of the image. Second, wavelet transform is applied to each of these subbands along the columns and then downsampled by 2. The result is in four M/2 x N/2 subbands. Figure 2-3 shows the steps of the 2-D separable transform. Figure 2-3: The filtering steps of 2-D separable wavelet transform For inverse 2-D separable wavelet transform, the process of decomposition can be reversed so it is possible to obtain a reconstructed image given that the wavelet coefficients of each subband are available. First, the wavelet coefficients are upsampled by 2 and then filtered along the columns. Second, the filtered column coefficients in two M x N/2 subbands are upsampled by 2 and filtered along the rows. 19 The result gives the MxN reconstructed image. Figure 2-4 shows the 2-D forward and inverse wavelet transform. Figure 2-4: 2-D Forward and Inverse Wavelet Transform The computations of Haar wavelet coefficients are shown below. The pixel values in a 2x2 sub-block are c 1 , c 2 , c 3 , c 4 . c ll is the average coefficient and c lh , c hl , c hh are detail coefficients and they all are stored in their corresponding subbands. The subscripts denote the filters that were applied. For instance, c lh specifies that the low pass filter is used horizontally and then the high pass filter is applied vertically. For reconstruction of the pixel values in the 2x2 sub-block, corresponding wavelet coefficients from each subband are calculated. 20 Two dimensional Haar decomposition is calculated as: c ll = (c 1 + c 2 + c 3 + c 4 )/2 c lh = (c 1 + c 2 − c 3 − c 4 )/2 c hl = (c 1 − c 2 + c 3 − c 4 )/2 c hh = (c 1 − c 2 − c 3 + c 4 )/2 Reconstruction can be written as: c 1 = (c ll + c lh + c hl + c hh )/2 c 2 = (c ll + c lh − c hl − c hh )/2 c 3 = (c ll − c lh + c hl − c hh )/2 c 4 = (c ll − c lh − c hl + c hh )/2 2.7 3-D Separable Wavelet Transform The 1-D wavelet theory explained in section 2.5 can be extended to three dimensions. An orthonormal wavelet basis is used as a primitive of a 3-D method 10 . Mallat 8 constructed a 2-D orthonormal wavelet using a tensor product of two subspaces of L 2 (R).The same method can be used to construct a 3-D orthonormal wavelet. Similar to Figure 2-4, a filtering operation can perform the 3-D wavelet decomposition for a volume data as shown in Figure 2-5. Figure 2-5 shows volume data decomposed into 8 blocks by a 3-D wavelet transform. The coarse approximation of the volume data is represented by a i j , while the detail information is represented in 21 seven subbands, d i j . The subscript i denotes the number of subbands and j is the decomposition level. By recursively applying this filtering operation to the coarse volume data a i j , a multi-resolution representation of the volume data is obtained. Figure 2-5: 3-D Wavelet Subbands A discrete Haar transform, which is the simplest wavelet basis, is applied. The Haar wavelet is computationally simple because it can be implemented by a few integer additions, subtractions and shift operations. The basis is very effective in fast decomposition and reconstruction though it does not achieve the compression ratios that are possible with other filters such as longer Daubechies filters. The Haar wavelets eight-band filter bank can be expressed as: c lll = (c 1 + c 2 + c 3 + c 4 + c 5 + c 6 + c 7 + c 8 )/8 c llh = (c 1 + c 2 + c 3 + c 4 − c 5 − c 6 − c 7 − c 8 )/8 c lhl = (c 1 + c 2 − c 3 − c 4 + c 5 + c 6 − c 7 − c 8 )/8 c lhh = (c 1 + c 2 − c 3 − c 4 − c 5 − c 6 + c 7 + c 8 )/8 c hll = (c 1 − c 2 + c 3 − c 4 + c 5 − c 6 + c 7 − c 8 )/8 22 c hlh = (c 1 − c 2 + c 3 − c 4 − c 5 + c 6 − c 7 + c 8 )/8 c hhl = (c 1 − c 2 − c 3 + c 4 + c 5 − c 6 − c 7 + c 8 )/8 c hhh = (c 1 − c 2 − c 3 + c 4 − c 5 + c 6 + c 7 − c 8 )/8 The reconstruction process is the reverse of decomposition. The reconstruction formulae of Haar wavelets are: c 1 = (c lll + c llh + c lhl + c lhh + c hll + c hlh + c hhl + c hhh )/8 c 2 = (c lll + c llh + c lhl + c lhh − c hll − c hlh − c hhl − c hhh )/8 c 3 = (c lll + c llh − c lhl − c lhh + c hll + c hlh − c hhl − c hhh )/8 c 4 = (c lll + c llh − c lhl − c lhh − c hll − c hlh + c hhl + c hhh )/8 c 5 = (c lll − c llh + c lhl − c lhh + c hll − c hlh + c hhl − c hhh )/8 c 6 = (c lll − c llh + c lhl − c lhh − c hll + c hlh − c hhl + c hhh )/8 c 7 = (c lll − c llh − c lhl + c lhh + c hll − c hlh − c hhl + c hhh )/8 c 8 = (c lll − c llh − c lhl + c lhh − c hll + c hlh + c hhl − c hhh )/8 The 3-D forward wavelet transform is performed by applying the filters along the three dimensions iteratively. For a volume data with a dimension of MxMxM, first M one dimensional voxels are filtered along the x axis and then downsampled by 2. The result is two MxMxM/2 subbands, one corresponding to the high pass filtered coefficients and the other corresponding to the low pass filtered coefficients of rows of the image. Second, filtering is applied to each of these subbands along the y axis and then downsampled by 2. Now there are four MxM/2xM/2 subbands corresponding to LL, LH, HL and HH. Last, the wavelet transform is computed along the z axis and 8 subbands are obtained as LLL, LLH, LHL, LHH, HLL, HLH, HHL and HHH. The 3- 23 D inverse wavelet transform is achieved by reversing the process of 3-D decomposition. First, the coefficients are upsampled then the filtering is applied. The 3-D forward and inverse wavelet transform is shown in Figure 2-6. Figure 2-6: 3-D Forward and Inverse Wavelet Transform 24 CHAPTER 3 IMAGE CODER: PROGRES 3.1 Introduction In this chapter, the PROGRES 3 algorithm is introduced. PROGRES is a fast version of SPIHT and has the ability to perform random access decoding. The hierarchical dynamic range coding enables fast random access decoding of any region of interest at different resolutions. The rate scalability feature of SPIHT is not found in PROGRES. This is because in SPIHT all wavelet coefficients in an image are processed several times at the encoding stage which limits the ability to support random access decoding and greatly increases the computational complexity. In addition, in medical applications, all master files must be in lossless format which requires that all the coefficients must be encoded. The PROGRES algorithm views wavelet coefficients as a collection of spatial orientation trees. Each tree consists of coefficients from all wavelet subbands corresponding to the same spatial area in an image. The energy in each subband tends to decrease with frequency and the statistics in a local neighborhood are similar. Thus, the spatial orientation tree which is a group of wavelet coefficients having the same frequency orientation is designed. Bit plane coding is not used because it places limit on random access decoding. Each coefficient is represented with its sign and magnitude. The dynamic range of a 25 tree of wavelet coefficients is represented as a dynamic range value which gives the number of bits required to represent the magnitude of every coefficient in the tree. The magnitudes of coefficients are decreasing through the tree as the dynamical range decreases. Thus, PROGRES algorithm is designed based on the decaying energy along successive lower to higher frequency of wavelet subbands. The magnitudes of coefficients will be binary coded and the dynamic range values are coded using lossless differential coding. 3.2 Hierarchical Coefficient Dynamic Ranges The dynamic range number of the magnitude of a set of coefficients is the number of bits required to represent all coefficients in the set. c i,j and r i,j are defined as a wavelet coefficient and a dynamic range, respectively. s i,j denotes the set of all coefficients in an orientation tree. The dynamic range of the set is: r i,j = [log 2 ( cp,qЄcs,j max |c p,q | + 1)] (3-1) For example if the maximum magnitude is 14, then the dynamic range value is 4. If the maximum magnitude is 5, the dynamic range value is 3, and so on. If the maximum magnitude is 0, then the dynamic range value is also zero. For a non-zero coefficient one bit is needed to represent its sign. When a tree is divided into subtrees, each subtree will most likely have a decreasing dynamic range number because the root coefficient of the tree which is called parent tree will likely have the highest magnitude in the tree. Thus, at any position in the tree, when considering a subtree and its parent tree, the dynamic range 26 of each subtree will be represented using as a reference of the dynamic range of its parent tree, denoted by r parent . Figure 3-1: Coefficients of trees The dynamic range value of all subtrees, r children , is: r children = max (m,n)Є I(i,j) (r m,n ) (3-2) Using the dynamic range of the parent as a reference: d base = r parent - r children (3-3) d base is encoded instead of r children for dynamic range encoding because P(d base ≤ r children ) > 0.5 for any wavelet transformed image. This value is almost 1 for lower bit rates. Thus, for encoding d base a lesser number of bits are required as compared to encoding r children directly. Then, in the decoder stage, given d base and r parent , r children can be reconstructed and this value will be used as the dynamic range value for the 27 subtrees. Each tree is encoded independently of other trees and d base is calculated for each subtree. Assuming that a tree has 3 subtrees, then there will be 3 different r children , one r children for every subtree. There is one r parent which is the root coefficient of the tree. Computing only d base for 4 children coefficients is an efficient model. However, if the resolution is increased, then the number of coefficients is 16 in this new subtree. Calculating d base is no longer sufficient because the dynamic range of energy in this subtree will vary. This subtree includes 4 sets where each set has 4 coefficients. Each set can have different energy levels, so the dynamic range for each set can be calculated in two stages. First, r children is calculated by d base and second d local is used for calculating each set. r local is the highest coefficient value of the corresponding set: r local = r parents - d base - d local (3-4) The same process can be applied if the resolution is required to be increased. As mentioned before, the differential dynamic range values, d base are coded by lossless unary coding. With the unary code, an integer number is the number of 1’s equal to that integer followed by a single 0. For instance, the integer value 4 corresponds to 11110 (4 ones followed by 0). 8 is the maximum number of bits when the maximum coefficient value is 255. However, depending on the parameters of exponential, generally the average coding rate of unary coding is around 2 bits/symbol. 28 3.3 Coding Algorithm In the beginning of the algorithm, d base which is the difference between the root coefficient in LL subband and its corresponding 2x2 children coefficients in each of LH, HL, HH subbands, is coded independently. For each resolution level, d base , which is the difference of the two dynamic range numbers between 2 adjacent resolutions, is encoded. If any subset has coefficients, then d local is the difference of the children coefficients and the four coefficients in the subset. Figure 3-3 shows the encoding algorithm in detail. The coefficient scanning order is shown in Figure 3-2. The coefficients in the next higher resolution level will not be processed until all the coefficients in the current resolution level have been processed. In any resolution, LL, LH, HL and HH subbands are processed consecutively. Figure 3-2: Scanning order of coefficients 29 Table 3-1: PROGRES encoder algorithm 1. Find the maximum dynamic range value, r parents and perform unary encoding 2. If r parents = 0, then return (There is no coefficient to code) 3. Initialize a list L from the lowest resolution, LL subband and binary encode the root coefficients in the list L by using r parents bits 4. Find the maximum dynamic range value of the children coefficients of the root coefficient 5. Calculate d base by subtracting r children from r parents (d base = r parents - r children ) 6. Unary encode d base 7. Binary encode the value of coefficients of the set in LL by using r children bits. 8. For each resolution level: a. For each set j in current resolution level i. The dynamic range value, r parents , of the children coefficients in current set j is coded ii. The maximum dynamic range value of children coefficients of the subsets, r children is found iii. d base is calculated as d base = r parents - r children iv. Unary encode d base v. If r children = 0, then return to a (there is no coefficients to code) vi. For each subset i, if subset i has children coefficients, 1. Find the maximum dynamic range value of the children of the current subset i 2. Calculate d local by: d local = r children - r subset 3. Unary encode d local 4. If r local = 0, then jump to vii. (there is no coefficients to code) 5. Binary encode the children coefficients of subset i by using r local bits 6. Move subset i to the list L for the next resolution coding vii. Remove the current set j from the list 30 3.4 An Example of PROGRES ENCODER An example of PROGRES 3 coding for the wavelet coefficients of Figure 3-3 is described step by step in Table 3-1. The parenthesized numbers represent the scanning order of wavelet coefficients. For instance, 16 (1) means that 16 is the first coefficient in scanning, 33 (2) is the second coefficient and -15 (3) is the third one in the scanning order. Figure 3-3: Wavelet Coefficients of a 16x16 block First, the initial dynamic range value r parents of the block (16x16) is 8. This means that the maximum coefficient magnitude can be 2 8 -1 which is 255. The coefficient range with sign is [-255,255]. All the coefficients in this block can be represented by 8 bits, but then compression will not be achieved. However, the dynamic range prediction scheme of PROGRES will further reduce the dynamic range throughout the block. The maximum coefficient is 142 as seen in Figure 3-3, which is 31 located at LL subband at the resolution level 0. Thus, the magnitude 142 and its sign which is positive are coded and the total number of bits is 8 +1 which is 9 bits. At resolution level 1, each of the three coefficients, 16 (1), 33 (2), -15 (3) are coded with 6 bits to accommodate the range maximum for 33 (2). Note that the root coefficient (0) has three children, (1), (2) and (3), although all other coefficients have 4 children coefficients. If the current dynamic range value is 0, it means that all the coefficients are zeros, so coding is not required. Considering the resolution level 3 in Table 3-3, the dynamic range value is 0 so encoding is not proceeded. Also if the parent coefficient is equal to zero, processing the algorithm is stopped because all the coefficients in the tree are also zero. Thus, no coefficients will be coded. Considering the resolution level 4 in Table 3-5, it can be seen that r parents is equal to zero and none of the local coefficients are coded. Thus, the algorithm automatically moves to the next block. When r children is equal to zero, the same process is applied. All the coefficients, including the local coefficients, are zero so none of the zero coefficients are encoded. Considering the resolution level 4 in Table 3-3, it can be seen that r children is zero and therefore the coefficients in the block are not coded. Most of the wavelet coefficients are zero so by implementing this algorithm good compression can be achieved. 32 Resolution Level r parents r children d base d local Dynamic Range Coefficient Values (scan order) 0 8 - - - 8 142(0) 1 8 6 2 - 6 16(1), 33(2), -15(3) 2 6 5 1 0 5 4(4), 9(5), 5(6), -19(7) 2 6 5 1 1 4 8(8), -5(9), 6(10), -1(11) 2 6 5 1 1 3 1(12), 3(13), 6(14), 4(15) 3 5 4 1 1 3 6(16), 0(17), 0(18), 2(19) 3 5 4 1 3 1 1(20), -1(21), 0(22), 0(23) 3 5 4 1 2 2 0(24), 1(25), 3(26), 1(27) 3 5 4 1 0 4 -1(28), -4(29), 0(30), -10(31) 3 4 3 1 1 2 3(32), 0(33), 0(34), 0(35) 3 4 3 1 2 1 1(36), 1(37), 1(38), 0(39) Table 3-2: Encoding of wavelet coefficients by PROGRES-1 33 Resolution Level r parents r children d base d local Dynamic Range Coefficient Values (scan order) 3 4 3 1 2 1 0(40), 1(41), -1(42), 0(43) 3 4 3 1 0 3 1(44), -1(45), 6(46), -4(47) 3 3 3 0 3 0 - 3 3 3 0 0 3 0(52), 1(53), 4(54), 0(55) 3 3 3 0 2 1 1(56), 1(57), 0(58), -1(59) 3 3 3 0 1 2 2(60), -1(61), 0(62), -3(63) 4 3 0 3 0 0 - 4 1 1 0 1 0 - 4 1 1 0 1 0 - 4 1 1 0 1 0 - 4 1 1 0 0 1 0(92), 0(93), 0(94), -1(95) 4 2 2 0 1 1 0(96), 1(97), 0(98), -1(99) Table 3-3: Encoding of wavelet coefficients by PROGRES-2 34 Resolution Level r parents r children d base d local Dynamic Range Coefficient Values (scan order) 4 2 2 0 1 1 0(100), 0(101), -1(102), 0(103) 4 2 2 0 1 1 0(104), 0(105), -1(106), 0(107) 4 2 2 0 0 2 0(108), 0(109), -2(110), 0(111) 4 4 2 2 0 2 0(112), 1(113), 0(114), -2(115) 4 4 2 2 0 2 -1(116), 0(117), 0(118), -2(119) 4 4 2 2 2 0 - 4 4 2 2 2 0 - 4 2 2 0 1 1 1(128), 0(129), 0(130), 0(131) 4 2 2 0 2 0 - 4 2 2 0 0 2 2(136),0(137), 0(138),0(139) 4 2 2 0 1 1 -1(140),0(141), 0(142),-1(143) 4 1 1 0 1 0 - Table 3-4: Encoding of wavelet coefficients by PROGRES-3 35 Resolution Level r parents r children d base d local Dynamic Range Coefficient Values (scan order) 4 1 1 0 1 0 - 4 1 1 0 1 0 - 4 1 1 0 0 1 1(156), 0(157), 0(158), 1(159) 4 1 0 1 1 0 - 4 3 2 1 2 0 - 4 3 2 1 0 2 1(180), 1(181), -2(182), 0(183) 4 3 2 1 1 1 -1(184), 0(185), 0(186), 1(187) 4 3 2 1 1 1 1(184), 0(185), 0(186), -1(187) 4 0 0 0 0 0 - 4 3 0 1 0 0 - 4 1 0 2 0 0 - 4 2 0 2 0 0 - Table 3-5: Encoding of wavelet coefficients by PROGRES-4 In the next chapter, we show how to determine the wavelet coefficients necessary for reconstruction in 2-D and 3-D wavelet domain. The following chapter presents the 3- D encoder which is our proposed improved version of 2-D PROGRES. 36 CHAPTER 4 DETERMINING WAVELET COEFFICIENTS FOR VISUALIZING A 2-D PLANE 4.1 Introduction Medical volumetric image data generated by magnetic resonance (MR), computer tomography (CT) or positron emission tomography (PET) produces a three dimensional view of the object being imaged. A typical representation is as a series of two dimensional slices, which require very large storage capacity. Due to their large volume, efficient compression techniques are necessary to reduce the storage space required. For example, clients may access these volume data, which are held in compressed form over the network and request individual subvolumes of interest for visualization. The requested subvolume, such as a section of the brain, lung or heart should be extracted directly from the compressed data, transmitted in compressed form over the network and then decoded and visualized on the client’s side. For instance, if the requested subvolume’s uncompressed size is 256 Mb (512x512x512x2) it would require a long transmission time and this data set would run very slowly in almost all current workstations. Compression would improve this problem significantly. For the visualization of a slice of the requested volume, only the requested slice would be in uncompressed form. Our proposed system enables the requested slice to be any diagonal cut from the compressed volumetric data, so that 37 only that slice will be decoded. In contrast, previous algorithms, which enable random access decoding, are limited to use only existing slices as their elementary units. From the example given above, the requested subvolume has 512 slices, each compressed independently. Therefore, at the decoder, decoding any selected slice is easier because these slices are already predefined by the encoder. However, enabling decoding of a diagonal cut from a subvolume is a more complicated process. For decoding a diagonal cut, the following steps are required: 1. User selects any cut of interest including potentially diagonal ones from a compressed volumetric data. 2. The coordinates of this diagonal cut are put into the system to calculate the plane of the straight cut. 3. The voxels on this plane are identified. 4. The wavelet coefficients needed to reconstruct those voxels are located. 5. Only these located wavelet coefficients are decoded. 6. The inverse wavelet transform is applied to these coefficients to reconstruct the diagonal cut. 7. For visualization, each row of the reconstructed image represents the voxels of each slice which is being cut. In this chapter, the first 4 steps will be explained. The last 3 steps will be discussed in the next chapter. Enabling user selection of any diagonal cut from a compressed volumetric data can enhance the utility of the data set, as compared to allowing only horizontal slices to be displayed. This method allows users to quickly 38 browse any data of interest through a compressed volume data without requiring decompression of the whole data. 4.2 Region of Interest The first step is to prompt the user to identify a cut within the volumetric data. The coordinates representing the cut are entered in the system. Now, the voxels on this cut have to be identified because only those voxels are required to be visualized. The goal is to avoid errors in identifying the voxels, since identifying more voxels than needed would require additional rate and complexity. Each straight cut generates a plane. Using the plane equation, the voxels on the plane can be found. Equation of a plane One of the methods for determining a plane is by specifying three points. The three points are expressed as: p 1 =           1 1 1 z y x p 2 =           2 2 2 z y x p 3 =           3 3 3 z y x Normal of the plane is calculated by: n = | ) ( ) ( | ) ( ) ( 1 2 1 3 1 2 1 3 x x x x x x x x − × − − × − and n =           c b a (4-1) If n is the normal to the plane, all points p on the plane satisfy the following: p • n = d (4-2) 39 When given any 3 points, we derive from the equation above: 0 ) ( ) ( ) ( 1 1 1 = − + − + − c z z b y y a x x 0 ) ( 1 1 1 = + + − + + cz by ax cz by ax ) ( 1 1 1 cz by ax d + + − = 0 = + + + d cz by ax (4-3) Equation 4-3 is the plane equation and all its parameters are calculated using any 3 points from the plane. This equation is used to find the voxels on the plane. 4.3 Discovering the Voxels on the Plane of Region of Interest 4.3.1 Algorithm Description So far, the three points on the plane and the plane equation are known and now all the voxels on this plane of region of interest will be found. In this algorithm, from the three points, the first two points are the entrance points to the volume data and the last point is the exit point out of the volume data. The slices which are being cut first and last are known. By using the plane equation, the voxels which are being cut on each slice can be determined. Figure 4-1 shows a volumetric data which is cut by a plane. The blue plane is the diagonal cut. p 1 and p 2 are two points on the slice which is cut first. p 3 is on the slice which is cut last. Table 4-1 describes how to identify the voxels on the plane 40 Figure 4-1: A 2-D diagonal plane cuts the volumetric data Table 4-1: Steps of determining the voxels on a given 2-D plane from a volume data 1. Three points are entered by a user (two points are on the same first cut slice and the third point on the same last slice). 2. Plane equation is calculated (described above in detail). 3. The plane equation is applied on each slice independently. a. The slope of the line on the slice is calculated. b. If the magnitude of the slope is less than or equal to 1 i. Set the image width value, x to 0 (x=0) ii. Add 0.001 to x (x=x+0.001) iii. Use the line equation to obtain the image height value, y iv. If value y is greater or equal to zero and less than the image height, by (y>=0 & y<by), then the integer values of x and y will represent the voxel’s location and this result will be stored. v. Now add 0.999 to x (x=x+0.999) vi. Use the line equation to obtain the image height, y vii. Check if integer values for both x and y are obtained. If this voxel is previously located, increment x by 1 and go to step ii. viii. If y is greater or equal to zero and less than the image height, the new voxel is located and the result is stored. Now, increment x by 1 and go to step ii. ix. Repeat this process until x is less than the image width, bx c. If the magnitude of the slope is greater than 1 i. Follow the same steps with y and x exchanged. 41 The slope of the line from each slice has to be calculated because when the magnitude of the slope is less than 1, the growth rate of x values is greater than that of y values. Therefore, when locating the voxels on a line, x values must be used to obtain y values by using the line equation. However, using y values does not give the correct voxels location on a line. In addition, using integer x values is not enough to identify all voxels on a line, because between two integer values more than one voxel can be covered by a line. Therefore, the best method is to check right before the border and after the border of a pixel block. This way, none of the voxels will be missed. If the magnitude of the slope is larger than 1, the opposite equation is applied; y values will be used to obtain x values on a line. Again, y values right before and after the border are checked. This process is continued until the shifted value is less than its axis size. Examples of determining voxels on a given plane are given in the Appendix. After all the voxels on the plane are discovered in the volume domain, the wavelet coefficients corresponding to these voxels must be decoded and the inverse wavelet transform can be applied to visualize the plane. How the wavelet coefficients are decoded will be explained in Chapter 6. In this chapter, we will discuss which wavelet coefficients are required for reconstructing the plane. 42 4.4 Determining the Wavelet Coefficients The 3-D volumetric data set is partitioned into blocks with a dimension of 16×16×16. For each block, a three dimensional wavelet transform is applied 4 times. The eight-band filter bank, which is based on Haar wavelets, includes 8 coefficients in a 2x2x2 subregion of a block. One coefficient is the average of all the voxel values in the subregion. The seven remaining coefficients are the detail coefficients which provide the detail information of the subregion. The eight coefficients on the first decomposition level are used to reconstruct the eight voxel values of the corresponding subregion. After the reconstruction process of the first level, the first reconstructed coefficient represents the average of all the voxel values in the subregion of the second level. Therefore, 7 other wavelet coefficients are required for the reconstruction of the second level. Again, 7 wavelet coefficients are needed for the third and fourth decomposition levels. Total, 29 wavelet coefficients must be decoded from the encoded unit block. Seven additions/subtractions are performed per formula and evaluation of a total of 28 additions/subtractions is necessary. Higher decomposition levels will require more wavelet coefficients. For level 3, 22 wavelet coefficients and 21 additions/subtractions are required. For level 2, 15 wavelet coefficients and 14 additions/subtractions are needed. The number of decomposition levels affects the decoding time directly. For example, higher levels increase the decoding time and lower levels reduce it. However, by using higher levels, more wavelet coefficients are zero, which increases the compression ratio 43 significantly. As a trade off between coding efficiency and complexity, a 4 level decomposition is chosen. 4.4.1 Extracting Coefficients in Two Dimensional Wavelet Domain First, given a random line on a 2-D image, the necessary wavelet coefficients will be determined for reconstructing the line. After this is discussed in detail, the method can be extended for reconstructing a plane from a volumetric data. In Figure 4-2, a 16x16 image is shown. This image data is cut by a diagonal cut which is shown in blue color. The colored pixels represent the pixels cut by the line. There are five different colors and each one is represented by a 2x2 block since 2-D Haar filter size is 2x2 (4 coefficients) in two dimensional wavelet domain. Again in this block, one coefficient is the average of all coefficients and the others are the detail coefficients. Considering each block individually, it can be seen that not all of the block coefficients are cut by the line. However, when 2-D inverse transform is applied, 4 wavelet coefficients are used to reconstruct 4 pixel values. As a matter of fact, even if only one pixel is needed, because of transformation, 4 pixels will be obtained. Now the wavelet coefficients should be determined for reconstruction. Table 4- 3 shows all necessary wavelet coefficients need to decode each of the 2x2 blocks. For every block, 13 wavelet coefficients are determined, 4 from the first level and 3 from the others. Table 4-3a shows all the wavelet coefficients needed to reconstruct the 44 pixels in the green block. For reconstruction, the 2-D inverse transform must be applied 4 times. The transformation starts from the last level, where magnitudes of the wavelet coefficients are largest. In Table 4-3a, 4 wavelet coefficients are transformed to obtain the next level’s wavelet coefficients. The obtained wavelet coefficients are located in the LL band which is used for the transformation of second level. After the second level transformation, again the resulting coefficients are located in the LL band. The third transformation is applied by using one of the resulting wavelet coefficients in the LL band. Last wavelet transformation is applied to reconstruct the pixels which are colored as green in Figure 4-2. The same idea is implemented to reconstruct other pixel blocks. Table 4-2 shows the steps for locating the coefficients of a 4-level 2-D wavelet transform. It will start from first level and continue to locate wavelet coefficients until the 4 th level. However, the decoder starts by decoding the 4 th level wavelet coefficients. Considering the image blocks in Table 4-3, it will be observed that some of the necessary wavelet coefficients are the same for reconstruction. In that case, only different wavelet coefficients will be stored in the system. This way, the same wavelet coefficients will not be retrieved multiple times thus avoiding redundancy. Table 4-3f shows the combined wavelet coefficients, which are all the wavelet coefficients that have to be retrieved for reconstruction of the colored pixels needed to visualize the line in Figure 4-2. 45 Figure 4-2: Blocks of pixels which are reconstructed for visualizing the line Table 4-2: Steps of locating coefficients in two dimensional wavelet domain 1. Given coordinates of a coefficient value (x,y) of 4 level 2-D wavelet transform 2. Divide the coordinates by 2 and obtain the result in integer format. This gives the coefficient in the LL 1 band, int( , ) x y 2 2 . 3. Add 8 to x and get the coefficient in LH 1 band, int ( , ) x y 2 8 2 + . 4. Add 8 to y and get the coefficient in HL 1 band, int ( , ) x y 2 2 8 + . 5. Add 8 to x and y and get the coefficient in HH 1 band, int( , ) x y 2 8 2 8 + + 6. Divide the coefficient in the LL 1 band by 2 and obtain the result in integer format which gives the coefficient in the LL 2 band, int( , ) x y 4 4 7. Add 4 to x and get the coefficient in LH 1 band, int ( , ) x y 4 4 4 + . 8. Add 4 to y and get the coefficient in HL 1 band, int ( , ) x y 4 4 4 + . 9. Add 4 to x and y and get the coefficient in HH 1 band, int( , ) x y 4 4 4 4 + + . 10. Follow the same process for obtaining the coefficients in level 3 and 4. 46 4-2a 4-2b 4-2c 4-2d 4-2e 4-2f Table 4-3: From 4-2a through 4-2e, wavelet coefficients for reconstructing each block of pixels are shown. 4-2f shows combined wavelet coefficients. 47 4.4.2 Extracting Coefficients in Three Dimensional Wavelet Domain Now the idea can be extended to a 3-D data to allow the visualization of a selected random plane from a volumetric image. As explained before, in three dimensional wavelet domain, the Haar filter has 8 coefficients (2x2x2). In Figure 4-3, multiple slices are combined together to generate a volumetric data. The decomposition level of the 2-D wavelet transform is 3. In volumetric data, the decomposition level of the 3-D wavelet transform is also 3. 2x2x2 filter is shifted throughout the volume, although 2x2 filter is shifted in the slice. Volumetric data with dimension 8x8x4 is shown in Table 4-4a. The volume is cut by a plane shown in blue. All the voxels that are needed to visualize the plane are shown in color. Every 2x2x2 block is shown in different colors representing each Haar filter. Figure 4-3: Generating a 3-D volume data from combining multiple slices 48 For the blocks on the front side (the first two slices), the wavelet coefficients needed for reconstruction are shown in Table 4-4b through 4-4e. The back side blocks are not illustrated since they are treated similarly and displaying the coefficients is not very easy. For determining the wavelet coefficients, the same process used in 2-D is followed. The decomposition level of the 3-D wavelet transform is 3 and for the first level 8 coefficients are used for transformation, with the transformed coefficients located in LLL band. One of the coefficients obtained in LLL band is used for the next level transformation. And again the new coefficients are generated in LLL band and the last wavelet transform is applied including an obtained coefficient to reconstruct the voxels in a 2x2x2 block. Extracting coefficients from a three dimensional data set is the same as in the two dimensional cases with only difference being the Haar filter size. The next chapter introduces the new 3-D encoder which is an extended version of 2-D PROGRES. In addition, a significance map is developed for enabling to recover individual wavelet coefficients so that only necessary wavelet coefficients have to be retrieved for visualizing a selected 2-D plane from volume data. 49 Table 4-4: The wavelet coefficients needed for reconstruction are shown in Table 4-5b through 4-5e. 4.5f shows all necessary coefficients for reconstruction of a 2-D plane. 4-5c 4-5a 4-5b 4-5d 4-5e 4-5f 50 CHAPTER 5 VOLUMETRIC IMAGE CODER 5.1 Introduction Muraki 10 introduced the idea of using wavelets to efficiently approximate 3-D data, as was explained in Chapter 2. Volumetric data sets can be viewed as a collection of slices. The slices represent cross sections of the subject at different positions along a third axis. It is possible to compress the slices independently. However, neighboring slices convey similar spatial details so that there is correlation along the third axis. Since applying a 2-D wavelet transform to a volume data will not remove all the redundancy, a 3-D separable wavelet transform can be used to remove inter-slice redundancy. Medical image applications impose very strict quality requirements on image compression algorithms. All digital master files from medical imaging devices should be archived in lossless form. For lossless compression of volumetric images, 3-D integer wavelets must be used since an invertible mapping is required from an integer input to an integer wavelet representation. The lifting scheme is used for constructing integer to integer wavelet transform as discussed in Chapter 2. Visualizing very large volumetric data has always been a problem for radiologists. This problem can be addressed with very high compression ratios but then the negative diagnostic and legal implications of lossy compression become an 51 issue. Lossless coding must be applied but its compression ratio varies 2:1 to 4:1 which is very small compared to lossy compression. A possible solution to this problem is to design a system that allows radiologists to browse rapidly through a volume data and then select a particular area for closer examination. This requires random access to a compressed volume and decoding only a particular area without decompression of the whole data. Many encoders are not well suited to random access of individual data items since some of them require decoding the whole data to visualize a region of interest. For instance, run-length coding has many benefits given the fact that most of the wavelet coefficients are usually zero. However, for accessing particular bits, all the previous bits in the bit stream must be decoded. Also, some encoders decode additional wavelet coefficients for visualization. For instance, if the bits are tile-based coded then, it is very possible that the decoder will decode based on the tile size. But if the region of interest is smaller than the unit tile, there is going to be redundancy. The designed system decodes only the necessary wavelet coefficients which are needed for inverse wavelet transform. Therefore, the designed system promises a very fast random access decoding, with reduced memory space for visualization of region of interest. The new 3-D encoder is introduced in this chapter. The uncoded 2-D PROGRES algorithm is extended for 3-D encoding. PROGRES algorithm is selected because it is very suitable to access independent wavelet coefficients. The new encoder does not include additional entropy coding since it makes the random access 52 very difficult and increases the computation complexity involved in enabling random access. 5.2 3-D PROGRES 3-D PROGRES is a straightforward extension of 2-D PROGRES 3 . The PROGRES algorithm views wavelet coefficients as a collection of spatial orientation trees. Each tree consists of coefficients from all wavelet subbands. In 2-D wavelet domain, the number of subbands is 4, while for 3-D wavelets, there are 8 subbands. One subband is the average of all coefficients and other 7 subbands include detail information of volume data. So the 3-D spatial orientation tree represents the parent- child relationship in seven orientations, three correspond to intra-slice information and four to inter-slice. 2-D and 3-D wavelet transforms are explained in Chapter 2. Again in 3-D PROGRES, the magnitudes of wavelet coefficients will be binary coded and the dynamic range values are unary coded. d base , r parent , r children , d local , r local are calculated in the same way as in 2-D PROGRES. The only difference will be the number of coefficients in the subtrees. Recall, from Chapter 3, that each tree is encoded independently of other trees and it includes r parent which is the root coefficient of the tree. Considering the subtrees, r children , which is the highest coefficient value of the corresponding subtree is computed. d base is calculated by taking the difference between r parent and r children . d base is computed based on 4 children coefficients. However, in 3-D PROGRES, d base is 53 calculated based on 8 coefficients due to 3-D Haar filter size. And if the resolution is increased, the number of coefficients will be 16 in the new subtree of 2-D data. However, in 3-D PROGRES, there will be 32 coefficients in that subtree. Calculating d base is not efficient due to variations in the dynamic range. Thus, this subtree is divided into 4 sets where each set has 8 coefficients. r local is the highest coefficient value of the corresponding set. Because of different energy level, the dynamic range for each set, d local can be calculated in two steps. 3-D spatial orientation tree is shown in Figure 5-1. Subbands and filter size are the main difference between 2-D and 3-D PROGRES. From the explanation given above, it can be said that a 3-D encoder is computationally more expensive than a 2-D encoder. However, a 3-D encoder can exploit inter-slice redundancy for better coding efficiency, as compared to a 2-D encoder. Also, more zero wavelet coefficients are obtained when a 3-D encoder is applied to the volume data instead of a 2-D encoder. Compressing a very large volume data directly may not be possible since the working memory space of typical encoders may be limited. 32 bit microprocessors support 2 32 (4 GB) physical memory space. And if the volume data is larger than the maximum memory space, virtual memory occurs and that causes overhead in encoder. Therefore, very large volume data requires dividing data into blocks and then encoding each block independently. Coding of an individual unit volume allows rapid compression and decompression. 54 Figure 5-1: 3-D Spatial Orientation Tree 5.3 Block-based Coding A volume data has a dimension V x xV y xV z and it is divided into blocks of dimensions 16x16x16. The original dimensions, V x ,V y ,V z are preferred if they are multiples of 16. Otherwise, the volume data requires zero padding. The linear unit size of 16 is proper for the purpose of random access of portions of the entire data that are needed for visualization. If the dimension is greater than 16, then decoding may produce a significant amount of volume data which may not be needed for a particular query. Thus, decoding time and storage space will be increased. However, compression efficiency will be improved since the block takes advantage of the correlation among more volumetric data samples. On the other hand, if the dimension is less than 16, the blocks get so small that encoding will not be very efficient. But decoding the particular subvolumes will be faster due to their small size, and will require less memory space. 3-D compression algorithm is applied to each unit block separately. 55 5.3.1 Encoding The encoder takes the wavelet coefficients coming from the quantizer or wavelet transform and attempts to represent the bit stream as efficiently as possible without any loss. The dynamic range values are unary coded. The magnitudes of wavelet coefficients are binary coded based on the selected dynamic range value. Coefficients are scanned following the order shown in Figure 3-3 and are coded using the steps described in Table 3-1. The wavelet coefficients of a 16x16x16 block are shown at Figure 5-2. The first coefficient, which is also expressed as r parents is 142 and its dynamic range is 8. First the dynamic range is unary coded which takes 9 bits (111111110). Later, first the sign then magnitude of the coefficient are binary encoded (110001110). If the coefficient is positive, the sign bit represents as 1 and if it is negative, then it is 0. Please refer Chapter 3 for the details of the PROGRES algorithm. The 3-D PROGRES is tested on several volume data and the lossless compression ratio is achieved as from 1.9:1 to 2.1:1. For extracting a bit stream of a region of interest, the positions of the wavelet coefficients must be known at the decoder so that it is possible to determine which wavelet coefficients need to be decoded. A significance map is designed to complement the 3-D PROGRES algorithm to determine the positions of the wavelet coefficients. With this system any individual coefficients can be located. 56 Figure 5-2: Wavelet coefficients of a first slice of a 16x16x16 block 5.4 Designing a Significance Map In 3-D PROGRES, d base and d local are the dynamic range values of coefficients and sets of 8 coefficients are encoded based on these dynamic range values. Designing a significance map based on the dynamic range values is sufficient to determine the positions of coefficients. The dynamic range is the number of bits to represent the coefficient. Thus, the magnitude of the coefficient does not have to be known by the decoder to locate the coefficient. Coding the dynamic range values in the right scanning order is sufficient. At the decoder, first, the entire significance map will be decoded then the necessary wavelet coefficients will be extracted from the coefficient bit stream. The decoder part will be explained in the next chapter. 57 The significance map derived from the wavelet coefficients in Figure 5-2 is shown in Figure 5-3. For understanding the design of significance map, let us assume that the first slice in Figure 5-2 includes all the highest coefficient values of their corresponding subtrees. The reason is that it is hard to show the coefficients of rear slices in a cube and the highest coefficient values play a significant role in determining dynamic range values. Figure 5-3: A significance map The significance map represents the coefficient set of 16x16x2 due to Haar filter (2x2x2) so 8 significance maps are needed for a 16x16x16 block. The significance map does not include only the first r parents value since it represents the fourth decomposition level of LLL band and it is the only coefficient in the band like the other bands in fourth level. So the dynamic range value of r parents is encoded separately. The first value in the map, 6, is the dynamic range value of other bands in fourth level including 7 coefficients. All other dynamic range values in the map are used for coding 8 wavelet coefficients. There are two different bit streams, one 58 corresponding to wavelet coefficients and the other one corresponding to values in significance map. Only 4 bits are needed for representing each value since the highest value can be 8. Thus, the significance map size is 32 bytes, with 4 bits for encoding r parents and 12 bits for offset. The offset value is the address of the first coefficient of a slice. 2 12 is larger than the number of coefficients in a given volume data. Thus, these additional 2 bytes are added to the significance map size, which gives a total bit budget of 34 bytes for representing the significance map of a slice. For the entire unit block (16x16x16), 272 bytes are needed for the significance maps. The unit block size is 4096 bytes (16x16x16x1) so designing a significance map for determining the positions of coefficients increases the bit budget by %6.6. This is the overhead of the system. The voxel is represented by one byte but in some cases, it may take up 2 bytes, then the percentage of overhead will be halved. As a matter of fact, having 272 bytes of overhead for a 16x16x16 unit block does not affect the compression ratio very significantly. Random accessing to significance maps is also required because otherwise, all the significance maps of the volume data must be decoded first for determining the individual wavelet coefficients. That would be a very time consuming operation. Therefore, the bits in the map bit stream will be binary encoded. At the decoder, it is known that representing the offset value requires 12 bits and all other values have 4 bits each. Based on their dynamic range values, any value can be extracted rapidly. 59 The next chapter presents the 3-D decoder where values of significance maps and wavelet coefficients are decoded. Fast random access decoding can be achieved. In addition, 3-D inverse wavelet transform with a decomposition level of 4 is applied to each 8 wavelet coefficients separately. 60 CHAPTER 6 3-D RANDOM ACCESS DECODER 6.1 Introduction The proposed system allows accessing to individual data without decoding the whole compressed data. Both encoder and decoder of the system play significant roles for enabling random access. At the 3-D encoder, first, the 3-D wavelet transform is applied to obtain the wavelet coefficients, then, significance maps for locating coefficients are binary encoded and coefficients are encoded by 3-D PROGRES. The 3-D encoder is explained extensively in Chapter 5. The 3-D decoder algorithm is applied to each block (16x16x16) independently. The process of extracting a voxel value from compressed volumetric data consists of four stages. First, all the wavelet coefficients necessary for reconstruction are determined. This part was explained in Chapter 4. From a given plane, which constitutes the region of interest, all necessary unit blocks and wavelet coefficients can be determined from the whole volumetric data for visualization. Second, the values of significance maps of the necessary coefficients are decoded. The significance map assists positioning the wavelet coefficients. Third, by using the significance maps, the necessary coefficients are located in the coefficient bit stream and retrieved. And last, the reconstruction formula is applied to the wavelet coefficients to obtain voxel values. 61 6.1. Decoding significance maps Significance maps are very effective for positioning the wavelet coefficients. Significance maps are binary encoded so that the decoder can use them for random access. The decoder will receive following bit size information from the header. Each value takes 4 bits so one significance map is represented by 32 bytes since there are 64 values in one map. In addition to that, 4 bits for r parents and 12 bits for offset are used in encoding each significance map. Therefore the total cost of encoding significance map is 34 bytes. One significance map corresponds to 16x16x2 so one unit block needs 8 significance maps. From a given plane, all necessary wavelet coefficients can be determined. By determining the wavelet coefficients, the significance maps are determined, too. Since only the wavelet coefficients’ slice and block information are sufficient for determining which significance map represents these wavelet coefficients. As an example, assume a necessary wavelet coefficient is located at position (x,y,z) = (9,0,6) in block 18. The direction z gives the slice information so this wavelet coefficient is positioned at the third slice. This means that the third significance map of this block needs to be decoded. It is highly possible that multiple coefficients are located in the same significance map. To avoid decoding the same map multiple times, first, all necessary significance maps are determined from all given wavelet coefficients. Each block takes 272 bytes and each significance map needs 34 bytes so the binary decoding will be started after (272x17) + (34x2) = 4692 bytes. First, the offset 62 and r parents are decoded, then all the values in the significance map are retrieved. The decoding process will be stopped after retrieving 34 bytes from the map bit stream. From these 34 bytes, the wavelet coefficient at (9, 0, 6) can be extracted. 6.2 Decoding wavelet coefficients Once all necessary wavelet coefficients are identified and the necessary significance maps for extracting wavelet coefficients are decoded, individual wavelet coefficients can be retrieved. Let us use the coefficient given in the example above. The coefficient is located at (9,4,6). Its dynamic range value, which is shown in blue in Figure 6-1, is located at (4,2,3) due to 3-D Haar filter (2x2x2). For locating the coefficient, the values up to this dynamic range value must be considered. First the offset value, which contains the address of the first wavelet coefficient of a slice then the dynamic range value of r parents and its magnitude coefficient are decoded. If r parents is 8 then the number of bits for retrieving the parent coefficient will be 8+1+1+8. The first 1 is representing the zero of the dynamic range value from unary coding and the second 1 is representing the sign bit. Recall that the dynamic range value is the number of bits to represent the magnitude of a coefficient. Now the map can be decoded based on scanning order of encoding coefficients. First value, 6, corresponds to 7 and all others correspond to 8 wavelet coefficients. Thus, the position of the wavelet coefficient at (9, 4, 6) is calculated as: p = (12+4) + (8+1) + (8+1) + (8-6+1) + ((6+1)x7) + (6-5+1) + (5-5+1) + ((5+1)x8) + 63 (5-4+1) + ((4+1)x8) + (5-3+1) + ((3+1)x8) + …………………...+ (3-0+1) + (3-3+1). The values in Figure 6-1 represent the dynamic range of values for each 2x2x2 subblock. For encoding these values, the difference between the parents and children values and also the difference between children and local values are encoded to save bits. At last, based on the local dynamic ranges, the wavelet coefficients are binary encoded. Refer Chapter 5 for more information on the 3-D encoder. Figure 6-1: A significance map representing 16x16x2 coefficients of a unit block and blue colored dynamic range value includes the seeking coefficient at (9, 4, 6). The calculation of the position of the local dynamic range shown above is not the final position of the coefficient. Since 8 coefficients are encoded with that local dynamic range value, then we need to access the specific coefficient we are seeking. Its location is (9,4,6) and it is represented as (4.5,2,3) in a 2x2x2 subblock. 4.5 represents the second coefficient on the x axis, 2 represents the first coefficient on the y axis and 3 represents the first slice on the z axis in the 2x2x2 subblock. From this result, coefficient to be looked-up is the second encoded coefficient in the subblock. Therefore, the number of bits of 5 coefficients are added to the equation above and it 64 is expressed as p + ((3+1)x5). This is the final result of the position of the coefficient at (9,4,6). This is illustrated in Figure 6-2. Figure 6-2: The 2x2x2 subblock of blue colored dynamic range value in Figure 6-1. The blue colored coefficient represents the seeking coefficient at (9, 4, 6). 6.3 3-D Inverse Haar Wavelet Transform After all necessary wavelet coefficients are retrieved, the 3-D inverse Haar wavelet transform is applied 4 times for decomposition to visualize the given plane. 8 wavelet coefficients one from each subband are needed for 3-D Haar wavelet transform. After the transformation, a 2x2x2 subblock containing 8 voxel values are reconstructed. Please refer Chapter 2 for 3-D Haar transform. 6.4 Experiments For experimental tests, 256x256x128 MRI head scan volume data set is used. Each voxel is represented with one byte, and the whole data set takes up 10 Mbytes. A lossless compression algorithm has been implemented on these volume data sets. The whole volume data is divided into 16x16x16 unit blocks and each unit block is 65 encoded separately. The lossless compression rates for this volumetric image sets are 3.72 bit/ voxel. In the Menegaz paper 9 , the performance of different compression schemes such as 3-D EZW, JPEG2000 16 and JPEG-LS 17 is compared on MRI head scan data with a dimension of 256x256x128. Table 6-1 summarizes the performance of these algorithms including Menegaz proposed algorithm, 2-D PROG 9 and our proposed scheme 3-D PROP which achieves the best lossless compression efficiency. Table 6-1: Lossless Compression Performances (Bit/Voxel) The new algorithm enables random access to wavelet coefficients in the compressed volume data that are needed to visualize an arbitrary plane. All previous algorithms provide random accessing only to individual slices, not allowing decoding of diagonal cuts. These cuts are visualized as a 2-D plane. In this plane, each row of the 2-D image plane represents the voxels from the cut of each slice. For instance, the volume data, which has dimensions 256x256x164, are completely cut with a 45° angle where 0° angle is the x axis and angles increase counterclockwise. One voxel from each row and 256 voxels from the whole slice are cut so that each row of the plane representing the cut has 256 voxels. There are 164 rows since 164 slices of the volume data are cut. Thus, the plane size is 256x164. This plane is shown in Figure 6-3b. Other cuts are at angles 0°, 90°, 135° shown also in Figure 6-3. It can be seen that these cuts have significant information from the volume data which can be very 66 important for medical examination. Thus, enabling random access to diagonal cuts from a compressed volume data is a great opportunity for obtaining new information about the data sets. The compression ratios of these diagonal cuts of 2-D planes are very similar and their results are shown in Figure 6-4. 6-3a 6-3b 6-3c 6-3d Figure 6-3: Diagonal cuts through the volume data with an angle of 0°, 45°, 90°, 135° shown from 6-3a to 6-3d. 67 3.5 3.6 3.7 3.8 3.9 4 4.1 4.2 0 45 90 135 Angles of 2-D planes Lossless rate(bit/voxel) Figure 6-4: Lossless bit rate (bit/voxel) for 2-D diagonal planes with different angles The bit-rate varies significantly depending on the plane location in the volume data set. Since some planes can be represented by so many zero wavelet coefficients which increases the compression ratio significantly. If the values in a significance map are zero, which means that all the coefficients in the corresponding subblock are zero, then it is not necessary to decode any wavelet coefficients. This method has two very important benefits. First, the number of bits per voxel will be very small and, second, fast random access is achieved since decoding is not taken a place. In addition, the size of the plane also affects the compression ratio significantly. Small planes with very few zero coefficients have very small compression ratios. 2-D diagonal planes with 45° and 135° from 256x164 to 64x164 are randomly selected from volumetric data and shown from Table 6-2 to Table 6-4. Overall, the bit-rate varies from 1.4 to 4.2 bit/voxel. The number of bits per voxel and the number of decoded wavelet coefficients are shown in Figure 6-5. 68 0 1 2 3 4 5 6 0 4096 8192 12288 16384 20480 24576 28672 32768 36864 Number of decoded voxels Lossless rate (bit/voxel) Figure 6-5: Lossless bit rate (bit/voxel) and number of decoded coefficients Resolution Angle 45° Angle 135° 256x164 224x164 Table 6-2: Random access decoding of cuts with 45° and 135° from resolution 256x164 to 224x164. 69 Resolution Angle 45° Angle 135° 192x164 160x164 Table 6-3: Random access decoding of diagonal cuts with 45° and 135° from resolution 192x164 to 160x164. 70 Resolution Angle 45° Angle 135° 128x164 92x164 Table 6-4: Random access decoding of diagonal cuts with 45° and 135° from resolution 128x164 to 92x164. 71 Resolution Angle 45° Angle 135° 64x164 Table 6-5: Random access decoding of diagonal cuts with 45° and 135° with resolution 64x164. 72 CHAPTER 7 CONCLUSION 7.1 Thesis Summary Many algorithms for compressing volumetric data have been proposed in the last decade. However, most of the research focus has been on efficient transmission and storage of volumes. However, this only solves one part of the problem, which is associated with handling such very large volume data. The problem rises when the volumes have to be used in uncompressed format for further medical examinations. Decoding the whole data sets is only possible on computers with extensive amounts of memory. Recently though, the need for keeping the volumes compressed while they are accessed by other applications has been recognized, thus the challenge of supporting random access to the compressed data is introduced. Fast reconstruction of an individual voxel in a 3-D wavelet domain is the central theme of this thesis. Enabling random access to an individual voxel from compressed data can improve the productivity of radiologists. The compression scheme is designed for extracting any regions of interest through the volume data rapidly. Thus, user can browse through a large volume data set very quickly and select any particular area for closer examination. The analysis of wavelet transform is performed. Perfect reconstruction is very significant for medical applications since all digital master files from medical imaging 73 devices should be archived in lossless format. Lifting scheme is described for lossless compression of volumetric images. 3-D separable Haar wavelet transform is applied since it is the simplest and computationally cheapest transform. It can be implemented by a few integer additions, subtractions and shift operations. This transform is very effective in fast decomposition and reconstruction and reduces inter-slice redundancy significantly. 2-D PROGRES algorithm is introduced. The algorithm is based on set- partitioning and non bit-plane coding is used, since bit plane coding causes overhead in computation. Instead of bit-plane coding, each wavelet coefficient, which is represented by sign and magnitude, is processed only once. However, in SPIHT each coefficient needs multiple processes to be completely reconstructed. Therefore, PROGRES decodes faster than SPIHT. The image in wavelet domain is divided into a number of non-overlapping spatial orientation trees, each one rooted by a coefficient in the lowest frequency subband. The dynamic range of coefficients is expressed as a dynamic range number which is the number of bits required to represent every magnitude of the coefficient in the tree. The dynamic ranges of descendant subtrees are most likely smaller than the parent tree. Thus, this decrease in dynamic range number is coded using lossless differential coding. Enabling decompression of a diagonal cut throughout the volume data without decoding any additional wavelet coefficients is explained. This diagonal cut is displayed as a 2-D plane for visualization where first, the voxels which are cut are 74 determined. Then, wavelet coefficients necessary for inverse wavelet transform are located. For 3-D separable wavelet transform, 8 wavelet coefficients from 8 different subbands are used for reconstructing 8 voxel values in a 2x2x2 subblock. 2-D PROGRES is improved for a 3-D encoder which supports fast reconstruction of an individual voxel. For retrieving individual voxels, a significance map representing the dynamic range values of wavelet coefficients must be designed for positioning the necessary wavelet coefficients. The significance maps are binary encoded and no additional coding is applied due to fast random access ability. Offset which gives the address of the first coefficient in a slice is encoded with each significance map. The overhead for encoding significance maps of a unit block 16x16x16 is 272 bytes. For 8bpp, the overhead is 6% of the bit budget. In the 3-D decoder, first the significance maps, corresponding to the selected voxels are decoded. In the significance map, the offset value and the necessary dynamic range values are used for determining the wavelet coefficients. After computing positions of the wavelet coefficients, these coefficients are retrieved from the coefficient bit stream. Then 3-D inverse Haar wavelet transform is applied to each 8 wavelet coefficients separately. After the inverse transform, only the necessary voxels values are reconstructed for visualizing any regions of interests including diagonal cuts from compressed volume data. 75 7.2 Future Work Fast reconstruction of an individual voxel is the central theme of this thesis. For higher coding efficiency, arithmetic coding can be considered. Significance maps are designed for decompression of individual wavelet coefficients but causes overhead in compression. The overhead can be reduced by applying additional entropy coding. The significance map consists of dynamic range values of coefficients which mean that unary encoding of dynamic range values in the coefficient bit stream is redundancy and can be removed. Dynamic ranges are coded by lossless differential coding in the coefficient bit stream. Thus, the coefficient stream will only include binary encoded wavelet coefficients which would increase the compression efficiency and decrease the computational complexity of random access decoding significantly. 76 REFERENCES 1. C.-L. Chang, A. Maleki and B. Girod, “Adaptive Wavelet Transform for Image Compression via Directional Quincunx Lifting,” Proc. IEEE International Workshop on Multimedia Signal Processing, MMSP-05, Shanghai, China, November 2005. 2. E. Christophe and W. A. Pearlman, “Three-dimensional SPIHT Coding of Volume Images with Random Access and Resolution Scalability,” IEEE Transactions on Image Processing, Dec. 13, 2006. 3. Y. Cho, A. Said and W. A. Pearlman, “Low complexity resolution progressive image coding algorithm: PROGRES (Progressive Resolution decompression),” in 2005 IEEE International Conference on Image Processing (ICIP '05), Sep 2005. 4. I. Daubechies and W.Sweldens, “Factoring wavelet transforms into lifting steps,” J. Fourier Anal. Appl., vol. 4, no. 3, pp. 245-267, 1998. 5. I. Ihm and S. Park, “Wavelet-based 3D compression scheme of interactive visualization of very large volume data,” Computer Graphics Forum, pp. 3-15, 1999. 6. G. Kaiser, “A friendly guide to wavelets,” Boston, Birkhäuser, 1994. 7. Y. Liu and W. A. Pearlman, “Scalable Three-Dimensional SBHP Algorithm with Region of Interest Access and Low Complexity,” Applications of Digital Image Processing XXIX , Proc. SPIE Vol. 6312, pp. 631209-1-11, Aug. 2006. 8. S. G. Mallat, “A theory for multiresolution signal decomposition: The wavelet representation,” IEEE Transactions on Patern Analysis and Machine Intelligence, vol.11, pp. 674-693, July 1989. 9. G. Menegaz and J.P. Thiran, “Three-dimensional encoding/two-dimensional decoding of medical data,” Medical Imaging, IEEE Transactions on Volume 22, Issue 3, pp. 424 – 440, March 2003. 10. S. Muraki, “Volume data and wavelet transforms,” IEEE Computer Graphics & Applications, 13(4), pp.50-56, July 1993. 11. F.F. Rodler, “Wavelet based 3D compression with fast random access for very large volume data,” Computer Graphics and Applications, 1999, Proceedings, 7. Pacific Conference, pp. 108 – 117, Oct. 1999. 77 12. A. Said and W. A. Pearlman, “A new fast and efficient image codec based on set partitioning in hierarchical trees,” IEEE Transactions on Circuits and Systems for Video Technology, vol. 6, pp. 243–250, June 1996. 13. J. M. Shapiro, “Embedded image coding using zerotrees of wavelet coefficients,” IEEE Transactions on Signal Processing, vol. 41, pp. 3445–3462, 1993. 14. Y. Sheng, “The transforms and applications handbook,” Ed. by A. D. Poularikas, P. 747-827. Boca Raton, Fl (USA): CRC Press, The Electrical Engineering Handbook Series, 1996. 15. W. Sweldens, “The lifting scheme: A custom-design construction of biorthogonal wavelets,” Appl. Comput. Harmon. Anal., vol. 3, no. 2, pp. 186-200, 1996. 16. B. E. Usewitch, “A tutorial on modern lossy wavelet image compression: Foundations of JPEG 2000, “IEEE Signal Processing Magazine, pp. 22-35, September 20001. 17. M. Weinberger, G. Seroussi, G. Sapiro, "The LOCO-I Lossless Image Compression Algorithm: Principles and Standardization into JPEG-LS", IEEE Trans. Image Processing, Vol. 9, pp.1309-1324, August 2000. 18. X. Wu and T. Qiu, “Wavelet coding of volumetric medical images for high throughput and operability,” Medical Imaging, IEEE Transactions on Volume 24, Issue 6, pp. 719 – 727, June 2005. 78 APPENDIX An Example of Determining Voxels on a Plane The image size is kept small (4x4) in order to make the example short and easy to follow. The three points are as follows: p 1 = (0,2.75,0), p 2 = (3,4,0),p 3 = (2,0,4) n r = | ) ( ) ( | ) ( ) ( 1 2 1 3 1 2 1 3 x x x x x x x x − × − − × − = (a, b, c) = (5,−12,−10.75) d = p ּ n r = (0 × 5) + (2.75 × −12) + (0 × −10.75) = −33 ax + by + cz = d → 5x −12y −10.75z = −33 m = |a / b| = |5 / 12| = 0.41 z=0 5x − 12y − 10.75z = −33 x = (12y/5) − (33/5) y = 0 → x = −6.6 y = 1 → x = −4.2 y = 2 → x = −1.8 y = 3 → x = 0.6 Given three points of a plane, the plane equation is determined. Then, the magnitude of the slope is calculated as 0.41, which is less than 1. This means that x values will be shifted throughout the image to calculate y values of a line. This is applied to every slice. The opposite equation, which uses y values to calculate x values, has incorrect results. For example, when y is equal to 2, x is supposed to be 0. However, a negative value is calculated. 79 z=0 5x − 12y − 10.75z = −33 y = (5x/12) + (33/12) x = 0 → y = 2.75 x = 1 → y = 3.16 x = 2 → y = 3.58 x = 3 → y = 4 z=1 5x − 12y −10.75z = −33 y = (5x/12) + (22.25/12) x = 0 → y = 1.85 x = 1 → y = 2.27 x = 2 → y = 2.68 x = 3 → y = 3.1 z=2 5x − 12y − 10.75z = −33 y = (5x/12) + (11.5/12) x = 0 → y = 0.95 x = 1 → y = 1.37 x = 2 → y = 1.79 x = 3 → y = 2.20 z=3 5x − 12y − 10.75z = −33 y = (5x/12) + (0.75/12) x = 0 → y = 0.06 x = 1 → y = 0.48 x = 2 → y = 0.89 x = 3 → y = 1.31 Figure 1: Determining the voxels of a line with a slope of less than 1 Table 1 shows all the voxel indices, which are located on the plane. The script z is the slice number and x represents the horizontal axis values. For instance, (0,2,0) is read as (x,y,z), which is the index value of a voxel. Considering Figure 1, for each slice, only integer x values are being used to calculate y values of a line. However, some voxels can be missed. The bolded voxel indexes in Table 1 are not being located. Using the method explained in Table 4-1, all the voxels will be discovered. Implementation of this method in this example is shown in Figure 2. 80 z=0 z=1 z=2 z=3 x=0 (0,2,0), (0,3,0) (0,1,1), (0,2,1) (0,0,2),(0,1,2) (0,0,3) x=1 (1,3,0) (1,2,1) (1,1,2) (1,0,3) x=2 (2,3,0) (2,2,1),(2,3,1) (2,1,2),(2,2,2) (2,0,3),(2,1,3) x=3 − (3,3,1) (3,2,2) (3,1,3) Table 1: The voxels indexes located on a line with a slope of smaller than 1 z=0 5x − 12y − 10.75z = −33 y = (5x/12) + (33/12) x=0.001→y = 2.75 x=0.999→y = 3.16 x=1.001→y = 3.16 x=1.999→y = 3.58 x=2.001→y = 3.58 x=2.999→y = 3.99 x=3.001→y = 4 x=3.999→y = 4.41 z=1 5x − 12y −10.75z = −33 y = (5x/12) + (22.25/12) x=0.001→y = 1.85 x=0.999→y = 2.26 x=1.001→y = 2.27 x=1.999→y = 2.68 x=2.001→y = 2.68 x=2.999→y = 3.09 x=3.001→y = 3.1 x=3.999→y = 3.51 z=2 5x − 12y − 10.75z = −33 y = (5x/12) + (11.5/12) x=0.001→y = 0.95 x=0.999→y = 1.36 x=1.001→y = 1.37 x=1.999→y = 1.78 x=2.001→y = 1.79 x=2.999→y = 2.20 x=3.001→y = 2.20 x=3.999→y = 2.61 z=3 5x − 12y − 10.75z = −33 y = (5x/12) + (0.75/12) x=0.001→y = 0.06 x=0.999→y = 0.47 x=1.001→y = 0.48 x=1.999→y = 0.89 x=2.001→y = 0.89 x=2.999→y = 1.31 x=3.001→y = 1.31 x=3.999→y = 1.72 Figure 2: Determining the voxels of a line with a slope of less than 1. Each side of the border is checked for catching different voxels on the line. 81 Each side of the border is checked to see whether any y values have changed. The bolded numbers in Figure 2 are the voxels that are missed by the previous method shown in Figure 1. It requires twice as much computation as the previous method, but it perfectly discovers all the voxels on a plane. If the index of the voxel was previously stored in the system, that index will not be stored a second time and will go to next step. There are no additional voxels which are not located on the plane, thus avoiding decompression of wavelet coefficients to reconstruct these voxels. Below is another example with a magnitude of slope higher than 1. The three points on a plane, the plane equation and the slope are shown as: p 1 = (3, 0, 0), p 2 = (4, 3, 0), p 3 = (1, 4, 4) → 6x − 2y + 5z = 18 m = |6/-2| = 3 Due to the magnitude being larger than 1, y values are shifted through the image to calculate x values of the line. This method for the calculation is applied and shown in detail in Figure 3. Table 2 shows all the indexes of voxels which are located on the plane. z=0 z=1 z=2 z=3 y=0 (3,0,0) (2,0,1) (1,0,2) (0,0,3) y=1 (3,1,0) (2,1,1) (1,1,2) (0,1,3), (1,1,3) y=2 (3,2,0) (2,2,1), (3,2,1) (2,2,2) (1,2,3) y=3 − (3,3,1) (2,3,2) (1,3,3) Table 2: The voxels indexes located on a line with a slope of larger than 1 82 z=0 6x − 2y + 5z = 18 x = (y/3) + 3 y=0.001→x = 3 y=0.999→x = 3.33 y=1.001→x = 3.33 y=1.999→x = 3.66 y=2.001→x = 3.66 y=2.999→x = 3.99 y=3.001→x = 4 y=3.999→x = 4.33 z=1 6x − 2y + 5z = 18 x = (y/3) + (13/6) y=0.001→x = 2.16 y=0.999→x = 2.49 y=1.001→x = 2.49 y=1.999→x = 2.82 y=2.001→x = 2.82 y=2.999→x = 3.15 y=3.001→x = 3.16 y=3.999→x = 3.49 z=2 6x − 2y + 5z = 18 x = (y/3) + (4/3) y=0.001→x = 1.33 y=0.999→x = 1.66 y=1.001→x = 1.66 y=1.999→x = 1.99 y=2.001→x = 2 y=2.999→x = 2.33 y=3.001→x = 2.33 y=3.999→x = 2.66 z=3 6x − 2y + 5z = 18 x = (y/3) + (1/2) y=0.001→x = 0.5 y=0.999→x = 0.83 y=1.001→x = 0.83 y=1.999→x = 1.16 y=2.001→x = 1.16 y=2.999→x = 1.49 y=3.001→x = 1.5 y=3.999→x = 1.83 Figure 3: Determining the voxels of a line with a slope of larger than 1

Abstract (if available)

Abstract Volumetric medical data is a three dimensional view of the object being captured, examples include imaging modalities such as magnetic resonance (MR), computer tomography (CT) and positron emission tomography (PET). For compression of these data, a new 3-D wavelet-based lossless coding scheme, which is an improved version of 2-D PROGRES, is developed. In this compression method, a significance map, which consists of dynamic range values of coefficients of a slice, is designed for reconstructing individual voxels. Significance maps allow any region of interest to be decoded by retrieving only necessary wavelet coefficients. This reduces the number of decoded voxels and the decoding time very significantly. A key novelty of the proposed compression scheme is that it achieves random access decoding of diagonal cuts through the compressed data. This allows flexible display, which can be useful for diagnostic, while not requiring manipulation of very large uncompressed data sets.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Parallel implementations of the discrete wavelet transform and hyperspectral data compression on reconfigurable platforms: approach, methodology and practical considerations

PDF

Complexity scalable and robust motion estimation for video compression

PDF

Distributed wavelet compression algorithms for wireless sensor networks

PDF

Lifting transforms on graphs: theory and applications

PDF

Distributed source coding for image and video applications

PDF

Joint routing and compression in sensor networks: from theory to practice

PDF

Compression algorithms for distributed classification with applications to distributed speech recognition

PDF

Mocap data compression: algorithms and performance evaluation

PDF

Robust video transmission in erasure networks with network coding

PDF

Advanced techniques for green image coding via hierarchical vector quantization

PDF

Efficient coding techniques for high definition video

PDF

Predictive coding tools in multi-view video compression

PDF

Statistical lesion detection in dynamic positron emission tomography

PDF

WOLAP: wavelet-based on-line analytical processing

PDF

Compression of signal on graphs with the application to image and video coding

PDF

Transmission tomography for high contrast media based on sparse data

PDF

H.264/AVC decoder complexity modeling and its applications

PDF

Error tolerant multimedia compression system

PDF

Power efficient multimedia applications on embedded systems

PDF

Critically sampled wavelet filterbanks on graphs

Asset Metadata

Creator Ayberk Ozturk (author)

Core Title Random access to compressed volumetric data

School Viterbi School of Engineering

Degree Master of Science

Degree Program Electrical Engineering

Publication Date 07/24/2007

Defense Date 05/29/2007

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag lossless compression,OAI-PMH Harvest,random access decoding,volumetric medical data,wavelet transform

Language English

Advisor Antonio Ortega (committee chair), Kuo, C.-C. Jay (committee member), Richard M. Leahy (committee member)

Creator Email ayberkoz@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-m653

Unique identifier UC168182

Identifier etd-Ozturk-20070724 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-520237 (legacy record id),usctheses-m653 (legacy record id)

Legacy Identifier etd-Ozturk-20070724.pdf

Dmrecord 520237

Document Type Thesis

Rights Ayberk Ozturk

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Repository Name Libraries, University of Southern California

Repository Location Los Angeles, California

Repository Email cisadmin@lib.usc.edu