IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 18, NO. 2, FEBRUARY 1999 181
Adaptive Predictive Multiplicative Autoregressive Model for Medical Image Compression Zuo-Dian Chen, Ruey-Feng Chang,* and Wen-Jia Kuo
Abstract—In this paper, an adaptive predictive multiplicative autore- gressive (APMAR) method is proposed for lossless medical image coding.
The adaptive predictor is used for improving the prediction accuracy of encoded image blocks in our proposed method. Each block is first adaptively predicted by one of the seven predictors of the JPEG lossless mode and a local mean predictor. It is clear that the prediction accuracy of an adaptive predictor is better than that of a fixed predictor. Then the residual values are processed by the MAR model with Huffman coding. Comparisons with other methods [MAR, SMAR, adaptive JPEG (AJPEG)] on a series of test images show that our method is suitable for reversible medical image compression.
Index Terms— Image coding, lossless compression, multiplicative au- toregressive model.
I. INTRODUCTION
There are many different kinds of medical digital images such as magnetic resonance (MR), computerized tomography (CT), ul- trasound (US), X-ray images, etc. Medical image compression can be classified into two categories, lossless and lossy methods. For example, differential pulse code modulation (DPCM) [1], hierarchical interpolation (HINT) [2], and multiplicative autoregressive models (MAR) [3]–[6] are lossless methods. The discrete-cosine transform (DCT) [7] and vector quantization (VQ) [8] are lossy methods.
The major advantage of lossless methods is that the image can be reconstructed exactly. This property is especially important for medical images because of legal reasons.
In the MAR approach, pixel values that subtract the block mean bmfrom the original pixel valuesff(i; j)g are predicted by the MAR predictor. Subsequently, the residual values are encoded by entropy coding. Although the correlation of the pixels can be reduced by subtracting the block meanbmbefore applying the MAR predictor, the performance of other predictors to reduce the correlation of the pixels is much better than that of using only the subtraction of block mean bm:
In this paper, a novel approach called adaptive predictive MAR (APMAR) is introduced. Two predictive methods are used in the decorrelation step of our method. First, adaptive prediction is used to process the image blocks. One out of eight predictors, viz., the seven predictors in the JPEG lossless standard [9] and the local mean predictor in the space-varying MAR (SMAR) method [10], is adaptively chosen as the predictor for the first decorrelation step.
Then we use the MAR predictor to further reduce the redundancy of the residual image. Entropy coding is used in the final coding step of APMAR.
The rest of this paper is organized as follows. In Section II, we introduce the JPEG lossless mode and the MAR method. The new compression method is described in Section III. In Section IV the test Manuscript received June 18, 1996; revised December 31, 1998. This work was supported by the National Science Council under Grant NSC85-2213-E- 194-016. The Associate Editor responsible for coordinating the review of this paper and recommending its publication was M. Vannier. Asterisk indicates corresponding author.
Z.-D. Chen, *R.-F. Chang, and W.-J. Kuo are with the Advanced System Integration Laboratory, Department of Computer Science and Information En- gineering, National Chung Cheng University, Chiayi, Taiwan 62107, R.O.C.
Publisher Item Identifier S 0278-0062(99)03153-5.
Fig. 1. The relational position of the four predictive samplesSa; Sb; Sc; andSdand the current predicted sampleSx:
TABLE I
PREDICTORS OF THEEIGHT-PREDICTORSELECTION
images and the experimental results are presented and conclusions are shown.
II. THESEQUENTIAL LOSSLESSMODE OF JPEGAND THE MAR METHOD
A. The Sequential Lossless Mode of JPEG
The JPEG lossless mode uses four neighboring pixelsSa; Sb; Sc; and Sd to predict the value of the current pixel indicated by Sx: Fig. 1 shows the relation of these pixels. The first seven predictors in Table I are the JPEG predictors. In order to reduce the bit rate and improve the image quality, we will select a suitable predictor for each encoded block.
B. Lossless Image Compression Using the MAR Method
There are many different kinds of MAR models. The consideration of choosing a MAR model mainly depends on the tradeoff among predictor performance, the computation complexity of the predictor, and the overhead of the parameters. We have chosen the 3 2 3 nonsymmetric half-plane support (NSHP) model shown in Fig. 2. In this model, the two-dimensional (2-D) polynomial operatora(pv; ph) is
a(pv; ph) = (1 + c1p01h )(1 + c2p01v )(1 + c3p01v ph) wherepmvpnhrepresents a shift ofm units along the vertical direction andn units along the horizontal direction and c1; c2; c3are constant coefficients.
Among the many different algorithms for finding the parameters of the MAR model, a simple, and thus attractive, estimation method is recursive pseudo-linear regression (RPLR). The RPLR method represents the pixel valueu(i; j) by the form u(i; j) = 888(i; j)T222 + e(i; j), where 222 is the unknown parameter matrix, 888(i; j) is the regression matrix, ande(i; j) represents zero-mean white noise. The recursive least squares (RLS) algorithm [11] is used to iteratively update the parameter matrix222 by the following steps:
0278–0062/99$10.001999 IEEE
182 IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 18, NO. 2, FEBRUARY 1999
Fig. 2. The support region of the 32 3 NSHP model, where x is the current predicted pixel.
1) 22(t) = 22 2(t 0 1) + K2 KK(t)(u(i; j) 0 888(i; j)T222(t 0 1));
2) KK(t) = PPK P (t 0 1)888(i; j)(1 + 888(i; j)TPPP (t 0 1)888(i; j))01; 3) PPP (t0) = [III 0 KKK(t)88(i; j)8 T]PPP (t 0 1)
whereKK(t) is a 3 2 1 matrix, PPK P (t) is a 3 2 3 matrix, 22(t) =2 [c1; c2; c3]T is the parameter matrix at iteration t, and 8(i; j) will affect the parameters in the next iteration
8(i; j) =
0u(i; j 0 1)
0u(i 0 1; j) 0 c1u(i 0 1; j 0 1) 0u(i 0 1; j + 1) 0 c1u(i 0 1; j)
0c2u(i 0 2; j + 1) 0 c1c2u(i 0 2; j) :
The recursive algorithm requires the initial values of222(0) and PPP (0) withPPP (0) invertible. A good choice for PPP (0) is generally PPP (0) = cIII, wherec is a small positive scalar. At each iteration, the parameters fc1; c2; c3g must be checked to see whether they are in the stable parameter range for the 3 2 3 NSHP model which is given by Sp = f(c1; c2; c3): jc1j; jc2j; jc3j < 1g:
Hence, the residual value ofu(i; j), which is obtained by subtract- ing the block meanbmfrom the original pixel valuef(i; j), will be predicted by the equation~u(i; j) = (1 0 a(pv; ph))u(i; j): For the 32 3 NSHP MAR model
~u(i; j) = (0c1c2c3pv020 c2c3p02v ph0 c1c2p01v p01h 0 c1c3p01v 0 c2p01v 0 c3pv01ph0 c1p01h )u(i; j)
= 0c1c2c3u(i 0 2; j) 0 c2c3u(i 0 2; j + 1) 0 c1c2u(i 0 1; j 0 1) 0 c1c3u(i 0 1; j) 0 c2u(i 0 1; j) 0 c3u(i 0 1; j + 1) 0 c1u(i; j 0 1):
III. APMAR METHOD
We propose a new compression method, APMAR, in this paper.
Both MAR and SMAR use a fixed prepredictor to reduce the correlation of the pixels before the actual prediction of the MAR model. In SMAR [10], the local mean(Sa+ Sb+ Sc+ Sd)=4 is first subtracted from the pixel value for each pixel. Then, the residual values are encoded by the MAR method. It is not appropriate to use a fixed predictor for different medical images because of the widely varying image characteristics. Hence, we admit all seven predictors of lossless JPEG to reduce the correlation of the pixels. The local mean predictor(Sa+ Sb+ Sc+ Sd)=4 is added as an option without increasing the additional overhead because three bits are already needed for the seven original JPEG predictors. We will refer to this prepredictive process as eight-predictor selection.
The first step of APMAR is to split the image into blocks.
Table I shows the eight predictors of eight-predictor selection. The one with the lowest total error between original pixel values and predicted values in a block is adopted as the predictor of the block. After choosing the predictor, MAR is used for predicting the residual values. Finally, the index of the eight-predictor selection, the parameters of the MAR model, and the final residual values for each
TABLE II
THESUMMATION OFENTROPY ANDOVERHEAD FORAPMAR
ANDOTHERLOSSLESSMETHODS FOR THETESTIMAGES
block are transmitted to the receiver. The encoder of our APMAR method is shown in Fig. 3.
The steps of the APMAR compression method are as follows.
1) Step 1: Splitting the image into blocks.
2) Step 2: Using the eight-predictor selection to reduce the cor- relation of pixels.
3) Step 3: Using MAR prediction to further reduce the correlation of the residual values.
4) Step 4: Encoding the final residual values by entropy coding.
IV. EXPERIMENTS AND RESULTS
Experiments to evaluate the proposed method have been performed on a number of 2-D images or slices of three-dimensional (3-D) images (see Fig. 4 for an overview).
The experiments showed that the best performance was obtained when the image was segmented into 64 2 64 blocks for a 512 2 512 image or 32 2 32 blocks for a 256 2 256 image in the 3 2 3 NSHP MAR model. Three bits are used for the eight-predictor selection in APMAR. The overhead of the MAR prediction is twenty- four bits for each block where each parameter of the MAR model is quantized with eight bits. The block mean for the original MAR method and the background pixel valuebgm are all quantized with eight bits.
The RLS algorithm requires initial values222(0) and PPP (0): The initial value of matrixPPP must be a small positive matrix. In our experiments, the initial parameters c1; c2; c3 of the 3 2 3 NSHP model for each block are set to be 0.01 andPPP (0) is set to be 0:01III whereIII is the identity matrix.
A. Experimental Results of APMAR
The summation of entropy and overhead for APMAR is compared with that of MAR, SMAR, and AJPEG. The entropy is defined as
H = 0 k
i=1
pi1 log2pi
wherepiis the probability of gray leveli in the image and k is the largest gray level.
The experimental results are listed in Table II. We find that APMAR is more suitable for reversibly compressing different medical images than the other lossless compression methods. Although in individual cases MAR or AJPEG may accomplish a slightly lower en- tropy than APMAR, the overall conclusion is that APMAR generally outperforms the other reversible compression techniques significantly.
REFERENCES
[1] K. Chen and T. V. Ramabadran, “Near-lossless compression of medical images through entropy-coded DPCM,” IEEE Trans. Med. Imag., vol.
13, pp. 538–548, Sept. 1994.
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 18, NO. 2, FEBRUARY 1999 183
Fig. 3. Encoder of the APMAR.
(a) (b) (c)
(d) (e)
(f) (g)
(h) (i)
Fig. 4. Overview of the test images. (a) Brain A: sagittal slice (2562 256) of an MR image. (b) Brain B: idem. (c) Brain C: idem. (d) Knee-A: idem.
(e) Knee-B: idem. (f) Chest: transaxial slice (5122 512) CT image. (g) Heart: 512 2 480 US image. (h) Angio: 512 2 480 angiographic image. (i) Chest: 5122 512 X-ray image. All images were quantized in 8 bits/pixel.
184 IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 18, NO. 2, FEBRUARY 1999
[2] M. A. Viergever and P. Roos, “Hierarchical interpolation,” IEEE Eng.
Med. Biol. Mag., vol. 12, pp. 48–55, Mar. 1993.
[3] M. Das and N. K. Loh, “New studies on adaptive predictive coding of images using multiplicative autoregressive models,” IEEE Trans. Image Processing, vol. 1, pp. 106–111, Jan. 1992.
[4] S. Burgett and M. Das, “Predictive image coding using multiresolution multiplicative autogressive models,” Proc. Int. Elect. Eng., vol. 140, no.
2, pp. 127–134, Apr. 1993.
[5] M. Das and S. Burgett, “Lossless compression of medical images using two-dimensional multiplicative autoregressive models,” IEEE Trans.
Med. Imag., vol. 12, pp. 721–726, Dec. 1993.
[6] M. Das and C. L. Li, “Simple space-varying least squares model for lossless medical image compression,” Electron. Lett., vol. 30, no. 11, pp. 849–850, May 1994.
[7] M. S. Rebelo, S. S. Furuie, A. C. L. Munhoz, L. Moura, and C. P.
Melo, “Assessment of a lossy image compression technique applied to cardiac images in nuclear medicine,” in Proc. Computers Cardiology, Sept. 1993, pp. 827–830.
[8] P. C. Cosman, C. Tseng, R. M. Gray, R. A. Olshen, L. E. Moses, H. C. Davidson, C. J. Bergin, and E. A. Riskin, “Tree-structured vector quantization of CT chest scans: Image quality and diagnostic accuracy,”
IEEE Trans. Med. Imag., vol. 12, pp. 727–739, Dec. 1993.
[9] W. B. Pennebaker and J. L. Mitchell, JPEG Still Image Data Compres- sion Standard. New York: Van Nostrand Reinhold, 1993.
[10] M. Das, D. L. Neuhoff, and C. L. Lin, “Near-lossless compression of medical images,” in Proc. ICASSP, 1995, pp. 2347–2350.
[11] L. Ljung and T. Soderstrom, Theory and Practice of Recursive Identifi- cation. Cambridge, MA: MIT, 1983.
Correction to “Performance Characteristics of a Compact Position-Sensitive LSO Detector Module”
J. J. Vaquero, J. Seidel, S. Siegel, W. R. Gandler, and M. V. Green*
In Section III (pp. 971–972) of the above paper,1the three italicized paragraphs were incorrectly moved to Section IV, p. 973. The corrected text is reprinted in the following.
We apologize to the authors and readers for this error.
III. RESULTS
Single photon field flood images acquired at 30, 140, and 511 keV, and count profiles along the central row of each image, are shown in Fig. 3. These images were created with the centroid algorithm (1) and are corrected for the LSO background. The least restrictive event selection criteria were applied during the creation of these images.
A version of Fig. 3(c), enhanced to reveal faint structures, is shown in Fig. 4(a). Each crystal in this image appears to be joined to its eight neighbors by faint horizontal, vertical, and diagonal straight lines in a “connect-the-dot” pattern. The details of this pattern change depending on whether the LSO array is illuminated from the front [Fig. 4(a)] or from the left side [Fig. 4(b)].
Manuscript received March 4, 1999. Asterisk indicates corresponding au- thor.
J. J. Vaquero, J. Seidel, S. Siegel, and W. R. Gandler are with the Nuclear Medicine Department, National Institutes of Health, Bethesda, MD 20892 USA.
*M. V. Green is with the Nuclear Medicine Department, National Institutes of Health, Building 10, Room 1C401 MSC 1180, Bethesda, MD 20892 USA (e-mail: mike@nmdhst.cc.nih.gov).
Publisher Item Identifier S 0278-0062(99)04761-8.
1J. J. Vaquero, J. Seidel, S. Siegel, W. R. Gandler, and M. V. Green, IEEE Trans. Med. Imag., vol. 17, pp. 967–978, Dec. 1998.
An LSO background image acquired with no external sources present, and a count profile along the central row, are shown in Fig. 5. The least restrictive event selection criteria were applied during creation of this image.
Average apparent crystal widths at the half maximum and tenth maximum levels, expressed as a fraction of the actual crystal width (2 mm), are listed in Table I for crystals in the central row of Fig. 3.
Average widths, and the largest and smallest widths in the row, are included for images created from all data within the dynamic range of the ADC’s. Average widths calculated from these same image data when subjected to energy windowing (W) are also shown. At 140 keV the energy window was 140 keV 620 keV, and at 511 keV the window was 511 keV 6 70 keV. Energy windowing was not applied to the I-125 field flood data since nearly all 30-keV photon interactions in LSO are photoelectric and windowing is unnecessary.
The fraction of detected coincidence events that occur within each of three adjacent crystals scanned by the collimated F-18 line source is shown as a function of source position for different exclusion conditions in Fig. 6(a). The dotted curve portrays the fraction of events that occur in ROI’s the size of an LSO crystal (L) when no energy exclusion condition is applied. The dashed curve portrays this fraction for events within the L-ROI and in the 511-keV6 70-keV energy range. The continuous curve portrays the fraction of events assigned to a crystal when events outside small (S) ROI’s drawn around each apparent crystal are excluded, as well as events outside the energy window.
Accuracy and sensitivity are portrayed in the bar graph shown in Fig. 6(b) for different energy and spatial exclusion conditions when the scanning source is positioned directly over the center of the middle crystal in Fig. 6(a). For L-ROI’s, energy exclusion conditions alone reduce the fraction of accepted events from 100% to 55%, while the S spatial exclusion condition combined with the energy exclusion condition reduces the fraction of accepted events further to 34%. The accuracy of event positioning increases from 72% to 95% as these increasingly selective criteria are imposed.
The apparent separation between peaks along the main diagonal of Fig. 3(c) is plotted in Fig. 7(a) where the actual diagonal crystal separation is 3.07 mm. The largest deviation of peak separation from the mean peak separation was 0.82 mm, while the average deviation from the mean peak separation was 0.6 mm. The sinusoidal- like variation in spacing between peaks in Fig. 7(a) is detectable in Figs. 3, 4, and 5(a).
Fig. 7(b) shows apparent module gain before correction as a func- tion of crystal location. Gain is defined as the normalized channel number in which the photopeak maximum occurs for each crystal.
Energy spectra for three different incident photon energies are shown in Fig. 8 for the same central crystal. These spectra have been corrected for the LSO background, whose spectrum is shown in Fig. 8(c) The fraction of LSO background events occurring in the energy range 1406 20 keV [Fig. 8(c)] was 1.2% of the total LSO background rate or about 12 counts/s.
Energy resolution and its variation within the UFOV are listed in Table II for several different incident photon energies. For illustrative purposes, energy resolution at 511 keV in an identical detector module composed of 2 mm2 2 mm 2 10-mm-long BGO crystals rather than LSO crystals is also included.
The TAC spectrum obtained with the detector module in time coincidence with a second LSO detector is shown in Fig. 9. The second timing peak was acquired with an interposed time delay of 4 ns. The average FWHM of these peaks was 1.2 ns.
0278–0062/99$10.001999 IEEE
