![Fig. 10. Extension of the spatial support in the previous ( 1) and the next ( + 1) subband images. Squares with same pattern represent voxels whose significance states are combined in the definition of the corresponding [k; l; j].](/figures/fig-10-extension-of-the-spatial-support-in-the-previous-1-2gocvgr3.webp)
424 IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 22, NO. 3, MARCH 2003
Three-Dimensional Encoding/Two-Dimensional
Decoding of Medical Data
Gloria Menegaz*, Member, IEEE, and Jean-Philippe Thiran, Member, IEEE
Abstract—We propose a fully three-dimensional (3-D)
wavelet-based coding system featuring 3-D encoding/two-dimen-
sional (2-D) decoding functionalities. A fully 3-D transform is
combined with context adaptive arithmetic coding; 2-D decoding
is enabled by encoding every 2-D subband image independently.
The system allows a finely graded up to lossless quality scalability
on any 2-D image of the dataset. Fast access to 2-D images is
obtained by decoding only the corresponding information thus
avoiding the reconstruction of the entire volume. The performance
has been evaluated on a set of volumetric data and compared to
that provided by other 3-D as well as 2-D coding systems. Results
show a substantial improvement in coding efficiency (up to 33%)
on volumes featuring good correlation properties along the
axis.
Even though we did not address the complexity issue, we expect a
decoding time of the order of one second/image after optimization.
In summary, the proposed 3-D/2-D multidimensional layered zero
coding system provides the improvement in compression efficiency
attainable with 3-D systems without sacrificing the effectiveness
in accessing the single images characteristic of 2-D ones.
Index Terms—3-D/2-D, compression, lossless, volumetric data,
wavelets.
I. INTRODUCTION
M
OST of the current medical imaging techniques produce
three-dimensional (3-D) data distributions. Some of
them are intrinsically volumetric, like magnetic resonance
(MR), computerized tomography (CT), positron emission to-
mography (PET), and 3-D ultrasound, while others describe the
temporal evolution of a dynamic phenomenon as a sequence of
two-dimensional (2-D) images, so that they are more properly
labeled as 2-D+time. The huge amount of data generated every
day in the clinical environment has triggered considerable
research in the field of volumetric data compression for their
efficient storage and transmission. The basic idea is to take
advantage of the correlation among the data samples in the 3-D
space to improve compression efficiency. The most widespread
approach combines a 3-D decorrelating transform with the
extension of a coding algorithm that has proved to be effective
on 2-D images. In [1], the 3-D version of the set partitioning
in hierarchical trees (SPIHT) [2] algorithm for image com-
pression is applied to volumetric medical images. The same
guideline is followed in [3], where the authors also address the
Manuscript received December 10, 2000; revised August 30, 2002. Asterisk
indicates corresponding author.
*G. Menegaz is with the Department of Computer Science, University of
Fribourg, CH-1700 Fribourg, Switzerland (e-mail: gloria@ieee.org).
J.- P. Thiran is with the Signal Processing Institute, School of Engineering
Techniques and Sciences, Swiss Federal Institute of Technology, CH-1015
Lausanne, Switzerland.
Digital Object Identifier 10.1109/TMI.2003.809689
problem of context modeling for efficient entropy coding. The
performance of the 3-D extension of the embedded zerotree
wavelet (EZW)-based coding algorithm [4] is analyzed in
[5]–[7]. A slightly different approach is described in [8], where
a 3-D-DCT is followed by quantization, adaptive bit allocation
and Huffman encoding. In [9]–[11], a 3-D separable wavelet
transform is used to remove interslice redundancy, while in
[12] different sets of wavelet filters are used in in the (
)
plane and
direction, respectively, to account for the difference
between the intraslice and and interslice resolution.
This led to the common consensus that the exploitation of the
full 3-D data correlation potentially improves compression. The
main drawback of 3-D systems is computational complexity. If
an increase in the encoding time might be tolerated, a swift de-
coding is of prime importance for the efficient access to the data.
A possible solution has been proposed in [1] and [5]. It consists
in splitting the volume in coding units of 8 or 16 images each
and processing those independently in order to save memory
and reduce the coding time. Coding units are fixed a priori,as
well as the number of images which are decoded at one time.
Our solutionthatisbasedonthe observationthat it is common
practice to analyze 3-D data distributions one image at a time
for medical examination. Accordingly, in order to be suitable
within a picture archiving and communication system (PACS)
a coding system must provide a fast access to the single 2-D
images. In the proposed solution, the decoding time is kept low
by minimizing the amount of information to be decoded to re-
construct any 2-D image (or, more in general, subset of images)
of the dataset. This is accomplished by independently encoding
each subband image, and making the corresponding information
accessible through the introduction of some special characters
(i.e., markers) into the bitstream. Once the user has specified
the position of the image of interest along the
axis, the set
of subband images that are needed for its reconstruction is de-
termined and the related information is decoded. The inverse
discrete wavelet transform (IDWT) is performed locally and the
single image is recovered. The coding scheme is based on the
multirate 3-D subband coding of video described in [13]. What
we retain is the strategy used for entropy coding, namely the
multidimensional context-adaptive arithmetic coding [14]. The
subtended subband structure is nevertheless different. We per-
form a 3-D-DWT on the volume instead of treating differently
the spatial and temporal dimensions.
The paper is organized as follows. Section II gives an
overview on the global system. In Section III, the lifting
scheme and integer wavelet transform are revisited. Section IV
describes the procedure followed to determine the set of sub-
band images needed to reconstruct a given image of interest.
0278-0062/03$17.00 © 2003 IEEE
MENEGAZ AND THIRAN: THREE-DIMENSIONAL ENCODING/TW-DIMENSIONAL DECODING OF MEDICAL DATA 425
Fig. 1. Volumetric data. We call
z
the third dimension, and assume that the
images are the intersections of the volume with a plan orthogonal to
z
axis.
The coding principle is presented in Section V and Section VI
illustrates the different working modalities. The compression
performance is analyzed in Section VII, and Section VIII
derives conclusions.
II. T
HE 3-D/2-D MULTIDIMENSIONAL LAYERED
ZERO CODING (MLZC) SYSTEM
The combination of the 3-D wavelettransform with an ad-hoc
coding strategy provides high coding efficiency and fast access
to any 2-D image of the dataset. Given the index of the image
of interest along the
axis ( coordinate in Fig. 1), the corre-
sponding portion of the bitstream is accessed and decoded to
recover it at the desired quality. At the encoder, the data are first
decorrelated by a 3-D DWT and then encoded via the MLZC
technique. At the decoder, the set of wavelet coefficients nec-
essary to reconstruct an image of index
is automatically de-
termined and only the corresponding parts of the bitstream are
decoded. The IDWT is performed locally, reducing the memory
requirements and the computational cost.
The wavelet transform has many features that make it suitable
for our application. The approximation properties of reasonably
smooth signals have determined the success of wavelet-based
techniques for image compression. Noteworthy, the JPEG2000
standard [15] follows the same approach. The implementation
via the lifting steps scheme [16] is particularly advantageous
in this framework. First, it provides a very simple way of con-
structing nonlinear wavelet transforms mapping integer-to-in-
tegervalues[17].This is veryimportant for medical applications
because it enables lossless coding. Second, perfect reconstruc-
tion is guaranteed by construction for any kind of signal exten-
sion along borders. This greatly simplifies the management of
the boundary conditions and facilitates the selection of the coef-
ficientsneeded to reconstruct an image. Third, it is computation-
ally efficient. It can be shown that the lifting steps implementa-
tion asymptotically reduces the computational complexity by a
factor 4 with respect to the classical filter-bank implementation
[18]. Finally, the transformation can be implemented in-place,
namely progressively updating the values of the original sam-
ples, without allocating auxiliary memory.
The 3-D-DWT is followed by successive approximation
quantization and context adaptive arithmetic coding. Some
markers are placed in the codestream for the random access to
the encoded information. By combining the 3-D-DWT with 2-D
spatial neighborhoods for entropy coding, the resulting MLZC
algorithm features 3-D encoding/2-D decoding capabilities.
However, many degrees of freedom are left for the design of
the system. The shape of the spatial support of the neighborhood
defining the context and the placement rule of the markers in the
bitstream lead to different working modes. The global-progres-
sive (G-PROG) mode is obtained by encoding the volume as a
whole and without putting any marker. This mode provides the
best compression efficiency. Both 2-D and 3-D contexts can be
used. The resulting bitstream is fully embedded, supporting a
finely-graded range of bit-rates ensuring scalable quality on the
volume, but 2-D decoding is not possible. The layer-per-layer
(LPL), and LPL progressive (LPL-PROG) modes are obtained
by adding some markers in order to enable random access to
the information of interest in the bitstream. More specifically,
the LPL mode provides random access to every subband image.
The idea is to decode the entire information concerning the set
of subband images needed to reconstruct the image of interest
at full quality (i.e., lossless). To achieve quality scalability on
the final 2-D image, other markers must be added, leading to
the LPL-PROG mode. Direct access is possible to every quan-
tization layer of every subband image. Scalable quality is ob-
tained by successively decoding the quantization layers, i.e.,
the bitplanes, of the concerned subband images. The drawback
is the bitstream over-heading due to the additional informa-
tion needed for data addressing, which reduces compression ef-
ficiency. Fig. 2 summarizes the three working modalities and
illustrates the position of the markers in the bitstream. In the
figure, H is for the Header of the bitstream, and
represents
the quantization layer
of the subband image at position in a
given 3-D subband. In the G-PROG mode, the whole informa-
tion concerning the quantization layer
is encoded for all the
subband images. In the LPL mode, all quantization layers are
located in the same segment and markers are placed only be-
tween
and , being the number of quantization steps.
This reduces the number of markers while preserving 2-D de-
coding capabilities, improving the compression efficiency at the
expense of the SNR scalability in lossy regime. Indeed, such a
mode is intended for recovering the 2-D image of interest at full
quality. Finally, in the LPL-PROG mode, the order is the same
but markers are put between
and , , .
III. I
NTEGER WAVELET TRANSFORM VIA LIFTING
The spatial correlation among data samples is exploited by
a fully 3-D separable wavelet transform. The signal is suc-
cessively filtered and down-sampled in all spatial dimensions.
The decomposition is iterated on the approximation low-pass
band, which contains most of the energy [19]. Fig. 3 shows the
classical filter-bank implementation of the DWT. The forward
transform uses two analysis filters,
(low-pass) and
(bandpass), followed by subsampling, while the inverse trans-
form first up-samples and then applies two synthesis filters,
426 IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 22, NO. 3, MARCH 2003
Fig. 2. MLZC working modalities. H: bitstream header;
L
quantization layer
i
of the subband image at position
j
in a given 3-D subband. In the G-PROG mode,
the whole information concerning the quantization layer
i
is encoded for all the subband images. In LPL mode, all the quantization layers are located in the same
segment and markers are only between
L
and
L
,
n
being the number of quantization steps. In the LPL-PROG mode, the order is the same but markers are
put between
L
and
L
,
8
i
,
j
.
Fig. 3. DWT.
(low-pass) and (bandpass). Fig. 4 shows a two levels
DWT on a natural image. The approximation subband is a
coarser version of the original, while the other subbands rep-
resent the high frequencies (details) in the horizontal, vertical
and diagonal direction, respectively.
In the proposed system, the DWT is implemented according
to the recently developed lifting steps scheme [16]. The lifting
scheme provides a way to perform any DWT with finite filters
with a finite number of lifting steps. The lifting steps represen-
tation of a given filter is obtained by the Euclidean factoriza-
tion of the polyphase matrix (see Fig. 5) of the filter bank into
a sequence of 2
2 upper and lower triangular matrices. The
polyphase matrix
is defined as
(1)
where
(2)
and , respectively, and , are the even and
odd polyphase components of the synthesis filter
, respec-
tively,
. If the determinant of is equal to one, then the
filter pair (
)iscomplementary. In this case, the following
theorem holds [16]:
Theorem 1: Given a complementary filter pair (
), then
there always exist Laurent polynomials
and for
and a nonzero constant so that
(3)
The block diagrams for the forward and inverse transforms are
illustrated in Figs. 6 and 7, respectively. Each triangular matrix
corresponds to one lifting step. The number of lifting steps
depends on both the length of the filters and the factorization.
It is worth noticing that the result of the Euclidean factoriza-
tion is not unique, so many lifting representations are possible
for the same
. From Figs. 6 and 7 it is easy to realize that
the synthesis chain can be obtained by mirroring the filter-bank
from the analysis counterpart and changing the sign of the fil-
ters. The global system can be seen as a sequence of do/undo
steps, for which the perfect reconstruction property is ensured
by construction. This provides additional degrees of freedom in
the design of the filters, allowing any nonlinear operations into
thebasic blocks and any kindof signal extension outside the bor-
ders. In particular, the integer DWT is obtained by introducing
a rounding operation after each lifting step [17]. As mentioned
in Section II, the availability of an integer version of the trans-
form enables lossless coding and makes the algorithm suitable
for the implementation on a device. However, the integer coef-
ficients are approximations of those that would be obtained by
projecting the signal on the original wavelet basis. This can be
modeled by an equivalent noise which becomes noticeablewhen
the hypothesis of high-resolution quantization holds. It can be
shownthat it introduces an additional contribution to the quanti-
MENEGAZ AND THIRAN: THREE-DIMENSIONAL ENCODING/TW-DIMENSIONAL DECODING OF MEDICAL DATA 427
Fig. 4. DWT of a natural image. (a) Original; (b) DWT subbands for a two level decomposition. The approximation subband is a coarser version of the original,
while the other subbands represent the high frequencies (details) in the horizontal, vertical and diagonal direction.
Fig. 5. Polyphase representation of the wavelet transform.
zation noise, which degrades the rate/distortion performance of
the coding system [20]. Furthermore, it is responsible for an os-
cillatory trend of the PSNR along the
axis, making the quality
of the reconstructed image dependent on its position within the
volume. The analysis of such a phenomenon is out of the scope
of this paper. We refer to [21] for more details. What it is impor-
tant to mention here is that the amount of such noise is propor-
tional to the number of rounding operations, which in turn de-
pends on the decomposition depth and the lifting chain length.
Accordingly,we haverestrictedthe choice of filters to the family
of the interpolating filters [22] admitting a two-steps chain
(4)
As the choice of the filter-bank is not critical for compression
performances, we choose the 5
3 [22] filter. Being extremely
short – two steps of length two each – it minimizes the number
of subband images to decode for recovering the 2-D image of
interest, as will be discussed in Section IV.
The particular structure of the lifting chain facilitates the
determination of the set of subband images needed for the
point-wise IDWT (PW-IDWT). The separability of the trans-
form allows to map such a task to the one-dimensional (1-D)
case. The core of the problem consists in finding the set of
subband coefficients needed to recover one signal sample.
Then, results can be easily extended to intervals (i.e., signal
segments) and, eventually, multiple dimensions.
IV. P
OINT-WISE IDWT
In this section, we formalize the PW-IDWT. It is basically a
1-D problem: each pixel of the image to recover is regarded as
the sample
in position of the 1-D signal observed along
the parallel to the
axis passing trough it. Correspondingly, the
set of subband coefficients that are needed for its reconstruction
by IDWT maps to the coordinates of the subband images along
the
axis.
The proposed solution exploits the inherent recursive nature
of the wavelet transform. The IDWT is an iterative process
starting at the coarsest scale: the approximation subband at
the finer (
) level is reconstructed by filtering the set of
coefficients at the coarser
level according to [19]
(5)
428 IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 22, NO. 3, MARCH 2003
Fig. 6. Forward wavelet transform using lifting.
Fig. 7. Inverse wavelet transform using lifting.
where and are the approximation and detail subbands,
respectively, and
is the decomposition level, which
increases with the depth of the decomposition. The signal
is
reconstructed by iterating such a procedure for
. The
number of coefficients taking part to the convolution in a given
subband depends on the length of the filter and on the number
of decomposition levels. The method used to determine the po-
sitions of the involved coefficients in each subband consists in
climbing back the synthesis filter-bank and keeping track of the
positions of the subband coefficients that get involved step by
step. Given the position
of the sample of interest in the signal
domain, we start by identifying the set of coefficients
that are needed at the finest resolution (i.e., ). Here, is
the subband index and takes the values
for approximation and
for details, respectively. For doing this, we look into the syn-
thesis chain from its output, and follow it step by step, keeping
track of the samples needed by the lifting steps filters. Due to
the recursiveness of the IDWT, given
the procedure is
iterated to get
, at the next finer resolution (i.e., ).
The only difference is that now there is a set of samples to be
recovered {
} instead of a single one. The iteration of
such a procedure for
results in the complete
set of necessary subband coefficients.
The procedure can be easily generalized to sets {
} of samples in the signal space. Let identify the
coefficients in subband
at level needed to reconstruct the
signal sample in position
. Then, the solution for the set of
samples {
}is
(6)
Formula (6) also applies to subband intervals. It is worth men-
tioning here that
depends on being even or odd. In
general, with the usual structure of the lifting scheme starting
with an
-type step, odd indexed samples correspond to larger
. We refer to the Appendix A for the details.
TABLE I
N
UMBER OF SAMPLES
GP(
l
;j
)
IN SUBBAND (
l; j
)NEEDED TO RECOVER
THE
SAMPLE AT POSITION
k
USING THE 5
2
3FILTER AND FOR
L
=3
LEVELS OF DECOMPOSITION.
TABLE II
N
UMBER OF SAMPLES
GP(
l;j
)
IN SUBBAND (
l; j
)NEEDED TO RECOVER
THE
SAMPLE AT POSITION
k
USING THE 9/7 FILTER AND FOR
L
=3
LEVELS OF DECOMPOSITION.
Tables I and IIgive as a function of the sample posi-
tion
for , . The number of samples required
in each subband turns out to be a periodic function of
with pe-
riod
2 . To outline the dependency of from ,
results are provided for 2
successive values of .Asthe
filter used is very short, the number of wavelet coefficients in-
volved in the PW-IDWT is very small. For comparison, Table II