ABSTRACT
Image Coding and Data Compression
High spatial resolution and ne grayscale quantization are often required
in biomedical imaging Digital mammograms are typically represented in
arrays of pixels with bpixel leading to rawdata les of
the order of MB per image Volumetric data obtained by CT and MRI
could be of size voxels with bvoxel occupying MB
per examination Patients with undetermined or multiple complications may
undergo several examinations via dierent modalities such as Xray imaging
ultrasound scanning CT scanning and nuclear medicine imaging resulting
in large collections of image les
Most healthcare jurisdictions require medical records including images of
adults to be stored for durations of the order of seven years from the date
of acquisition Children s records and images are required to be maintained
until at least the time they reach adulthood
With the view to improve the eciency of storage and access several imag
ing centers and hospitals have moved away from lmbased storage toward
electronic storage Furthermore most medical imaging systems have moved
to direct digital image acquisition with adequate resolution putting aside the
debate on the quality of an original lmbased image versus that of its scanned
digitized representation Since an entire series of conferences has been
dedicated to PACS see the PACS volumes of the SPIE Medical Imaging con
ference series Networks and systems for PACS are integrated into the
infrastructure of most modern hospitals The major advantages and disad
vantages of digital and lmbased archival systems are listed below
Films deteriorate with age and handling Digital images are unaected
by these factors
Despite elaborate indexing schemes lms tend to get lost or misplaced
Digital image les are less likely to face these problems
Digital image les may be accessed simultaneously by several users
Although multiple copies of lmbased images may be made it would
be an expensive option that adds storage and handling complexities
With the proliferation of computers digital images may be viewed and
manipulated at several convenient locations including a surgical suite
a patient s bedside and one s home or oce Viewing lmbased im
ages with detailed attention requires specialized viewing consoles under
controlled lighting conditions
Digital PACS require signicant initial capital outlay as well as routine
maintenance and upgrading of the computer storage and communi
cation systems However these costs may be oset by the savings in
the continuing costs of lm as well as the associated chemical process
ing systems and disposal The environmental concerns related to lm
processing are also removed by digital PACS
Digital images may be compressed via image coding and data compres
sion techniques so as to occupy less storage space
The nal point above forms the topic of the present chapter
Although the discussion above has been in the context of image storage or
archival similar concerns regarding the size of image les and the desirability
of compression arise in the communication of image data In this chapter
we shall study the basic concepts of information theory that apply to image
coding compression and communication We shall investigate several tech
niques for encoding image data including decorrelation procedures to modify
the statistical characteristics of the data so as to permit ecient representa
tion coding and compression
The representation of the signicant aspects of an image in terms of a small
number of numerical features for the purpose of pattern classication may
also be viewed as image coding or data compression however we shall treat
this topic separately see Chapter
Considerations Based on Information Theory
Image data compression is possible due to the following basic characteristics
Code redundancy all code words pixel values do not occur with
equal probability
Spatial redundancy the values of neighboring pixels tend to lie within
a small dynamic range and exhibit a high level of correlation
Psychovisual redundancy human analysts can recognize the essential
nature and components of an image from severely reduced versions such
as caricatures edges and regions and need not or do not pay attention
to precise numerical values
Informationtheoretic considerations are based upon the notion of informa
tion as related to the statistical uncertainty of the occurrence of an event
such as a signal an image or a pixel value rather than the structural sym
bolic pictorial semantic or diagnostic content of the entity The measure of
entropy is based upon the probabilities of occurrence of the various symbols
involved in the representation of a message or image see Section Despite
the mathematical and theoretical powers of measures such as entropy the
standpoint of viewing an image as being composed of discrete and indepen
dent symbols numerical values removes the analyst from the realworld and
physical properties of the image The use of the underlying assumptions also
lead to severe limitations in entropybased source coding with lossless com
pression factors often limited to the order of Additional techniques based
upon decorrelation of the image data via the identication and modeling of
the underlying imagegeneration phenomena or the use of pattern recognition
techniques could assist in improving the performance of image compression
procedures
Noiseless coding theorem for binary transmission
Given a code with an alphabet of two symbols and a source A with an alphabet
of two symbols the average length of the code words per source symbol may
be made arbitrarily close to the lower bound entropy HA by encoding
sequences of source symbols instead of encoding individual symbols
The average length Ln of encoded nsymbol sequences is bounded by
HA
Ln
n
HA
n
Diculties exist in estimating the true entropy of a source due to the fact
that pixels are statistically dependent that is correlated from pixel to pixel
row to row and frame to frame of reallife images The computation of the
true entropy requires that symbols be considered in blocks over which the
statistical dependence is negligible In practice this would translate to esti
mating joint PDFs of excessively long vectors Values of entropy estimated
with single pixels or small blocks of pixels would result in overestimates of
the source entropy If blocks of pixels are chosen such that the sequence
entropy estimates converge rapidly to the limit then blockcoding methods
may provide results close to the minimum length given by Equation
Runlength coding may be viewed as an adaptive blockcoding technique see
Section
Lossy versus lossless compression
A coding or compression method is considered to be lossless if the original
image data can be recovered with no error from the coded and compressed
data Such a technique may also be referred to as a reversible bitpreserving
or errorfree compression technique
A compression technique becomes lossy or irreversible if the original data
cannot be recovered with complete pixelbypixel numerical accuracy from
the compressed data In the case of images the human visual system can tol
erate signicant numerical dierences or error in the sense that the degraded
image recovered from the compressed data is perceived to be essentially the
same as the original image This arises from the fact that a human ob
server will typically not examine the numerical values of individual pixels
but instead assess the semantic or pictorial information conveyed by the data
Furthermore a human analyst may tolerate more error noise or distortion in
the uniform areas of an image than around its edges that attract visual atten
tion Data compression techniques may be designed to exploit these aspects
to gain signicant advantages in terms of highly compressed representation
with high levels of loss of numerical accuracy while remaining perceptually
lossless On the same token in medical imaging if the numerical errors in the
retrieved and reconstructed images do not cause any change in the diagnostic
results obtained by using the degraded images one could achieve high levels
of numerically lossy compression while remaining diagnostically lossless
In the quest to push the limits of numerically lossy compression techniques
while remaining practically lossless under some criterion the question arises as
to the worth of such practice Medical practice in the present highly litigious
society could face large nancial penalties and loss due to errors Radiologi
cal diagnosis is often based upon the detection of minor deviations from the
normal or average patterns expected in medical images If a lossy data com
pression technique were to cause such a faint deviation to be less perceptible
in the compressed and reconstructed image than in the original image and
the diagnosis based upon the reconstructed image were to be in error the
nancial compensation to be paid would cost several times the amount saved
in data storage the loss in professional standing and public condence could
be even more damaging In addition dening the delity of representation in
terms of the closeness to the original image or distortion measures is a di
cult and evasive activity Given the high levels of the professional care and
concern as well as the scal and emotional investment that are part of medi
cal image acquisition procedures it would be undesirable to use a subsequent
procedure that could cause any degradation of the image In this spirit only
lossless coding and compression techniques will be described in the present
chapter Regardless it should be noted that any lossy compression technique
may be made lossless by providing the numerical error between the original
image and the degraded image reconstructed from the compressed data Al
though this step will lead to additional storage or transmission requirements
the approach can facilitate the rapid retrieval or transmission of an initial
lowquality image followed by completely lossless recovery such a procedure
is known as progressive transmission especially when performed over multiple
stages of image quality or delity
Distortion measures and delity criteria
Although we have stated our interest in lossless coding of biomedical images
other processes such as the transmission of large quantities of data over noisy
channels may lead to some errors in the received images Hence it would be
relevant to consider the characterization of the distortion so introduced and
analyze the delity of the received image with respect to the original
The binary symmetric channel is characterized by a single parameter the
biterror probability p see Figure The channel capacity is given by
C p log p q log q
where q p
FIGURE
Transmission error probabilities in a binary symmetric channel
The leastsquares singleletter delity criterion is dened as
n
xy
n
n
X
l
x
l
y
l
l
where x and y are the transmitted and received nbit vectors blocks or
words respectively
The Hamming distance between the vectors x and y is dened as
D
H
xy
n
n
X
l
x
l
y
l
Measures of delity may also be dened based upon entire images by den
ing an error image as
emn gmn fmn
where gmn is the received degraded version of the original transmitted
image fmn and then dening the RMS value of the error as
e
RMS
v
u
u
t
N
N
X
m
N
X
n
gmn fmn
or SNR as
SNR
P
N
m
P
N
n
g
mn
P
N
m
P
N
n
e
mn
See Section for more details on measures of SNR
Fundamental Concepts of Coding
In general coding could be dened as the use of symbols to represent infor
mation The following list provides the denitions of a few basic terms and
concepts related to coding
An alphabet is a predened set of symbols
A word is a nite sequence of symbols from an alphabet
A code is a mapping of words from a source alphabet into the words of
a code alphabet
A code is said to be distinct if each code word is distinguishable from
the other code words
A distinct code is uniquely decodable if every code word is identiable
when immersed in a sequence of code words with no separators between
the words
A desirable property of a uniquely decodable code is that it be decodable
on a wordtoword basis This is ensured if no code word may be a prex
to another the code is then instantaneously decodable
A code is said to be optimal if it is instantaneously decodable and has
the minimum average length for a given source PDF
Examples of symbols are f g in the binary alphabet f g
in the octal system f g in the decimal system f
ABCDE Fg in the hexadecimal system fI V X LC D
Mg in the Roman system with the decimal equivalents of
and respectively and fA Z a zg in the English alphabet not
considering punctuation marks and special symbols An example of a word
in the context of image coding is in b binary coding standing
for the gray level in the decimal system Table lists the codes for
integers in the range in the Roman decimal binary Gray octal
and hexadecimal codes The Gray code has the advantageous feature
that only one digit is changed from one number to the next Observe that in
general all of the codes described here including the English language fail
the conditions dened above for an optimal code
Direct Source Coding
Pixels generated by reallife sources of images bear limitations in dynamic
range and variability within a small spatial neighborhood Therefore codes
used to represent pixel data at the source may be expected to demonstrate
certain patterns of limited variation and high correlation Furthermore real
life sources of images do not generate random uncorrelated values that are
equally likely instead it is common to encounter PDFs of gray levels that
are nonuniform Some of these characteristics may be exploited to achieve
ecient representation of images by designing coding systems tuned to specic
properties of the source Because the coding method is applied directly to pixel
values generated by the source without processing them by an algorithm to
generate a dierent series of values such techniques are categorized as direct
source coding procedures
Human coding
Human proposed a coding system to exploit the occurrence of some
pixel values with higher probabilities than other pixels The basic idea in
Human coding is to use short code words for values with high probabilities of
occurrence and longer code words to represent values with lower probabilities
of occurrence This implies that the code words used will be of variable length
the method also presumes prior knowledge of the PDF of the source symbols
gray levels It is required that the code words be uniquely decodable on
a wordbyword basis which implies that no code word may be a prex to
another Human devised a coding scheme to meet these requirements and
lead to average codeword lengths lower than that provided by xedlength
codes Human coding provides an average codeword length L that is limited
by the zerothorder entropy of the source H
see Equation and H
H
L H
The procedure to generate the Human code is as follows
TABLE
Integers in the Range in Several Alphabets or Codes
English Portuguese Roman Decimal Binary Gray Octal Hex
Zero Zero
One UnUma I
Two DoisDuas II
Three Tres III
Four Quatro IV
Five Cinco V
Six Seis VI
Seven Sete VII
Eight Oito VIII
Nine Nove IX
Ten Dez X A
Eleven Onze XI B
Twelve Doze XII C
Thirteen Treze XIII D
Fourteen Catorze XIV E
Fifteen Quinze XV F
Sixteen Dezesseis XVI
Seventeen Dezessete XVII
Eighteen Dezoito XVIII
Nineteen Dezenove XIX
Twenty Vinte XX
Leading zeros have been removed in the decimal and hexadecimal Hex codes
but retained in the binary Gray and octal codes
Prepare a table listing the symbols gray levels in the source image
sorted in decreasing order of the probabilities of their occurrence
Combine the last two probabilities The list of probabilities now has
one less entry than before
Copy the reduced list over to a new column rearranging as necessary
such that the probabilities are in decreasing order
Repeat the procedure above until the list of probabilities is reduced to
only two entries
Assign the code digits and to the two entries in the nal column
of probabilities Note There are two possibilities of this assignment
that will lead to two dierent codes however their performance will be
identical
Working backwards through the columns of probabilities assign addi
tional bits of and to the two entries that resulted in the last com
pounded entry in the column
Repeat the procedure until the rst column of probabilities is reached
and all symbols have been assigned a code word
It should be noted that a Human code is optimal for only the given source
PDF a change in the source PDF would require the design of a dierent code
in order to be optimal A disadvantage of the Human code is the increasing
length of its code words especially for sources with several symbols The
method does not perform any decorrelation of the data and is limited in
average codeword length by the zerothorder entropy of the source
Example Figure shows a part of the image in Figure a
quantized to bpixel The gray levels in the image are in the range
and would require bpixel with straight binary coding The histogram of
the image is shown in Figure it is evident that some of the pixel values
occur with low probabilities
The procedure for accumulating the probabilities of occurrence of the source
symbols is illustrated in Figure The Human coding process is shown in
Figure Note that a dierent code with equivalent performance may be
generated by reversing the order of assignment of the code symbols and
at each step The average codeword length is bpixel which is slightly
above the zerothorder entropy of b of the image The advantage is
relatively small due to the fact that the source in the example uses only eight
symbols with bpixel and has a relatively wellspread histogram PDF
However simple representation of the data using ASCII coding would require
a minimum of bpixel the savings with reference to this requirement are
signicant Larger advantages may be gained by Human coding of sources
with more symbols and narrow PDFs
FIGURE
Top to bottom A part of the image in Figure a quantized to
bpixel shown as an image a D array and as a string of integers with
the graylevel values of every pixel The line breaks in the string format have
been included only for the sake of printing within the width of the page
FIGURE
Graylevel histogram of the image in Figure Zerothorder entropy H
b
B i o m e d i c a l I m a g e A n a l y s i s
FIGURE
Accumulation of probabilities Prob of occurrence of gray levels in the derivation of the Human code for the image in
Figure The probabilities of occurrence of the symbols have been rounded to two decimal places and add up to unity
I m a g e C o d i n g a n d D a t a C o m p r e s s i o n
FIGURE
Steps in the derivation of the Human code for the image in Figure Prob probabilities of occurrence of the gray
levels The binary words in bold italics are the Human code words at the various stages of their derivation See also
Figure
Interpixel correlation may be taken into account in Human coding by
considering combinations of pixels gray levels as symbols If we were to
consider pairs of gray levels in the example above with gray levels quantized to
bpixel we would have a total of possibilities see Table The
rstorder entropy of the image considering pairs of gray levels isH
b
an average codeword length close to this value may be expected if Human
coding is applied to pairs of gray levels
TABLE
Counts of Occurrence of Pairs of Pixels in the Image in Figure
Current pixel Next pixel in the same row
For example the pair occurs times in the image The last pixel in
each row was not paired with any pixel The rstorder entropy of the image
considering the probabilities of occurrence of pairs of graylevel values as in
Equation is H
b The zerothorder entropy is H
b
Although the performance of Human coding is limited when applied di
rectly to source symbols the method may be applied to decorrelated data with
signicant advantage due to the highly nonuniform or concentrated PDFs of
such data The performance of Human coding as a postencoder following
decorrelation methods is discussed in several sections to follow
Runlength coding
Images with high levels of correlation may be expected to contain strings of
repeated occurrences of the same gray level such strings are known as runs
Data compression may be achieved by coding such runs of gray levels For
example the rst three rows of the image in Figure may be represented
as follows
Row
Row
Row
In the code above each pair of values represents a run with the rst value
standing for the gray level and the second value giving the number of times
the value has occurred in the run The coding procedure is interrupted at the
end of each row to permit synchronization in case of errors in the reception
and decoding of run values Direct coding of the pixels in the
three rows considered at bpixel leads to a total code length of b In
the runlength code given above if we were to use three bits per gray level and
four bits per runlength value we get a total code length of b
The savings are small in this case due to the busy nature of the image
which represents an eye
Runlength coding is best suited for the compression of bilevel images where
long runs may be expected of the two symbols and Images with ne
details intricate texture and highresolution quantization with large numbers
of bits per pixel may not present long runs of the same gray level runlength
coding in such cases may lead to data expansion rather than compression
This is the case with the image in Figure past the third row
Runlength coding may be advantageously applied to bit planes of gray
level and color images The use of Gray coding see Table improves the
chances of long runs in the bit planes due to the feature that the Gray code
changes in only one bit from one numerical value to the next
Errors in run length could cause severe degradation of the reconstructed
image due to the loss of pixel position Synchronization at the end of each
row can avoid the carrying over of such errors beyond the aected row
Runs may also be dened over D areas However images with ne details
do not present such uniform areas with large numbers of occurrence to lend
much coding advantage
Arithmetic coding
Arithmetic coding is a family of codes that treat input symbols as mag
nitudes Shannon presented the basic idea of representing a string of
symbols as the sum of their scaled probabilities Most of the development
towards practical arithmetic coding has been due to Langdon and Rissa
nen The basic advantage of arithmetic coding over Human
coding is that it does not suer by the limitation that each symbol should
have a unique code word that is at least one bit in length
The mechanism of arithmetic coding is illustrated in Figure
The symbols of the source string are represented by their individual probabili
ties p
l
and cumulative probabilities sum of the probabilities of all symbols up
to but not including the current symbol P
l
At any given stage in coding
the source string is represented by a code point C
k
and an interval A
k
The
code point C
k
represents the cumulative probability of the current symbol
on a scale of interval size A
k
A new symbol being appended to the source
string is encoded by scaling the interval by the probability of the current
symbol as
A
k
A
k
p
l
and dening the new code point as
C
k
C
k
A
k
P
l
Decoding is performed by the reverse of the procedure described above The
interval A
is initialized to unity and the current code point is determined by
the range into which the nal code point C
nal
falls The scaling of the interval
and code point is performed as during encoding The encoding procedure
ensures that no future code point C
k
exceeds the current value C
k
A
k
Thus
a carry over to a given bit position in the binary representation of C
nal
occurs at most once during encoding This fact is made use of for incremental
coding and for using niteprecision arithmetic Finite precision is used by
employing a technique known as bit stung where if a series of ones
longer than the specied precision occurs in the binaryfraction representation
of C
k
a zero is inserted this ensures that further carries do not propagate into
the series of ones Witten et al provide an implementation of arithmetic
coding using integer arithmetic
Direct arithmetic coding of an image consists of an initial estimation of the
probabilities of the gray values in the image followed by rowwise arithmetic
coding of pixels Direct coding does not take into account the correlation
between adjacent pixels Arithmetic coding can be modied to make use of the
correlation between pixels to some extent by using conditional probabilities
of occurrence of gray levels
In a version of arithmetic coding known as Qcoding the individual
bit planes of an image are coded using probabilities conditioned on the sur
rounding bits in the same plane as the context A more ecient procedure is
to perform decorrelation of the pixels of the image separately and to use the
basic arithmetic coder as a postencoder on the decorrelated set of symbols
see for example Rabbani and Jones The performance of arithmetic
coding as a postencoder after the application of decorrelation methods is
discussed in several sections to follow
FIGURE
Arithmetic coding procedure The range A
is initialized to unity Each sym
bol is represented by its individual probability p
l
and cumulative probability
P
l
The string being encoded is represented by the code point C
k
on the cur
rent range A
k
The range is scaled down by the probability p
l
of the current
symbol and the process is repeated One symbol is reserved for the end of
the string Figure courtesy of GR Kuduvalli
Example The symbols used to the represent the image in Figure and
their individual as well as cumulative probabilities are listed in Table
The intervals representing the symbols are also provided in the table Let us
consider the rst three symbols f g in the fth row of the image The
procedure to derive the arithmetic code for this string of symbols is shown in
Figure
TABLE
The Symbols Used in the Image in
Figure Along with Their Individual
Probabilities of Occurrence p
l
Cumulative
Probabilities P
l
and Intervals Used in
Arithmetic Coding
Symbol l Count p
l
P
l
Interval
The initial code point is C
and the initial interval is A
When
the rst symbol is encountered the code point and interval are updated
as C
C
A
P
A
A
p
For the
next symbol we get C
C
A
P
and
A
A
p
With the third symbol appended to
the string we have C
C
A
P
and A
A
p
The code points have been given in decimal code to full precision as re
quired in this example the individual code probabilities have been rounded
to two decimal places In actual application the code points need to be rep
resented in binary code with nite precision The average codeword length
FIGURE
The arithmetic coding procedure applied to the string f g formed by the
rst three symbols in the fth row of the image in Figure See Table
for the related probabilities and intervals see also Figure All intervals
are shown mapped to the same physical length although their true values
decrease from the interval at the top of the gure to that at the bottom
The numerals in italics indicate the symbols gray levels in the
image The values in bold at the ends of each interval give the values of C
k
and C
k
A
k
at the corresponding stage of coding
per symbol is reduced by encoding long strings of symbols such as an entire
row of pixels in the given image
LempelZiv coding
Ziv and Lempel proposed a universal coding scheme for encoding symbols
from a discrete source when their probabilities of occurrence are not known
a priori The coding scheme consists of a rule for parsing strings of symbols
from the source into substrings or words and mapping the substrings into
uniquely decipherable code words of xed length Thus unlike the Human
code where codes of xed length are mapped into variablelength codes the
LempelZiv code maps codes of variable length corresponding to symbol
strings of variable length into codes of xed length
The LempelZiv coding scheme is illustrated in Figure The
coding procedure starts with a buer of length
n L
s
L
s
where L
s
is the maximum length of the input symbol strings being parsed
is the cardinality of the symbol source in the case of image coding the
number of possible gray levels and is chosen such that The
buer is initially lled with n L
s
zeros and the rst L
s
symbols from the
source The buer is then parsed for the string whose length l
k
is less than
L
s
but is the maximum of all such strings from to n L
s
and which
has an identical string in the buer starting at position n L
s
The code to
be mapped consists of the beginning position b
k
of this string in the buer
from position to nL
s
the length of the string l
k
and the last symbol
s
k
following the end of the string The total length of the code for a straight
binary representation is
l dlog
n L
s
log
L
s
log
e
where dxe is the smallest integer x After coding the string the buer is
advanced by l
k
number of symbols Ziv and Lempel showed that as the
total length of the input symbols tends to the average bit rate for coding
the string approaches that of an optimal code with complete knowledge of the
statistics of the source
The LempelZiv coding procedure may be viewed as a search through a
xedsize variablecontent dictionary for words that match the current string
A modication of this procedure known as LempelZivWelch LZW cod
ing consists of using a variablesized dictionary with every new string
encountered in the source string added to the dictionary The dictionary is
initialized to singlesymbol strings made up of the entire symbol set This
eliminates the need for including the symbol s
k
in the code words The LZW
string table has the prex property for every string of symbols in the table
its prex is also present in the table
The LempelZiv coding procedure At each iteration the buer is scanned
for strings of length l
k
L
s
for a match in the substring of length n L
s
within the buer The matched string location b
k
is encoded and transmitted
Figure courtesy of GR Kuduvalli
Kuduvalli implemented a slight variation of the LZW code in which
the rst symbol of the current string is appended as the last symbol of the
previously parsed string and the new string is added to the string table
With this method the decoded strings are generated in the same order as
the encoded strings The string table itself is addressed during the encoding
procedure as a linklist Each string contains the address of every other string
of which it is a prex Such a linklist is not necessary during decoding
because the addresses of the strings are directly available to the decoder
LZW coding may be applied directly for source coding or applied to decor
related data The following sections provide examples of application of the
LZW code
Example Let us consider again the image in Figure The image
has eight symbols in the range each of which will be an item in the
LZW table shown in Table the eight basic symbols may be represented
by their own code Index and Index represent two possibilities of coding
Consider the string f g which occurs ve times in the image In order
to exploit this feature we need to add the strings f g and f g to the
table we could use the codes for the former and or for the latter
represents the symbol being appended to the string represented by
the code The string f g present at the beginning of the fourth
row in the image may be represented as In this manner long strings of
symbols are encoded with short code words The code index in the dictionary
table is used to represent the symbol string A predened limit may be
applied to the length of the dictionary
TABLE
Development of the
LempelZivWelch LZW
Code Table for the Image
in Figure
String Index Index
Contour coding
Given that a digital image includes only a nite number of gray levels we
could expect strings of the same values to occur in some form of D contours or
patterns in the image The same expectation may be arrived at if we were
to consider the gray level to represent height an image may then be viewed
as a relief map with isointensity contours representing steps or plateaus as
in elevation maps of mountains Information related to all such contours
may then be used to encode the image Although the idea is appealing in
principle ne quantization could result in low probabilities of occurrence of
large contours with simple patterns
Each isointensity contour would require the encoding of the coordinates of
the starting point of the contour the associated gray level and the sequence
of steps movement needed to trace the contour A consistent rule is required
for repeatable tracing of contours The leftmostlooking rule for tracing
a contour is as follows
Select a pixel that is not already a member of a contour
Look at the pixel to the left relative to the direction of entry to the current
pixel on the contour
If the pixel has the same gray level as the current pixel move to the pixel
and encode the type of movement
If not look at the pixel straight ahead
If the pixel has the same gray level as the current pixel move to the pixel
and encode the type of movement
If not look at the pixel to the right
If the pixel has the same gray level as the current pixel move to the pixel
and encode the type of movement
If not move to the pixel behind the current pixel
Repeat the procedure will trace a closed loop and return to the starting
point
Repeat the procedure until every pixel in the image belongs to a contour
The movements allowed are only to the left straight ahead to the right
and back the four possibilities may be encoded using the Freeman chain
code illustrated in Figure
Example The image in Figure is shown in Figure along with the
tracings of three isointensity contours With reference to the contour with
the value at the topleft corner of the image observe that several spatially
connected pixels with the same value and lying within the contour have not
been included in the contour these pixels require additional contours The
contourcoding procedure needs to be applied repeatedly until every pixel in
the image belongs to a contour The encoding of short strings or isolated
occurrences of pixel values could require several coding steps and lead to
increased code length
The data required to represent the contour with the gray level starting
with the pixel in the rst row and rst column would be as follows
Initial point Coordinates Gray level
Freeman code
Contour code requirement four bits for each coordinate three bits for the
pixel value two bits per Freeman code Total b
Direct binary code requirement for the pixels on the contour at bpixel
b at bpixel b
The data required to represent the contour with the gray level starting
with the pixel in the fth row and eighth column would be as follows
Initial point Coordinates Gray level
Freeman code
Contour code requirement Total b Direct binary code requirement for
the eight pixels on the contour at bpixel b at bpixel b
The data required to represent the contour with the gray level starting
with the pixel in the ninth row and rst column would be as follows
Initial point Coordinates Gray level
Freeman code
Contour code requirement b Direct binary code requirement for the
pixels on the contour at bpixel b at bpixel b
It is evident that higher advantages may be gained if a number of long
contours with simple patterns are present in the image
Application Source Coding of Digitized Mammo
grams
Kuduvalli et al applied several coding techniques for the com
pression of digitized mammograms In their work lm mammograms were
scanned using an Eikonix digitizing camera with a linear CCD array to
provide a horizontal scan line of pixels A vertical array size of
pixels was achieved by stepping the array over scan lines A Gordon
Instruments Plannar light box was used to illuminate the Xray lm be
ing digitized The gain and oset variations between the CCD elements were
corrected for in the camera and the data were transferred to the host com
puter over an IEEE bus Corrections were applied for the lightintensity
variations in the Plannar light source used to illuminate the lms and
the digitized image was stored in a Megascan FDP frame buer with a
capacity of b Several b and b transformations
were developed to test the dynamic range light intensity and focus settings
The eective dynamic range of an imaging system depends upon the scaling
factors used for correcting lightintensity variations and the SNR of the imag
ing system Kuduvalli et al analyzed the intensity proles of the Plannar
light box and observed that a scaling factor of about was required
FIGURE
Contour coding applied to the image in Figure Three contours are
shown for the sake of illustration of the procedure The pixels included in
the contours are shown in bold italics The initial point of each contour is
underlined Doubleheaded arrows represent two separate moves in the two
directions of the arrows
to correct the values of the pixels at the edges of the image Furthermore
the local standard deviation of the intensity levels measured by the camera
with respect to a movingaverage window of counts was estimated to be
about counts The eective noise level at the edge of the imaging eld
was estimated at counts For these reasons it was concluded that two of
the leastsignicant bits in the b data would only contribute to noise and
aect the performance of data compression algorithms In addition by
scanning a standard calibrated grayscale step pattern the eective dynamic
range of the digitizer was observed to be about OD In consideration of
all of the above factors it was determined that truncating the b pixel val
ues to b values would be adequate for representing the digitized image
This procedure reduced the eective noise level by a factor of four to about
counts
Kuduvalli et al also estimated the MTF of the digitization system from
measurements of the ESF see Sections and with a view to demon
strate the resolving capability of the system to capture submillimeter details
on Xray lms such as microcalcications on mammograms The normalized
value of the MTF at onehalf of the sampling frequency was estimated to be
which is considered to be adequate for resolving objects and details at
the same frequency which is the highest frequency component retained in the
digitized image
The average number of bits per pixel obtained for ten Xray images using the
Human arithmetic and LZW coding techniques are listed in Table the
zerothorder entropy values of the images are also listed The high values of the
zerothorder entropy indicate limits on the performance of the Human coding
technique The arithmetic coding method has given bit rates comparable
with those provided by the Human code but at considerably higher level of
complexity Figure shows plots of the average bit rate as a function of
the buer length in LZW coding for four of the ten images listed in Table
The maximum length of the symbol strings scanned was xed at L
s
The LZW code provided the best compression rates among the three methods
considered to the extent of about of the initial number of bits per pixel
in the images The average bit rate provided by LZW coding is well below the
zerothorder entropy values of the images indicating that ecient encoding
of strings of pixel values can exploit the redundancy and correlation present
in the data without performing explicit decorrelation
The Need for Decorrelation
The results of direct source encoding discussed in Section indicate the
limitations of direct encoding of the symbols generated by an image source
TABLE
Average Number of Bits Per Pixel with Direct Source Coding for Data
Compression of Digitized Mammograms and Chest Xray
Images
Image Type Size pixels Entropy Human Arith LZW
Mammo
Mammo
Mammo
Mammo
Chest
Chest
Chest
Chest
Mammo
Chest
Average
The entropy listed is the zerothorder entropy Pixel values in the original
images were quantized at bpixel See also Table Note Arith
Arithmetic coding LZW LempelZivWelch coding Mammo Mam
mogram Reproduced with permission from GR Kuduvalli and RM Ran
gayyan Performance analysis of reversible image compression techniques for
highresolution digital teleradiology IEEE Transactions on Medical Imaging
c
IEEE
FIGURE
Average bit rate as a function of the buer length using LempelZivWelch
coding for four of the ten images number and listed in Table
The abscissa indicates the value of B with the buer length given by n
B
The maximum length of the symbol strings scanned was xed at L
s
Re
produced with permission from GR Kuduvalli and RM Rangayyan Perfor
mance analysis of reversible image compression techniques for highresolution
digital teleradiology IEEE Transactions on Medical Imaging
c
IEEE
Although some of the methods described such as the LempelZiv and run
length coding methods have the potential to exploit the redundancy present
in images their eciency in this task is limited
The term decorrelation indicates a procedure that can remove or reduce
the redundancy or correlation present between the elements of a data stream
such as the pixels in an image The most commonly used decorrelation tech
niques are the following
dierentiation which can remove the commonality present between ad
jacent elements
transformation to another domain where the energy of the image is con
ned to a narrow range such as the Fourier KarhunenLo eve discrete
cosine or WalshHadamard orthogonal transform domains
modelbased prediction where the error of prediction would have re
duced information content and
interpolation where a subsampled image is transmitted the pixels in
between the preceding data are obtained by interpolation and the error
of interpolation which has reduced information content is transmitted
Observe that the decorrelated data transform coecients prediction error
etc need to be encoded and transmitted the techniques described in Sec
tion for direct source encoding may also be applied to decorrelated data
In addition further information regarding initial values and the procedures
for the management of the transform coecients or the model parameters will
also have to be sent to facilitate complete reconstruction of the original image
The advantages of decorrelating image data by dierentiation are demon
strated by the following simple example The use of transforms for image
data compression is discussed in Section Interpolative coding is brie!y
described in Section Methods for predictionbased data compression are
described in Section Techniques based upon dierent scanning strategies
to improve the performance of decorrelation by dierentiation are discussed in
Sections and Strategies for combining several decorrelation steps
are discussed in Section
Example Consider the image in Figure a the histogram of the
image is shown in part b of the gure The image has a good spread of
gray levels over its spatial extent and the histogram while not uniform does
exhibit a good spread over the dynamic range of the image The zerothorder
entropy at b is close to the maximum possible value of b for the image
with gray levels These characteristics suggest limited potential for direct
encoding methods
The image in Figure a was subjected to a simple rstorder partial
dierentiation procedure given by
f
mn fmn fm n
The result shown in Figure a has an extremely limited range of de
tails the histogram of the image shown in part b of the gure indicates
that although the image has more gray levels than the original image most of
the gray levels occur with negligible probability The concentrated histogram
leads to a lower value of entropy at b Observe that the histogram
of the dierence image is close to a Laplacian PDF see Section and
Figure The simple operation of dierentiation has reduced the entropy
of the image by about The reduced entropy suggests that the coding
requirement may be reduced signicantly Observe that the additional infor
mation required to recover the original image from its derivative as above is
just the rst row of pixels in the original image Data compression techniques
based upon dierentiation are referred to as dierential pulse code modulation
DPCM techniques DPCM techniques vary in terms of the reference value
used for subtraction in the dierentiation process The reference value may
be derived as a combination of a few neighboring pixels in which case the
method approaches linear prediction in concept
Transform Coding
The main premise of transformdomain coding of images is that when orthog
onal transforms are used the related coecients represent elements that are
mutually uncorrelated See Section for the basics of orthogonal trans
forms Furthermore because most natural images have limitations on the
rate of change of their elemental values that is they are generally smooth
their energy is conned to a narrow lowfrequency range in the transform do
main These properties lead to two characteristics of orthogonal transforms
that are of relevance and importance in data compression
orthogonal transforms perform decorrelation and
orthogonal transforms compress the energy of the given image into a
narrow region
The second property listed above is commonly referred to as energy com
paction
Example The logmagnitude Fourier spectrum of the image in Figure
a is shown in Figure a It is evident that most of the energy
of the image is concentrated in a small number of DFT coecients around
the origin in the D Fourier plane at the center of the spectrum displayed
In order to demonstrate the energycompacting nature of the DFT the cu
mulative percentage of the total energy of the image present at the
coordinate of the DFT and contained within concentric square regions of
a
b
FIGURE
a A test image of size pixels with gray levels b Graylevel
histogram of the test image Dynamic range Zerothorder entropy
b
a
b
FIGURE
a Result of dierentiation of the test image in Figure a b Gray
level histogram of the image in a Dynamic range Zerothorder
entropy b
halfwidth DFT coecients were computed The result is plot
ted in Figure b which shows that of the energy of the image is
present in the DC component and of the energy is contained within the
central DFT components around the DC point only of the energy
lies in the highfrequency region beyond the central region of the
DFT array Regardless of the small fraction of the total energy
present at higher frequencies it should be observed that highfrequency com
ponents bear important information related to the edges and sharpness of the
image see Sections and
The DFT is the most commonly used orthogonal transform in the analysis
of systems signals and images However due to the complex nature of the
basis functions the DFT has high computational requirements In spite of
its symmetry the DFT could lead to increased direct coding requirements
due to the need for large numbers of bits for quantization of the transform
coecients Regardless the discrete nature of most images at the outset could
be used to advantage in lossless recovery from transform coecients that have
been quantized to low levels of accuracy see Section
We have already studied the DFT Sections and and the WHT
Section The WHT has a major computational advantage due to the
fact that its basis functions are composed of only and and has been ad
vantageously applied in data compression In the following sections we shall
study two other transforms that are popular and relevant to data compres
sion the discrete cosine transform DCT and the KarhunenLo eve transform
KLT
The discrete cosine transform
The DCT is a modication of the DFT that overcomes the eects of discon
tinuities at the edges of the latter and is dened as
F k l
ak l
N
N
X
m
N
X
n
fmn cos
h
m
N
k
i
cos
h
n
N
l
i
for k N and l N where
ak l
if k l
otherwise
The inverse transformation is given by
F k l
N
N
X
k
N
X
l
ak lF k l cos
h
m
N
k
i
cos
h
n
N
l
i
for m N and n N The basis vectors
of the DCT closely approximate the eigenvectors of a Toeplitz matrix see
a
b
FIGURE
a Logmagnitude Fourier spectrum of the image in Figure a com
puted over a DFT array b Distribution of the total image energy
The values represent the cumulative percentage of the energy of the im
age present at the or DC position contained within square boxes of
halfwidth pixels centered in the DFT array and in the entire
DFT array The numbers of DFT coecients corresponding to the
regions are and
Section whose elements can be expressed by increasing powers of a
constant The ACF matrices of most natural images can be closely
modeled by such a matrix An N N DCT may be computed from the
results of a N N DFT thus FFT algorithms may be used for ecient
computation of the DCT
The KarhunenLoeve transform
The KLT also known as the principalcomponent transform the Hotelling
transform or the eigenvector transform is a datadependent transform
based upon the statistical properties of the given image The image is treated
as a vector f that is a realization of an imagegenerating stochastic process
If the image is of size M N the vector is of size P MN see Section
The image f may be represented without error by a deterministic linear
transformation of the form
f A g
P
X
m
g
m
A
m
A A
A
A
P
where jAj and A
m
are row vectors that make up the P P matrix
A The matrix A needs to be formulated such that the vector g leads to an
ecient representation of the original image f
The matrix A may be considered to be made up of P linearly indepen
dent row vectors that span the P dimensional space containing f Let A be
orthonormal that is
A
T
m
A
n
m n
m n
It follows that
A
T
A I or A
A
T
Then the row vectors of A may be considered to form the set of orthonormal
basis vectors of a linear transformation This formulation leads also to the
inverse relationship
g A
T
f
P
X
m
A
T
m
f
m
In the procedure described above each component of g contributes to the
representation of f Given the formulation of A as a reversible linear trans
formation g provides a complete lossless representation of f if all of its
P MN elements are made available
Suppose that in the interest of ecient representation of images via the
extraction of the most signicant information contained we wish to use only
Q P components of g The omitted components of g may be replaced
with other values b
m
m Q P Then we have an approximate
representation of f given as
"
f
Q
X
m
g
m
A
m
P
X
mQ
b
m
A
m
The error in the approximate representation as above is
f
"
f
P
X
mQ
g
m
b
m
A
m
The MSE is given by
E
T
E
P
P
mQ
P
P
nQ
g
m
b
m
g
n
b
n
A
T
m
A
n
P
P
mQ
Eg
m
b
m
The last step above follows from the orthonormality of A
Taking the derivative of the MSE with respect to b
m
and setting the result
to zero we get
b
m
Eg
m
b
m
The optimal MMSE choice for b
m
is therefore given by
b
m
Eg
m
g
m
A
T
m
Ef m Q P
that is the omitted components are replaced by their means The MMSE is
given by
min
P
P
mQ
Eg
m
g
m
P
P
mQ
EA
T
m
f f f f
T
A
m
P
P
mQ
A
T
m
f
A
m
where
f
is the covariance matrix of f
Now if the basis vectors A
m
are selected as the eigenvectors of
f
that is
f
A
m
m
A
m
and
m
A
T
m
f
A
m
because A
T
m
A
m
where
m
are the corresponding eigenvalues then we
have
min
P
X
mQ
m
Therefore the MSE may be minimized by ordering the eigenvectors the rows
of A such that the corresponding eigenvalues are arranged in decreasing
order that is
P
Then if a component g
m
of g is replaced
by b
m
g
m
the MSE increases by
m
By replacing the components of g
corresponding to the eigenvalues at the lower end of the list the MSE is kept
at its lowestpossible value for a chosen number of components Q
From the above formulation and properties it follows that the components
of g are mutually uncorrelated
g
A
T
f
A
P
where is a diagonal matrix with the eigenvalues
m
placed along its diago
nal Because the eigenvalues
m
are equal to the variances of g
m
a selection
of the larger eigenvalues implies the selection of the transform components
with the higher variance or information content across the ensemble of the
images considered
The KLT has major applications in principalcomponent analysis PCA
image coding data compression and feature extraction for pattern classi
cation Diculties could exist in the computation of the eigenvectors and
eigenvalues of the large covariance matrices of even reasonably sized images
It should be noted that a KLT is optimal only for the images represented
by the statistical parameters used to derive the transformation New trans
formations will need to be derived if changes occur in the statistics of the
imagegenerating process being considered or if images of dierent statistical
characteristics need to be analyzed
Because the KLT is a datadependent transform the transformation vectors
the matrix A need to be transmitted however if a large number of images
generated by the same underlying process are to be transmitted the same
optimal transform is applicable and the transformation vectors need to be
transmitted only once Note that the error between the original image and
the image reconstituted from the KLT components needs to be transmitted
in order to facilitate lossless recovery of the image
See Section for a discussion on the application of the KLT for the
selection of the principal components in multiscale directional ltering
Encoding of transform coecients
Regardless of the transform used such as the DFT DCT or KLT the trans
form coecients form a set of continuous random variables and have to be
quantized for encoding This introduces quantization errors in the transform
coecients that are transmitted and hence errors arise in the image recon
structed from the transformcoded image In the following paragraphs a rela
tionship is derived between the quantization error in the transform domain and
the error in the reconstructed image Kuduvalli and Rangayyan
used such a relationship to develop a method for errorfree transform coding
of images
Consider the general D linear transformation with the forward and inverse
transform kernels consisting of orthogonal basis functions amn k l and
bmn k l respectively such that the forward and inverse transforms are
given by
F k l
N
X
m
N
X
n
fmn amn k l
k l N and
fmn
N
X
k
N
X
l
F k l bmn k l
mn N See Section Now let the transform coecient
F k l be quantized to
"
F k l with a quantization error q
F
k l such that
F k l
"
F k l q
F
k l
The reconstructed image from the quantized transform coecients is given by
"
fmn
N
X
k
N
X
l
"
F k l bmn k l
The error in the reconstructed image is
q
f
mn fmn
"
fmn
N
X
k
N
X
l
F k l
"
F k l bmn k l
N
X
k
N
X
l
q
F
k l bmn k l
The sum of the squared errors in the reconstructed image is
Q
f
N
X
m
N
X
n
jq
f
mnj
N
X
m
N
X
n
N
X
k
N
X
l
q
F
k l bmn k l
N
X
k
N
X
l
q
F
k
l
bmn k
l
N
X
k
N
X
l
N
X
k
N
X
l
q
F
k l q
F
k
l
N
X
m
N
X
n
bmn k l b
mn k
l
N
X
k
N
X
l
jq
F
k lj
Q
F
where the last line follows from the orthogonality of the basis functions bmn
k l and Q
F
is the sum of the squared quantization errors in the transform
domain This result is related to Parseval s theorem in D see Equation
Applying the expectation operator to the rst and the last expressions above
we get
N
X
m
N
X
n
Ejq
f
mnj
N
X
k
N
X
l
Ejq
F
k lj
or
N
X
m
N
X
n
q
f
mn
N
X
k
N
X
l
q
F
k l Q
where the symbol indicates the expected average values of the correspond
ing variables and Q
is the expected total squared error of quantization in
either the image domain or the transform domain
It is possible to derive a condition for the minimum average number of bits
required for encoding the transform coecients for a given total distortion in
the image domain Let us assume that the transform coecients are normally
distributed If the variance of the transform coecient F k l is
F
k l the
average number of bits required to encode the coecient F k l with the MSE
q
F
k l is given by its ratedistortion function
Rk l
log
F
k l
q
F
k l
The overall average number of bits required to encode the transform coe
cients with a total squared error Q
is
R
av
N
N
X
k
N
X
l
Rk l
N
N
X
k
N
X
l
log
F
k l
q
F
k l
We now need to minimize R
av
subject to the condition given by Equa
tion Using the method of Lagrange s multiplier the minimum occurs
when
q
F
kl
n
N
P
N
k
P
N
l
log
h
F
kl
q
F
kl
i
h
Q
P
N
k
P
N
l
q
F
k l
io
k l N where is the Lagrange multiplier It follows that
N
ln q
F
k l
or
q
F
k l
N
ln
q
k l N where q
is the average MSE which is a constant for
all of the transform coecients Thus the average number of bits required
to encode the transform coecients R
av
is minimum when the total squared
error is equally distributed among all of the transform coecients
Maximum error limited encoding of transform coecients Kudu
valli and Rangayyan derived the following condition for en
coding transform coecients subject to a maximum error limit Consider a
uniform quantizer with a quantization step size S for encoding the transform
coecients such that the maximum quantization error is limited to
S
It
may be assumed that the quantization error is uniformly distributed over the
set of transform coecients in the range
S
S
see Figure Then
the average squared error in the transform domain is
q
S
From the result in Equation it is seen that the errors in the recon
structed image will also have a variance equal to q
We now wish to estimate
the fraction of the total number of pixels in the reconstructed image that are
in error by more than S This is given by the area under the tail of the PDF
of the reconstruction error shown in Figure The worst case occurs
when the entropy of the reconstruction errors is at its maximum under the
constraint that the variance of the reconstruction errors is bounded by q
this occurs when the error is normally distributed Therefore the upper
bound on the estimated fraction of the pixels in error by more than S is
ES
Z
S
q
q
exp
x
q
dx erfc
p
FIGURE
Schematic PDFs of the transformcoecient quantization error uniform PDF
solid line and the image reconstruction error Gaussian PDF dashed line
The gure represents the case with the quantization step size S
where erfcx is the error function integral dened as
erfcx
Z
x
p
exp
x
dx
Thus only a negligible number of pixels in the reconstructed image will
be in error by more than the quantization step size S Conversely if the
maximum error is desired to be S a quantization step size of S could be
used to encode the transform coecients with only a negligibly small number
of the reconstructed pixels exceeding the error limit The pixels in error by
more than the specied maximum could be encoded separately with a small
overhead When the maximum allowed error is errorfree reconstruction
of the image is possible by simply rounding o the reconstructed values to
integers
Variable length encoding and bit allocation The lower bound on
the average number of bits required for encoding normally distributed trans
form coecients F k l with the MSE q
F
k l is given by Equation
Goodnessoft studies of PDFs of transform coecients have shown
that transform coecients tend to follow the Laplacian PDF The PDF of a
transform coecient F k l may be modeled as
pF k l
k l
expk l jF k lj
where k l
p
F
k l is the constant parameter of the Laplacian PDF
A shift encoder could be used to encode the transform coecients such that
the maximum quantization error is S The shiftencoding procedure is
shown in Figure In a shift encoder k l levels are nominally allo
cated to encode a transform coecient F k l covering the range fk l
gS fk lgS with the codes k l The code k l
indicates that the coecient is out of the range fk l gS fk l
gS For the outofrange coecients an additional k l levels are allo
cated to cover the ranges fk lgSk lS and k lS fk l
gS The process is repeated with the allocation of additional levels until
the actual value of the transform coecient to be encoded is reached If the
code value is represented by a simple binary code at each level the average
number of bits required to encode the transform coecient F k l is given
by
Rk l
log
k l
expk l k l S
and the nominal number of bits allocated to encode the transform coecient
is
bk l dlog
k le
It is now required to nd the bk l that minimizes the average number of
bits Rk l required to encode F k l This can be done by using a nonlinear
optimization technique such as the NewtonRaphson method However be
cause only integral values of bk l need to be searched it is computationally
less expensive to search the space of bk l for the corresponding
minimum value of Rk l which is the range of the nominal number of bits
allocated to encode the transform coecient F k l
FIGURE
Schematic representation of shift coding of transform coecients with a Lapla
cian PDF With reference to the discussion in the text the gure represents
k lS and k l
This allocation requires an estimate of the variance
F
k l of the transform
coecients or a model of the energycompacting property of the transform
used Most of the linear orthogonal transforms used in practice result in the
concentration of the variance in the lowerorder transform coecients The
variance of the transform coecients may be modeled as
F
k l F exp
k
l
where F is the lowestorder transform coecient and and are the
constants of the model For most transforms except the KLT F is the
average value or a scaled version of the average value of the image pixels The
parameters and may be estimated using a leastsquares t to the rst few
transform coecients of the image
In an alternative coding procedure a xed total number of bits may be
allocated for encoding the transform coecients and the dierence image
encoded by a lossless encoding method In such a procedure no attempt is
made to allocate additional bits for transform coecients that result in errors
that fall out of the quantization range The total number of bits allocated to
encode the transform coecients is varied until the total average bit rate is at
its minimum Using such a procedure Cox et al found an optimal com
bination of bit rates between the error image and the transform coecients
Wang and Goldberg used a method of requantization of the quantization
errors in addition to encoding the error image they too observed a minimum
total average bit rate after a number of iterations of requantization Exper
iments conducted by Kuduvalli and Rangayyan showed that
the lowest bit rates obtained by such methods for reversible compression can
also be obtained by allocating additional bits to quantize the outofrange
transform coecients as described earlier This is due to the fact that a large
quantization error in the transform domain while needing only a few addi
tional bits for encoding will get redistributed over a large number of pixels
in the image domain thereby increasing the entropy of the error image
The large sizes of image arrays used for highresolution representation of
medical images preclude the use of fullframe transforms Partitioning an im
age into blocks not only leads to computational advantages but also permits
adaptation to the changing statistics of the image In the coding procedure
used by Kuduvalli and Rangayyan the images listed in Ta
ble were partitioned into blocks of size pixels The model
parameters and in Equation were computed for each block by us
ing a leastsquares t to the corresponding set of transform coecients The
parameters were stored or transmitted along with the encoded transform co
ecients in order to allow the decoder to reconstruct the model and the bit
allocation table Blocks at the boundaries of the images that were not squares
were encoded with D DPCM and the Human code
Figure shows the average bit rate for one of the images listed in
Table as a function of the maximum allowed error using four transforms
The KLT with the ACF estimated from the image and the the DCT show
the best performance among the four transforms The performance of the
KLT is only slightly superior to and in some cases slightly worse than that
of the DCT this is to be expected because of the general nonstationarity of
medical images and due to the problems associated with the estimation of
the ACF matrix from a nite image
Figure shows the average bit rate obtained by using the DCT as a
function of the maximum allowed error for four of the ten images listed in Ta
ble When the maximum allowed error is errorfree reconstruction
of the original image is possible otherwise the compression is irreversible
The average bit rate for errorfree reconstruction is seen to be in the range of
bpixel for the images considered down from bpixel in their original
format
FIGURE
Average bit rate as a function of the maximum allowed error using four trans
forms KLT DCT DFT and WHT with the rst image listed in Table
A block size of pixels was used for each transform The compression
is lossless if the maximum allowed error is otherwise it is irreversible or
lossy See also Table Reproduced with permission from GR Kuduvalli
and RM Rangayyan Performance analysis of reversible image compression
techniques for highresolution digital teleradiology IEEE Transactions on
Medical Imaging
c
IEEE
FIGURE
Average bit rate as a function of the maximum allowed error using the DCT
for the rst four images listed in Table A block size of pixels
was used The compression is lossless if the maximum allowed error is
otherwise it is irreversible or lossy See also Table Reproduced
with permission from GR Kuduvalli and RM Rangayyan Performance
analysis of reversible image compression techniques for highresolution digital
teleradiology IEEE Transactions on Medical Imaging
c
IEEE
Interpolative Coding
Interpolative coding consists of encoding a subsampled image using a re
versible compression technique deriving the values of the remaining pixels
via interpolation with respect to their neighboring pixels that have already
been processed and then encoding the dierence between the actual pixels
and the interpolated pixels in successive stages using discrete symbol coding
techniques This technique is also referred to as hierarchical interpolation
HINT and is illustrated in Figure In the image shown in
the gure the pixels marked correspond to the original image decimated
by a factor of The decimated image could be encoded using any coding
technique The pixels marked are estimated from those marked by
bilinear interpolation and rounded to ensure reversibility This completes one
iteration of interpolation Next the pixels marked are interpolated from
the pixels marked and and the process is repeated to interpolate the
pixels marked and The dierences between the actual pixel values
marked and the corresponding interpolated values form a discrete
symbol set with a small dynamic range and may be encoded eciently using
Human arithmetic or LZW coding
The illustration in Figure corresponds to interpolation of order four
higherorder interpolation may also be used In general interpolative coding
of order
P
where P is an integer involves P iterations of interpolation
In the work of Kuduvalli and Rangayyan the digitized radio
graphic images listed in Table were partitioned into blocks of size
pixels for interpolative coding The initial subsampled images were decorre
lated using D DPCM of order the dierence data were compressed
by Human coding It was observed that the dierences between the inter
polated and actual pixel values could be modeled by Laplacian PDFs see
Figure The variance of the interpolation errors was seen to decrease
with increasing resolution because pixels are correlated more to their immedi
ate neighbors than to pixels farther away The interpolation errors at dierent
iterations were modeled by Laplacian PDFs with the variance equal to the cor
responding meansquared interpolation errors The PDFs were then used in
compressing the interpolation errors via arithmetic or Human coding LZW
coding does not need modeling of the error distribution but performed con
siderably worse than the Human and arithmetic codes Figure shows
the average bit rate for eight images as a function of the order of interpo
lation using Human coding as the postencoding technique It is observed
that increasing the order of interpolation has only a small eect on the overall
average bit rate
For examples on the performance of HINT see Tables
and
FIGURE
Stages of interpolative coding The pixels marked are coded and trans
mitted rst Next the pixels marked are estimated from those marked
and the dierences are transmitted The procedure continues iteratively
with the pixels marked and Reproduced with permission from
GR Kuduvalli and RM Rangayyan Performance analysis of reversible im
age compression techniques for highresolution digital teleradiology IEEE
Transactions on Medical Imaging
c
IEEE
FIGURE
Results of compression of eight of the images listed in Table by inter
polative and Human coding Order zero corresponds to D DPCM
coding Reproduced with permission from GR Kuduvalli and RM Ran
gayyan Performance analysis of reversible image compression techniques for
highresolution digital teleradiology IEEE Transactions on Medical Imaging
c
IEEE
Predictive Coding
Samples of reallife signals and images bear a high degree of correlation espe
cially over small intervals of time or space The correlation between sam
ples may also be viewed as statistical redundancy An outcome of these
characteristics is that a sample of a temporal signal may be predicted from
a small number of its preceding samples When such prediction is per
formed in a linear manner we have a linear prediction LP model given
by
"
fn
P
X
p
ap fn p Gdn
where fn is the signal being modeled
"
fn is the predicted value of fn
ap p P are the coecients of the LP model and P is the order
of the model The signal fn is considered to be the output of a linear
system or lter with dn as the input or driving signal G is the gain of
the system Because the prediction model uses only the past values of the
output signal fn it is known as an autoregressive AR model The need for
causality in physically realizable lters and signal processing systems dictates
the requirement to use only the past samples of f and the present value of
d in predicting the current value fn If the past values of d are also used
the model will include a movingaverage or MA component
The model represented in Equation indicates that given the initial
set of P values of the signal f and the input or driving signal d any future
value of f may be approximately computed with a knowledge of the set
of coecients ap Therefore the model coecients ap p P
represent the signalgenerating process The model coecients may be used
to predict the values of f or to analyze the signalgenerating process Several
methods exist to derive the LP model coecients for a given signal subject
to certain conditions on the error of prediction
In the context of image data compression we have a few considerations that
dier from the temporal signal application described above In most cases the
input or driving signal d is not known the omission of the related component
in the model represented in Equation will cause only a small change in
the error of prediction Furthermore causality is not a matter of concern in
image processing however a certain sequence of accessing or processing of
the pixels in the given image needs to be dened which could imply past
and future samples of the image see Figure In the context of image
processing the samples used to predict the current pixel could be labeled as
the ROS of the model Then we could express the basic D LP model as
"
fmn
X X
pqROS
ap q fm p n q
In the application to image coding the ROS needs to be dened such that
in the decoding process only those pixels that have already been decoded are
included in the ROS for the current pixel being processed
The error of prediction is given by
emn fmn
"
fmn
and the MSE between the original image and the predicted image is given by
Ee
mn
The coecients ap q need to be chosen or derived so as to minimize the MSE
between the original image and the predicted image Several approaches are
available for the derivation of optimal predictor coecients
Exact reconstruction of the image requires that in addition to the initial
conditions and the model coecients the error of prediction be also trans
mitted and made available to the decoder The advantage of LPbased image
compression lies in the fact that the error image tends to have a more concen
trated PDF than the original image close to a Laplacian PDF in most cases
see Figure b which lends well to ecient data compression
In the simplest model of prediction the current pixel fmn may be mod
eled as being equal to the preceding pixel fmn or fm n let
"
fmn fm n
Then the error of prediction is given by
emn fmn fm n
which represents simple DPCM Note Any data coding method based upon
the dierence between a data sample and a predicted value of the same using
any scheme is referred to as DPCM hence all LPbased methods fall under
the category of DPCM The dierentiation procedure given by Equation
and illustrated in Figure is equivalent to LP with
"
fmn fm n
as above Several simple combinations of the immediate neighbors of the
current pixel may also be used in the prediction model see Section
Ecient modeling requires the use of the optimal model order ROS size
and the derivation of the optimal coecients subject to conditions related to
the minimization of the MSE Several methods for the derivation of the model
coecients are described in the following sections
Twodimensional linear prediction
The error of prediction in the D LP model is given by
emn fmn
"
fmn
fmn
X X
pqROS
ap q fm p n q
The squared error is given by
e
mn f
mn
X X
pqROS
ap q fmn fm p n q
X X
pqROS
X X
rsROS
ap q ar s fm p n q fm r n s
Applying the statistical expectation operator we get
Ee
mn
f
X X
pqROS
ap q
f
p q
X X
pqROS
X X
rsROS
ap q ar s
f
r p s q
where
f
p q is the ACF of f and the imagegenerating process is assumed
to be widesense stationary
The coecients that minimize the MSE may be derived by setting to zero
the derivative of
with respect to ap q p q ROS which leads to
f
r s
X X
pqROS
ap q
f
r p s q r s ROS
Using this result in Equation we get
f
X X
pqROS
ap q
f
p q
Combining Equation and we get
f
r s
X X
pqROS
ap q
f
r p s q
r s ROS
r s
The equations represented by the expressions above are known as the D
normal or YuleWalker equations and may be solved to derive the prediction
coecients Because the method uses the ACF of the image to derive the pre
diction coecients it is known as the autocorrelation method The ACF
may be estimated from the given nite image fmn m M
n N as
"
f
p q
MN
M
X
mp
N
X
nq
fmn fm p n q
The prediction coecients may also be estimated by using leastsquares
methods to minimize the prediction error averaged over the entire image
indicated as mn IMG as follows
MN
X X
mnIMG
e
mn
MN
X X
mnIMG
f
mn
X X
pqROS
ap q fmn fm p n q
X X
pqROS
X X
rsROS
ap q ar s fm p n q fm r n s
f
X X
pqROS
ap q
f
p q
X X
pqROS
X X
rsROS
ap q ar s
f
p q r s
Here
f
is the covariance of the image f dened as
f
p q r s
MN
X X
mnIMG
fm p n q fm r n s
The coecients that minimize the averaged error may be derived by setting
to zero the derivative of
with respect to ap q which leads to
f
r s
X X
pqROS
ap q
f
p q r s
It follows that
f
X X
pqROS
ap q
f
p q
The D normal equations for this condition are given by
f
r s
X X
pqROS
ap q
f
p q r s
r s ROS
r s
which may be solved to obtain the prediction coecients Because the covari
ance of the image is used to derive the prediction coecients this method is
known as the covariance method
Now if the region IMG is dened so as to span the entire image array
of size M N and the image is assumed to be zero outside the range m
M n N we have
f
p q r s
MN
M
X
m
N
X
n
fm p n q fm r n s
"
f
r p s q
where
"
f
is an estimate of the ACF of the image f Then the covariance
method yields results that are identical to the results given by the autocorre
lation method
Equation may be expressed in matrix form as
f
a
where for the case of a QP ROS of size P
P
the matrices
f
a and
are of size P
P
P
P
P
P
and
P
P
respectively The extended ACF matrix
f
is given by
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
P
where the subscript f has been dropped from the entries within the matrix
for the sake of brevity The matrices composed by the prediction coecients
and the error are given by
a
a
a P
a
a P
a
aP
P
and
The matrix
f
is ToeplitzblockToeplitz in nature ecient algorithms are
available for the inversion of such matrices
The methods described above to compute the prediction coecients assume
the image to be stationary over the entire frame available In practice most
images are nonstationary and may be assumed to be locally stationary only
over relatively small segments or ROIs In order to maintain the optimality
of the model over the entire image and in order to maintain the error of
prediction at low levels we could follow one of two procedures
partition the image into small blocks over which stationarity may be
assumed and compute the prediction coecients independently for each
block or
adapt the model to the changing statistics of the image
Encoding the prediction error In order to facilitate errorfree recon
struction of the image the prediction error has to be transmitted and made
available at the decoder in addition to the prediction coecients and the ini
tial conditions For quantized original pixel values the prediction error may
be rounded o to integers that may be encoded without error using a source
coding technique such as the Human code The prediction error has been
observed to possess a Laplacian PDF which lends well to ecient
compression by the Human code
Results of application to medical images The average bit rates ob
tained by the application of the autocorrelation method of computing the
prediction coecients using blocks of size pixels to six of the im
ages listed in Table are shown in Figure for various model orders
For the images used comparable performance was obtained using NSHP and
QP ROSs of similar extent Good compression performance was obtained
with average bit rates in the range bpixel for most of the images
with the original pixel values at bpixel using QP ROSs of size or
See Section for a method for the detection of calcications based
upon the error of prediction
Multichannel linear prediction
The dierence between LP in D and D may be bridged by multichannel
LP where a certain number of rows of the given image may be viewed as a
collection of multichannel signals Kudu
valli and Rangayyan proposed the following procedures for the
application of multichannel LP to predictive coding and compression of D
images
Consider a multichannel signal with P
channels with the channels
indexed as q P
and the individual signals labeled as fm
m N The collection of the signal s values at a position m
given by f
q
m q P
may be viewed as a multichannel signal or
vector or a matrix of size P
see Figures and If we
were to use a multichannel linear predictor of order P
we could predict the
vector fm as a linear combination of the vectors fm p p P
"
fm
P
X
p
ap fm p
where ap p P
are multichannel LP coecient matrices each of
size P
P
FIGURE
Results of compression of six of the images listed in Table by D LP
autocorrelation method and Human coding The method was applied on
a blockbyblock basis using blocks of size pixels In the case of
modeling using the NSHP ROS a model order of indicates a
ROS see Figure indicates a ROS etc The orders of
models using the QP ROS are indicated by integers indicates a ROS
indicates a ROS etc Reproduced with permission from GR Kuduvalli
and RM Rangayyan Performance analysis of reversible image compression
techniques for highresolution digital teleradiology IEEE Transactions on
Medical Imaging
c
IEEE
FIGURE
Multichannel linear prediction Each row of the image is viewed as a channel
or component of a multichannel signal or vector The column index of the
image may be considered to be equivalent to a temporal index
The indices shown correspond to Equations and See also
Figure
The error of prediction is given by
em fm
P
X
p
ap fm p
The covariance matrix of the error of prediction is given by
e
Eem e
T
m
For optimal prediction we need to derive the prediction coecient matrices
ap that minimize the trace of the covariance matrix of the error of prediction
From Equation we can write
em e
T
m fm f
T
m
P
X
p
ap fm p f
T
m
P
X
q
fm f
T
m q a
T
q
P
X
p
P
X
q
ap fm p f
T
m q a
T
q
Applying the statistical expectation operator and assuming widesense sta
tionarity of the multichannel signalgenerating process we get
e
c
P
X
p
ap
c
p
FIGURE
Multichannel LP applied to a D image See also Figure
Reproduced with permission from GR Kuduvalli and RM Rangayyan An
algorithm for direct computation of D linear prediction coecients IEEE
Transactions on Signal Processing
c
IEEE
P
X
q
c
q a
T
q
P
X
p
P
X
q
ap
c
q p a
T
q
where
c
is the ACF of the image computed over the set of rows or channels
being used in the multichannel prediction model given by
c
r
r
P
r
P
r
P
P
r
and
pq
r Ef
p
m f
q
m r
In order to minimize the trace of the error covariance matrix
e
we could
dierentiate both the sides of Equation with respect to the prediction
coecient matrices ar r P
and equate the result to the null
matrix of size P
P
which leads to
e
ar
r P
c
r
P
X
p
T
c
r p a
T
p
c
r
P
X
p
c
p r a
T
p r P
Note
T
c
r p
c
p r
Now Equation may be rewritten as
e
c
P
X
p
ap
c
p
P
X
q
c
q
P
X
p
ap
c
q p
a
T
q
c
P
X
p
ap
T
c
p
c
P
X
p
c
p a
T
p
The relationships derived above may be summarized as
c
A
where the matrices are given in expanded form as
c
c
c
P
c
c
c
P
c
P
c
P
c
I
a
T
a
T
P
e
In the equation above the submatrices
c
are as dened in Equation
I is the identity matrix of size P
P
and is the null matrix
of size P
P
This system of equations may be referred to as
the multichannel version of the YuleWalker equations The solution to this
set of equations may be used to obtain the D LP coecients by making the
following associations
pq
r
f
r p q
compare Equations and
e
a
and
a
r
a
T
r a
where a
r
ar ar ar P
T
is composed by the elements of the
r
th
row of the matrix of prediction coecients a given in Equation
written as a column matrix
The Levinson Wiggins Robinson algorithm The multichannel pre
diction coecient matrix may be obtained by the application of the algorithms
due to Levinson and Wiggins and Robinson In the multichannel
version of this algorithm the prediction coecients for order P
are recursively related to those for order P
The prediction model given by
Equation is known as the forward predictor Going in the opposite di
rection the backward predictor is dened to predict the vector fm in terms
of the vectors fm p p P
as
"
fm
P
X
p
bp fm p
where bp p P
are the multichannel backward prediction coe
cient matrices each of size P
P
see Figure
In order to derive the multichannel version of Levinson s algorithm let us
rewrite the multichannel ACF matrix in Equations and as follows
P
P
P
P
P
where the subscript P
is used to indicate the order of the model and
the subscript c has been dropped from the submatrices for compact notation
The matrix
P
may be partitioned as
P
T
P
P
P
P
P
T
P
where
P
P
T
P
P
P
T
and the property that r
T
r has been used It follows that
P
P
P
and
P
P
P
Let us also dene partitions of the forward and backward prediction coecient
matrices as follows
A
P
I
"
A
P
and
B
P
"
B
P
I
where A
P
is the same as A in Equations and and B
P
is formed
in a similar manner for the backward predictor
Using the partitions as dened above we may rewrite the multichannel
YuleWalker equations given by Equation in two forms for forward
and backward prediction as follows
T
P
T
P
P
I
"
A
P
f
P
P
and
P
P
T
P
"
B
P
I
P
b
P
where
P
is the null matrix of size P
P
P
The matrix
f
P
is the covariance matrix of the error of forward prediction the same as
e
given by Equation with the matrix
b
P
being the counterpart for the
backward predictor It follows that
"
A
P
P
P
"
B
P
P
P
f
P
T
P
"
A
P
and
b
P
T
P
"
B
P
Applying the inversion theorem for partitioned matrices and making
use of the preceding six relationships we get
P
T
P
P
P
A
P
f
P
A
T
P
P
P
T
P
B
P
b
P
B
T
P
Multiplying both sides of Equation by
P
making use of the par
titioned form in Equation and using Equation we get
"
A
P
"
A
P
B
P
b
P
B
T
P
P
Extracting the lower P
P
matrix from both sides of Equa
tion in its partitioned form we have
a
T
P
P
b
P
B
T
P
P
which upon transposition yields
a
P
P
T
P
B
P
b
P
b
P
T
b
P
where
b
P
B
T
P
P
P
P
X
p
b
P
pP
p
Using Equation in Equation we get
"
A
P
"
A
P
B
P
a
T
P
P
which upon transposition and expansion of the matrix notation yields
a
P
p a
P
pa
P
P
b
P
P
p p P
Similarly multiplying both sides of Equation by
P
and using Equa
tion we get
b
P
P
f
P
T
f
P
where
f
P
P
P
X
p
a
P
pP
p
and
b
P
p b
P
pb
P
P
a
P
P
p p P
Substituting Equation in Equation and using the partitioned
form of
P
in Equation we get
f
P
T
P
"
A
P
T
P
B
P
a
T
P
P
f
P
b
P
T
a
T
P
P
f
P
a
P
P
b
P
a
T
P
P
where Equation has been used in the last step Similarly we can obtain
an expression for the covariance matrix of the backward prediction error as
b
P
b
P
b
P
P
f
P
b
T
P
P
Now consider the matrix product
B
T
P
P
A
P
B
T
P
T
P
P
P
A
P
B
T
P
P
B
T
P
P
A
P
b
P
b
P
I
"
A
P
b
P
Taking the transpose of the expression above and noting that
P
is sym
metric we get
b
P
T
A
T
P
P
B
P
f
P
Equations and consti
tute the LevinsonWigginsRobinson algorithm with the initialization
f
b
c
With the autocorrelation matrices p dened by the association given in
Equation the matrix
P
is a ToeplitzblockToeplitz matrix the
block elements submatrices p along the diagonals of
P
are mutually
identical and furthermore the elements along the diagonals of each subma
trix p are mutually identical Thus
P
and p are symmetrical
about their cross diagonals that is they are persymmetric A property of
ToeplitzblockToeplitz matrices that is of interest here is dened in terms of
the exchange matrices
J
of size P
P
and
J
P
J
J
J
of size P
P
P
P
such that JJ I and J
P
J
P
I
P
With these denitions we have
JpJ
T
p
and
J
P
P
J
P
P
Now premultiplying both sides of Equation by J
P
and postmultiplying
both sides by J we get
P
P
P
P
J a
T
P
P
J
J a
T
P
J
I
J
f
P
J
Equation is identical to the modied YuleWalker equations for com
puting the matrices b
P
p and
b
P
in Equation Comparing the terms
in the two equations we get
b
P
p J a
T
P
p J p P
and
b
P
J
f
P
J
P
With these simplications the recursive procedures in the multichannel
Levinson algorithm may be modied for the computation of D LP coecients
as follows
P
f
P
P
X
p
a
P
p
f
P
p
a
P
P
P
J
P
J
a
P
p a
P
p a
P
P
J a
P
P
p J p P
and
P
P
a
P
P
J
P
J a
T
P
P
with the initialization
f
where
f
p
f
p
f
pP
f
p P
f
p
Equations constitute the D Levinson algorithm for solving the
D YuleWalker equations for the case of a QP ROS The Levinson algorithm
provides results that are identical to those obtained by direct inversion of
f
Computation of the D LP coecients directly from the image
data The multichannel version of the Levinson algorithm may be used to
derive the multichannel version of the Burg algorithm as fol
lows Equation may be augmented using the partition shown in Equa
tion as
A
P
A
P
B
P
a
T
P
P
From Equation the forward prediction error vector for order P
is
e
f
P
m fm
P
X
p
a
P
p fm p
A
T
P
F
P
m
where F
P
m is the multichannel data matrix of order P
at m
which may be partitioned as
F
P
m
fm
F
P
m
F
P
m
fm P
Similarly the backward prediction error vector may be expressed as
e
b
P
m fm P
P
X
p
b
P
p fm P
p
B
T
P
F
P
m
Transposing both sides of Equation and multiplying by F
P
m
using the partitioned forms shown on the righthand side of Equation
as well as using Equations and we get
e
f
P
m e
f
P
m a
P
P
e
b
P
m
Similarly the backward prediction error vector is given by
e
b
P
m e
b
P
m b
P
P
e
f
P
m
The matrices a
P
P
and b
P
P
known as the re!ection
coecient matrices that minimize the the sum of the squared forward
and backward prediction errors over the entire multichannel set of data points
N given by
c
Tr
N
X
mP
n
e
f
P
m e
f
P
m
T
e
b
P
m e
b
P
m
T
o
are obtained by solving
b
P
b
P
E
b
P
a
P
P
a
P
P
f
P
E
f
P
f
P
b
P
b
P
E
fb
P
b
P
E
fb
P
f
P
where
E
f
P
N
X
mP
e
f
P
m e
f
P
m
T
E
b
P
N
X
mP
e
b
P
m e
b
P
m
T
and
E
fb
P
N
X
mP
e
f
P
m e
b
P
m
T
Equations and may be used to compute the
multichannel re!ection coecients directly from the image data without com
puting the ACF
In order to adapt the multichannel version of the Burg algorithm to the
D image case we could force the structure obtained by relating the D
and multichannel ACFs in Equation on to the expressions in Equations
and and redene the error covariance matrices E
f
P
E
b
P
and
E
fb
P
to span the entire M N image Then the D counterpart of the
re!ection coecient matrix a
P
P
is obtained by solving the following
equation
P
P
E
b
P
a
P
P
a
P
P
P
E
f
P
P
P
P
E
fb
P
P
E
fb
P
P
In order to compute the error covariance matrices E
f
P
E
b
P
and E
fb
P
a strip
of width P
is dened so as to span the top P
rows of the image
as shown in Figure and the strip is moved down one row at a time The
region over which the summations are performed includes only those parts of
the strip for which the forward and backward prediction operators do not run
out of data At the beginning of the recursive procedure the error values are
initialized to the actual values of the corresponding pixels Furthermore the
forward and backward prediction error vectors are computed by forcing the
relationship in Equation on to Equation resulting in
e
b
P
m e
b
P
m J a
P
P
J e
f
P
m
The D Burg algorithm for computing the D LP coecients directly from
the image data may be summarized as follows
The prediction error covariance matrix
is initialized to
f
The prediction error vectors are computed using Equations and
The prediction error covariance matrices E
f
P
E
b
P
and E
fb
P
are com
puted from the prediction error vectors using Equations
and summing over strips of width P
rows of the image
The re!ection coecient matrix a
P
P
is obtained from the
prediction error covariance matrices by solving Equation which
is of the form AX XB C and can be solved by using Kronecker
products
The remaining prediction coecient matrices a
P
p are computed by
using Equation and the expected value of the prediction error
covariance matrix
P
is updated using Equation
When the recursive procedure reaches the desired order P
the D LP
coecients are computed by solving Equations and
TABLE
Variables in the D LP Burg and Levinson Algorithms for LP
Variable Size Description
fmn M N D image array
fm P
Multichannel vector at column m
spanning P
rows
r P
P
Autocorrelation submatrix related
to fm
P
P
Extended D autocorrelation matrix
P
P
a P
P
D LP coecient matrix
a
P
p P
P
Multichannelequivalent prediction
coecient matrix
a
P
P
P
P
Multichannelequivalent re!ection
coecient matrix
P
P
P
Multichannel prediction error
covariance matrix
e
f
mn M N Forward prediction error array
e
b
mn M N Backward prediction error array
e
f
m P
Multichannel forward prediction
error vector
e
b
m P
Multichannel backward prediction
error vector
emn scalar n
th
element of the vector em
E
f
P
P
P
Forward prediction
error covariance matrix
E
b
P
P
P
Backward prediction
error covariance matrix
E
fb
P
P
P
Forwardbackward prediction
error covariance matrix
Bold characters represent vectors or matrices A QP ROS of size P
P
is
assumed
The variables involved in the D Burg and Levinson algorithms are summa
rized in Table
The modied multichannel version of the Burg algorithm oers advantages
similar to those of its D counterpart over the direct inversion method it is a
fast and ecient procedure to compute the prediction coecients and predic
tion errors without computing the autocorrelation function The optimization
of the prediction coecients does not make any assumptions about the im
age outside its nite dimensions and hence should result in lower prediction
errors and ecient coding Furthermore the forced D structure makes the
algorithm computationally more ecient than the direct application of the
multichannel Burg procedure
Computation of the prediction error In order to compute the predic
tion error for coding and transmission the trace of the covariance matrix in
Equation may be minimized using Equations and and
eliminating the covariance matric
P
as follows From Equation
the squared forward and backward prediction error vectors in the D Burg
algorithm are given as
E
P
N
X
mP
h
e
f
P
m e
f
P
m
T
e
b
P
m e
b
P
m
T
i
N
X
mP
h
e
f
P
m a
P
P
e
b
P
m
i
h
e
f
P
m a
P
P
e
b
P
m
i
T
h
e
b
P
m J a
P
P
J e
f
P
m
i
h
e
b
P
m J a
P
P
J e
f
P
m
i
T
E
f
P
E
b
P
E
fb
P
a
T
P
P
a
P
P
E
fb
P
T
a
P
P
E
b
P
a
T
P
P
E
fb
P
T
J a
T
P
P
J
J a
P
P
J E
fb
P
J a
P
P
J E
f
P
J a
T
P
P
J
The D Burg algorithm for the purpose of image compression consists of
determining the re!ection coecient matrix a
P
P
that minimizes the
trace of the error covariance matrix E
P
This is achieved by dierentiating
Equation with respect to a
P
P
and equating the result to the
null matrix which yields
E
fb
P
T
E
b
P
a
P
P
E
fb
P
J a
P
P
J E
f
P
or
J E
b
P
a
P
P
a
P
P
J E
f
P
J
n
E
fb
P
E
fb
P
T
o
If the D autocorrelation matrices
f
r are symmetric the matrix a
P
P
will also be symmetric which reduces Equation to
E
b
P
a
P
P
a
P
P
E
f
P
n
E
fb
P
E
fb
P
T
o
E
fb
P
When the recursive procedure reaches the desired order P
the multichannel
equivalent D prediction error image is obtained as
e
mP
p q e
f
P
mN qp
p P
q N m
M
P
Equation is in the form of the normal equations for D LP This
suggests that the D Burg algorithm for LP may be applied to the
multichannelequivalent D prediction error image to obtain the nal predic
tion error image in a recursive manner as follows
Compute the sum of the squared forward and backward prediction errors
as
f
N
X
mP
M
X
n
je
P
mnj
b
N
X
mP
M
X
n
jc
P
mnj
and
fb
N
X
mP
M
X
n
e
P
mn c
P
mn
where c
P
mn is theMN backward prediction error array initialized
as
c
mn e
mn m M n N
Compute the coecient a P
known as the re!ection coecient
as
a P
fb
f
b
Obtain the prediction errors at higher orders as
e
P
mn e
P
mn a P
c
P
mn
and
c
P
mn a P
e
P
mn e
P
mn
When the desired order P
is reached the prediction errors e
P
mn are
encoded using a method such as the Human code The re!ection coecient
matrices a
P
P
and the D re!ection coecients a P
are also
encoded and transmitted as overhead information
Error free reconstruction of the image from the forward and back
ward prediction errors In order to reconstruct the original image at
the decoder without any error the prediction coecients need to be re
computed from the re!ection coecients The prediction coecients a p
p P
may be computed recursively using the Burg algorithm as
a p a p a q a q p p q q P
The multichannelequivalent D prediction error image is given by
e
mn
P
X
p
a p e
m p n e
P
mn
n N m P
P
M
The multichannel prediction error vectors are related to the error data dened
above as
e
f
P
mN qp e
mP
p q
p P
q N m
M
P
The multichannel signal vectors may be reconstructed from the error vectors
via multichannel prediction as
"
fm
P
X
p
ap fmpe
f
P
m m P
P
N
Finally the original image is recovered from the multichannel signal vectors
as
f mP
p q fmN qp
p P
q N m
M
P
with rounding of the results to integers In a practical implementation for
values of e
P
mn exceeding a preset limit the true image pixel values
would be transmitted and made available directly at the decoder
Results of application to medical images Kuduvalli and Rangayyan
applied the D Levinson and Burg algorithms described above
to the highresolution digitized medical images listed in Table The
average bit rate with lossless compression of the test images using the D
blockwise LP method described in Section the D Levinson algorithm
and the D Burg algorithm were respectively and bpixel
with the original images having bpixel see also Table The multi
channel LP algorithms in particular the D Burg algorithm provided better
compression than the other methods described in the preceding sections in
this chapter The LP models described in this section are related to AR mod
eling for spectral estimation Kuduvalli and Rangayyan found the D
Burg algorithm to provide good D spectral estimates that were comparable
to those provided by other AR models
Adaptive D recursive leastsquares prediction
The LP model with constant prediction coecients given by Equation is
based on an inherent assumption of stationarity of the imagegenerating pro
cess The multichannelbased prediction methods described in Section
are twopass methods where an estimation of the statistical parameters of
the image is performed in the rst pass such as for example the autocor
relation matrix of the image in the D Levinson method and the parame
ters are then used to estimate the prediction coecients in the second pass
Once computed the same prediction coecients are used for prediction over
all of the image data from which the coecients were estimated However
this assumption of stationarity is rarely valid in the case of natural images
as well as biomedical images To overcome this problem in the case of the
multichannelbased methods the approach taken was that of partitioning the
image into blocks and computing the prediction coecients independently
for each block Another possible approach is to adapt the coecients recur
sively to the changing statistical characteristics of the image In this section
the basis for such adaptive algorithms is described and a D recursive least
squares D RLS algorithm for adaptive computation of the LP coecients
is formulated The procedures are based upon adaptive lter theory in
D and in multichannel signal ltering
With reference to the basic D LP model given in Equation several
approaches are available for adaptive computation of the coecients ap q
for each pixel being predicted at the location mn The approach
based on Wiener lter theory see Section leading to the D
LMS algorithm see Section although applicable to image
compression suers from the fact that the estimation of the coecients
ap q does not make use of all the image data available up to the current
location Adaptive estimation of the coecients based upon the Kalman
lter see Section where the prediction
coecients are represented as the state vector describing the current state
of the imagegenerating process has not been explored much However this
approach depends upon the statistics of the image represented in terms of
ensemble averages because only estimates of the ensemble averages can be
obtained this approach is likely to be suboptimal
The approach that is described in this section for adaptive prediction based
upon the work of Kuduvalli is founded upon the method of least squares
This approach is deterministic in its formulation and involves the minimiza
tion of a weighted sum of prediction errors In Section it was observed
that the estimation of the prediction coecients based on the direct mini
mization of the actual prediction errors the D Burg method yielded better
results in image compression than the method based on the estimation of an
ensemble image statistic the D ACF from the image data the D Levin
son method This result suggests that a deterministic approach could also
be appropriate for the adaptive computation of prediction coecients
In D RLS prediction the aim is to minimize a weighted sum of the squared
prediction errors computed up to the present location given by
mn
X X
pqROS
wmn p q ep q
where ep q is the prediction error at p q and wmn p q is a weight
ing factor chosen to selectively forget the errors from the preceding pixel
locations the past in order for the prediction coecients to adapt to the
changing statistical nature of the image at the current location Boutalis et
al used an exponential weighting factor whose magnitude reduces in the
direction opposite to the scanning model used in the generation of the image
With this weightingfactor model and special ROSs Boutalis et al used the
multichannel version of the RLS algorithm directly for adaptive estimation of
images however their weightingfactor model does not take into account the
D nature of images the weight assigned to the error at a location adjacent
in the row direction to the current location is higher than the weight assigned
to the error at a location adjacent in the column direction Kuduvalli
proposed a weightingfactor model that is truly D in its formulation In this
method using a rectangular region spanning the image up to the current loca
tion for minimizing the sum of the prediction errors as shown in Figure
and an exponential weighting factor dened as wmn p q
mpnq
where is a forgetting factor the weighted squared error is dened
as
mn
m
X
p
n
X
q
mpnq
ep q
Let us consider a QP ROS of order P Q for prediction as shown in
Figure and use the following notation for representing the prediction
FIGURE
ROSs in adaptive LP by the D RLS method While the image is scanned
from the position mn to mn a column of m pixels becomes available
as new information that may be used to update the forward and backward
predictors Observe that a part of the column of new information is hidden
by the ROS for both forward and backward prediction in the gure
coecients and the image data spanning the current ROS as vectors
amn
a
T
mn a
T
mn a
T
P
mn
T
where
a
p
mn amnp amnp amnpQ
T
with amn and
F
P
mn
f
T
m
n f
T
m
n f
T
mP
n
T
with
f
mp
n fm p n fm p n fm p nQ
T
Here the subscripts P and P represent the order size of the matrices
and vectors and the indices mn indicate that the values of the parameters
corresponding to the pixel location mn Observe that amn is a P
Q matrix with amnp q representing its element at p q With this
notation the prediction error may be written as
emn a
T
mn F
P
mn
The D RLS normal equations The coecients that minimize the
weighted sum of the squared prediction errors
mn given in Equation
are obtained as the solution to the D RLS normal equations which
are obtained as follows Let us perform partitioning of the matrices
amn and F
P
mn as
amn
"
amn
and
F
P
mn
fmn
"
F
P
mn
F
P
mn
fm P nQ
Observe that the coecient matrix
"
amn and the data matrix
"
F
P
mn
consist of all of the D RLS coecients amnp q and all of the image
pixels fm p n q such that p q QP ROS for a forward predictor
With partitioning as above the prediction error in Equation may be
written as
emn fmn
"
a
T
mn
"
F
P
mn
The sum of the squared prediction errors in Equation may now be
expressed as
mn
m
X
p
n
X
q
mpnq
ep q
m
X
p
n
X
q
mpnq
h
fp q
"
a
T
mn
"
F
P
p q
i
h
fp q
"
a
T
mn
"
F
P
p q
i
T
m
X
p
n
X
q
mpnq
h
f
p q
"
a
T
mn
"
F
P
p q fp q
"
a
T
mn
"
F
P
p q
"
F
T
P
p q
"
amn
i
In order to determine the coecients amnp q that minimize
mn
we could dierentiate the expression above for
mn with respect to the
coecient matrix
"
amn and equate the result to the null matrix of size
P Q which yields
mn
"
amn
m
X
p
n
X
q
mpnq
h
"
F
P
p q fp q
"
F
P
p q
"
F
T
P
p q
"
amn
i
Equation may be expressed in matrix notation as
m
X
p
n
X
q
mpnq
"
F
P
p q
h
fp q
"
F
T
P
p q
i
"
amn
In addition to the above using Equation in Equation we have
mn
m
X
p
n
X
q
mpnq
h
f
p q fp q
"
F
T
P
p q
"
amn
i
which may be written in matrix form as
mn
m
X
p
n
X
q
mpnq
fp q
h
fp q
"
F
T
P
p q
i
"
amn
Combining Equations and we get
m
X
p
n
X
q
mpnq
fp q
"
F
P
p q
h
fp q
"
F
T
P
p q
i
"
amn
mn
or
m
X
p
n
X
q
mpnq
F
P
p q F
T
P
p q amn
mn
which may be expressed as
P
mn amn mn
where
mn
mn
T
and
P
mn is the deterministic autocorrelation matrix of the weighted
image given by
P
mn
m
X
p
n
X
q
mpnq
F
P
p q F
T
P
p q
Equation represents the D RLS normal equations solving which we
can obtain the prediction coecients amnp q that adapt to the statistics
of the image at the location mn
Solving the D RLS normal equations Direct inversion of the auto
correlation matrix in Equation gives the desired matrix of prediction
coecients as
amn
P
mn mn
The matrix
P
mn is of size PQPQ the inversion
of such a matrix at every pixel mn of the image could be computationally
intensive Kuduvalli developed the following procedure to reduce the size
of the matrix to be inverted to Q Q The procedure starts with a
recursive relationship expressing the solution for the normal equations at the
pixel location mn in terms of that at m n Consider the expression
P
mn
m
X
p
n
X
q
mpnq
F
P
p q F
T
P
p q
mn
T
mn
T
P
mn
mn
mn
T
P
mn
P
mn
P
mn
PP
mn
where
rs
mn
m
X
p
n
X
q
mpnq
f
pr
q f
T
ps
q
Observe that
rs
mn
sr
m r s
which follows from the assumption that the image data have been windowed
such that fmn for m or n
The normal equations may now be expressed as
mn
T
mn
T
P
mn
mn
mn
T
P
mn
P
mn
P
mn
PP
mn
a
mn
a
mn
a
P
mn
mn
Q
Q
where
Q
is the null matrix of size Q
Equation may be solved in two steps First solve
mn
T
mn
T
P
mn
mn
mn
T
P
mn
P
mn
P
mn
PP
mn
I
Q
A
mn
A
P
mn
F
Q
mn
Q
Q
for the Q Q matrices A
p
mn p P and F
Q
mn
Here I
Q
is the identity matrix of size QQ and
Q
is the null
matrix of size Q Q Then obtain the solution to Equation
by solving
F
Q
mn a
mn mn
and using the relationship
A
p
mn a
mn a
p
mn p P
This approach is similar to the approach taken to solve the D YuleWalker
equations by the D Levinson method described in Section and leads
to a recursive algorithm that is computationally ecient the details of the
algorithm are given by Kuduvalli
Results of application to medical images Kuduvalli conducted
preliminary studies on the application of the D RLS algorithm to predictive
coding and compression of medical images In the application to coding the
value of the pixel at the current location mn is not available at the decoder
before the prediction coecient matrix amn is computed However the
prediction coecient matrix am n is available Thus for errorfree de
coding the a priori prediction error computed using the prediction coecient
matrix am n is encoded These error values which have a PDF that
is close to a Laplacian PDF may be eciently encoded using methods such
as the Human code Using a QP ROS of size and a forgetting factor
of Kuduvalli obtained an average bit rate of bpixel for
two of the images listed in Table this rate however is only marginally
lower than the bit rate of bpixel for the same two images obtained by
using the D Burg algorithm described in Section Although the D
RLS algorithm has the elegance of being a truly D algorithm that adapts
to the changing statistics of the image on a pixelbypixel basis the method
did not yield appreciable advantages in image data compression Regardless
the method has applications in other areas such as spectrum estimation and
ltering
Kuduvalli and Rangayyan performed a comparative analysis of several
image compression techniques including direct source coding transform cod
ing interpolative coding and predictive coding applied to the highresolution
digitized medical images listed in Table The average bit rates obtained
using several coding and compression techniques are listed in Table It
should be observed that decorrelation can provide signicant advantages over
direct source encoding of the original pixel data The adaptive predictive
coding techniques have performed better than the transform and interpola
tive coding techniques tested
In a study of the eect of sampling resolution on image data compression
Kuduvalli and Rangayyan prepared lowresolution versions of the images
listed in Table by smoothing and downsampling The results of the
application of the D Levinson predictive coding algorithm yielded average
bit rates of bpixel for images and bpixel for
images with the original images at bpixel This result
indicates that highresolution images possess more redundancy and hence
may be compressed by larger extents than their lowresolution counterparts
Therefore increasing the resolution of medical images does not increase the
amount of the related compressed data in direct proportion to the increase
in matrix size but by a lower factor This result could be a motivating
factor supporting the use of high resolution in medical imaging without undue
concerns related to signicant increases in datahandling requirements
See Aiazzi et al for a description of other methods for adaptive pre
diction and a comparative analysis of several methods for lossless image data
compression
Image Scanning Using the PeanoHilbert Curve
Peano scanning is a method of scanning an image by following the path de
scribed by a spacelling curve
Giuseppe Peano an Italian mathematician described the rst spacelling
curve in an attempt to map a line into a D space The term Peano
scanning is used to refer to such a scanning scheme irrespective of the space
lling curve used to dene the scan path Peano s curve was modied by
TABLE
Average Bit Rates Obtained in
the Lossless Compression of
the Medical Images Listed in
Table Using Several Image
Coding Techniques
Coding method Bits pixel
Original
Entropy H
Human
Arithmetic
LZW
DCT
Interpolative
D LP
D Levinson
D Burg
D RLS#
The Human code was used to encode the results of the transform interpola
tive and predictive coding methods #Only two images were compressed with
the D RLS method
Hilbert and the modied curve came to be known as the PeanoHilbert
curve
Moore studied the geometric and analytical interpretation of contin
uous spacelling curves Spacelling curves have aided in the development of
fractals see Section for a discussion on fractals The PeanoHilbert
curve has been applied to display continuoustone images in order to
eliminate deciencies of the ordered dithered technique such as Moir$e fringes
Lempel and Ziv used the PeanoHilbert curve to scan images and dene
the lowest bound of compressibility Zhang et al explored the
statistical characteristics of medical images using Peano scanning
Provine and Rangayyan studied the application of Peano scan
ning for image data compression with an additional step of decorrelation
using dierentiation orthogonal transforms or LP the following paragraphs
describe the basics of the methods involved and the results obtained in their
work
Denition of the Peanoscan path
If a physical scanner that can scan an image by following the Peano curve
is not available Peano scanning may be simulated by selecting pixels from
a rasterscanned image by traversing the D data along the path described
by the PeanoHilbert curve The reordered pixel data so obtained in a D
stream may be subjected to decorrelation and encoding operations as desired
An inverse Peanoscanning operation would be required at the receiving end
to reconstruct the original image A general image compression scheme as
above is summarized in Figure
FIGURE
Image compression using Peano scanning
The Peanoscan operation is recursive in nature and spans a D space
encountering a total of
i
i
points where i is an integer and i From the
perspective of processing a D array containing the pixel values of an image
this would require that the dimensions of the array be an integral power of
In the scanning procedure the given image of size
n
n
is divided into four
quadrants each of them forming a subimage see Figure Each of the
subimages is further divided into four quadrants and the procedure continues
The original image is divided into a total of T
i
ni
subimages each
of size
i
i
where i n In the following discussion the
subimages formed by the recursive subdivision procedure as above will be
referred to as s
i
k k T
i
where k increases along the direction of
the scan path the entire image will be referred to as s
n
Thus each of the four
quadrants formed by partitioning a subimage s
i
for any i is of size
i
i
The division of a given image into subimages is shown in Figure the
recursive division of subimages is performed until the s
subimages are formed
The four pixels within the smallest subimage are denoted as p p p
and p respectively in the order of being scanned As the scan path builds
recursively the path denition is based on the basic denitions for a
subimage as well as the recursive denitions for subimages of larger size until
the entire image is scanned The four basic denitions of the Peanoscanning
operation are given in Figure
The recursive denitions of the Peanoscanning operation are given in Fig
ure which inherently use the basic denitions shown in Figure
to obtain further pixels from the subimages The denitions go down recur
sively from i n to i that is from the image s
n
down to the subimages
s
the basic denitions are used to obtain the pixels from the s
subimages
The scanpath denition for an image or a subimage depends on i At a
higher level that is for an s
i
i the recursive denition is as follows
If i is odd the recursive denition is R see Figure
If i is even the recursive denition is D
From the recursive denitions shown in Figure the denitions of the
scan pattern in each of the subimages is obtained as follows
If s
i
k has the recursive denition R s
i
k will follow the path
given by D s
i
k and s
i
k the path R and s
i
k
the path U
If s
i
k has the recursive denition D s
i
k will follow the path
given by R s
i
k and s
i
k the path D and s
i
k
the path L
If s
i
k has the recursive denition L s
i
k will follow the path
given by U s
i
k and s
i
k the path L and s
i
k
the path D
If s
i
k has the recursive denition U s
i
k will follow the path
given by L s
i
k and s
i
k the path U and s
i
k
the path R
FIGURE
Division of a image into subimages during Peano scanning Repro
duced with permission from JA Provine and RM Rangayyan Lossless
compression of Peanoscanned images Journal of Electronic Imaging
c
SPIE and IS%T
FIGURE
Basic denitions of the Peanoscanning operation The points marked pp
represent the four pixels in a subimage Each scan pattern shown visits
four pixels in the order indicated by the arrows R right L left D down
U up Reproduced with permission from JA Provine and RM Rangayyan
Lossless compression of Peanoscanned images Journal of Electronic Imag
ing
c
SPIE and IS%T
FIGURE
Recursive denitions of the Peanoscanning operation Reproduced with per
mission from JA Provine and RM Rangayyan Lossless compression of
Peanoscanned images Journal of Electronic Imaging
c
SPIE and IS%T
The index k can take any value in the range T
i
for i n
From the recursive denitions it is evident that k cannot continuously increase
horizontally or vertically due to the nature of the scan Furthermore because
the scan path takes its course recursively except for i n the recursive
denitions L and U see Figure are also possible for lower values of
i All the subimages are divided and dened recursively until all the subimages
s
are dened The basic denitions of the Peano scan are then followed for
each of the s
subimages The Peanoscan pattern for a image is
illustrated in Figure where the heavy dots indicate the positions of
the rst pixels scanned Understanding the scan pattern is facilitated by
viewing Figure along with Figure
FIGURE
Peanoscan pattern for a image The positions of the pixels on the scan
pattern are shown by heavy dots in the rst subimage Reproduced with
permission from JA Provine and RM Rangayyan Lossless compression of
Peanoscanned images Journal of Electronic Imaging
c
SPIE and IS%T
Implementation of Peano scanning in software could use the recursive nature
of the scan path eciently to obtain the pixel stream Recursive functions
may call themselves within their body as the image is divided progressively
into subimages until the s
subimages are formed and the pixels are obtained
recursively as the function builds back from the s
subimages to the full image
s
n
Thus the D image is unwrapped into a D data stream by following a
continuous scan path
The inverse Peanoscanning operation accomplishes the task of lling up the
D array with the D data stream This operation corresponds to the original
works of Peano and Hilbert where the continuous mapping of a
straight line into a D plane was described Because the Peanoscanning
operation is reversible no loss of information is incurred
Properties of the PeanoHilbert curve
The PeanoHilbert curve has several interesting and useful properties The
curve is continuous but not dierentiable it does not have a tangent at any
point Moore gave an explanation of the PeanoHilbert curve adhering
to this property This property motivated the development of several other
curves with the same property which are used in the domain of fractals
The PeanoHilbert curve lls the D space continuously without passing
through any point more than once This feature enables the mapping of a D
array into a D data stream The recursive nature of the curve is useful in
ecient implementation of the path of the curve These two properties aid in
scanning an image recursively quadrant by quadrant leaving each quadrant
only after having obtained every pixel within that quadrant with each pixel
visited only once in the process see Figures and Preservation of
the local D context in the scanning path could be expected to increase the
correlation between successive elements in the D data stream This aspect
could facilitate improved image data compression
Two other aspects of the PeanoHilbert curve have proven to be useful
in the bilevel display of continuoustone images Linearizing a
D array along the path described by the PeanoHilbert curve reduces the
error between the sum of the bilevel values and the sum of the continuous
tone values of the original image because D locality is maintained by the scan
path unlike the D vector formed by concatenating the horizontal rasterscan
lines of the image The problem of long sections of scan lines running adjacent
to one another is eliminated by following the Peanoscan path instead of the
raster scan Thus Moir$e patterns can be eliminated in regions of uniform
intensity when presenting graylevel images on a bilevel display
Implementation of Peano scanning
A practical problem that could arise in implementing the Peanoscanning op
eration on large images is the diculty in allocating memory for the long linear
array used within the body of the recursive function for storing the scanned
data Provine suggested the following approach to address this problem
by using a symmetrical pattern exhibited by the PeanoHilbert curve
The PeanoHilbert curve exhibits a symmetrical pattern which may be
described as follows For any subimage s
i
k the scan paths for the subimages
s
i
k and s
i
k are the mirror re!ections of the paths for the
subimages s
i
k and s
i
k respectively Two types of symmetry
exist in the Peanoscan paths for a
i
i
subimage depending on whether i
is odd or even If i is odd the pattern for the upper half of the D space is
re!ected in the lower half if i is even the pattern in the lefthand half of the
D space is re!ected in the righthand half see Figure
A symmetrical scan pattern exists for any subimage formed by the recursive
division process Hence for the smallest subimage the symmetrical pattern
suggests that the basic denitions eectively obtain only two pixels one after
the other either horizontally or vertically In other words the basic scan
pattern eectively obtains only two pixels p and p out of the four pixels
in a subimage see Figure The manner in which the remaining
two pixels p and p are obtained follows the symmetry property stated above
substituting i and k The sequence in which the two pixels p and
p are obtained is shown in Figure in dashed lines
In scanning large images for any subimage s
i
k the pattern of the Peano
scan path from the rst pixel of s
i
k to the last pixel of s
i
k is
the same as that from the last pixel of s
i
k to the rst pixel of s
i
k
Because the Peanoscan path does not leave any quadrant without visiting all
the pixels within the quadrant two equal sections of the image can be un
wrapped independently without aecting each other into individual linear
arrays by following the same scan path but in opposite directions The re
sulting linear arrays when concatenated appropriately give the required D
data stream For the case illustrated in Figure a the pixels can be
obtained as a linear sequence by tracing the Peanoscan path on pixels
in the upper half followed by the pixels in the lower half Thus several
D arrays of a reasonable size may be used to hold the pixels obtained from
dierent sections of the image After all the subimages have been scanned in
parallel if desired the arrays may be concatenated appropriately to form
the long sequence containing the entire image data
Decorrelation of Peanoscanned data
The Peanoscanning operation scans the given picture recursively quadrant by
quadrant Therefore we could expect the D local statistics to be preserved
in the resulting D data stream Furthermore we could also expect a higher
correlation between pixels for larger lags in the Peanoscanned data than
between the pixels obtained by concatenating the rasterscanned lines into a
D array of pixels
a
b
FIGURE
Symmetrical patterns exhibited by the PeanoHilbert curve for a an
subimage and b a subimage The line AB indicates the axis of
symmetry in each case The pixels in the image in a are labeled in the
order of being scanned the pixels in the upper half followed by the
pixels in the lower half of the subimage in a Figure courtesy of JA
Provine
FIGURE
Symmetrical patterns in the basic denitions of the Peano scan Figure cour
tesy of JA Provine
Figure shows the ACFs obtained for the Lenna test image and a mam
mogram for rasterscanned and Peanoscanned pixel streams As expected
Peano scanning has maintained higher interpixel correlation than raster scan
ning Similar observations have been made by Zhang et al in their
study on the stochastic properties of medical images and data compression
with Peano scanning
The simplest method for decorrelating pixel data is to produce a data se
quence containing the dierences between successive pixels As shown earlier
in Section and Figure the dierentiated data may be expected to
have a Laplacian PDF which is useful in compressing the data Provine and
Rangayyan applied a simple rstdierence operation to decorre
late Peanoscanned image data and encoded the resulting values using the
Human arithmetic and LZW coding schemes In addition they applied the
D DCT and LP modeling procedures to rasterscanned and Peanoscanned
data streams as well as the D DCT and LP modeling procedures to the
original image data Some of the results obtained by Provine and Rangayyan
are summarized in Tables and The application of either
the Human or the arithmetic code to the dierentiated Peanoscanned data
stream resulted in the lowest average bit rate in the study
Image Coding and Compression Standards
Two highly recognized international standards for the compression of still
images are the Joint Bilevel Image experts Group JBIG standard
and the Joint Photographic Experts Group JPEG stan
dard JBIG and JPEG are sanctioned by the International Orga
nization for Standardization ISO and the Comit$e Consultatif International
T$el$ephonique et T$el$egraphique CCITT Although JBIG was initially pro
posed for bilevel image compression it may also be applied to continuous
a
b
FIGURE
ACF of rasterscanned and Peanoscanned pixels plotted as a function of the
distance lag between the scanned pixels a for the Lenna image and b for
a mammogram Figure courtesy of JA Provine
TABLE
Average Bit Rate with the Application of the Human Arithmetic and
LZW Coding Schemes to Rasterscanned and Peanoscanned Data
Obtained from Eight Test Images
Entropy Average number of bitspixel
Image Size H
LZW
bpixel pixels bits Human Arith Raster Peano
Airplane
Baboon
Cameraman
Lenna
Lenna
Peppers
Sailboat
Tiany
Mean
SD
See also Tables and Note Arith arithmetic coding
tone images by treating bit planes as independent bilevel images Note
The term continuoustone images is used to represent graylevel images
color images and multicomponent images whereas some authors use the
term mary for the same purpose the former is preferred as it is used
by JPEG The eciency of such an application depends upon preprocess
ing for bitplane decorrelation The Moving Picture Experts Group MPEG
standard applies to the compression of video images In the
context of medical image data handling and PACS the ACR and the US
National Electrical Manufacturers Association NEMA proposed standards
known as the ACR NEMA and DICOM Digital Imaging and Communica
tions in Medicine standards The following sections
provide brief reviews of the standards mentioned above
TABLE
Average Bit Rate with Dierentiated Peanoscanned Data PD
Compared with the Results of D and D DPCM Encoding of
Rasterscanned Data from Eight Test Images
Average number of bitspixel
Human Arithmetic LZW
DPCM DPCM DPCM
Image PD D D PD D D PD D D
Airplane
Baboon
Cameraman
Lenna
Lenna
Peppers
Sailboat
Tiany
Mean
SD
See also Tables and
The JBIG Standard
JBIG is a standard for progressive coding of bilevel images that supports three
coding modes progressive coding compatible progressive sequential coding
and singlelayer coding A review of the singlelayer coding mode is presented
in the following paragraphs
For a bilevel image bmn m M n N a typical JBIG coding
scheme includes four main functional blocks as shown in Figure The
following items form important steps in the JBIG coding procedure
The typical prediction step is a lineskipping algorithm A given line is
marked as typical if it is identical to the preceding line The encoder
adds a special label to the encoded data stream for each typical line
instead of encoding the line The decoder generates the pixels of typical
lines by line duplication
TABLE
Average bit rates for eight test images with several decorrelation and
encoding methods
Average bit rate bitspixel
Direct D linear D
arith PD PD predictive DCT
Peano Human arith coding coding
Image or raster encoding encoding raster raster
Airplane
Baboon
Cameraman
Lenna
Lenna
Peppers
Sailboat
Tiany
Mean
SD
See also Tables and Note arith arithmetic coding PD
Dierentiated Peanoscanned data
The adaptive templates block provides substantial coding gain by look
ing for horizontal periodicity in the bilevel image When a periodic
template is changed the encoder multiplexes a control sequence into
the output data stream
The model templates block is a context arithmetic coder The context is
determined by ten particular neighboring pixels that are dened by two
model templates the threeline template and the twoline template as
shown in Figure labeled as &X and &A where &A is the adaptive
pixel whose position could be varied during the coding process
The adaptive arithmetic encoder is an entropy coder that determines
the necessity of coding a given pixel based upon the outputs of the
typical prediction block and the model templates block If necessary
the encoder notes the context and uses its internal probability estimator
to estimate the conditional probability that the current pixel will be of
a given value
It should be noted that the JBIG coding algorithm includes at least three
decorrelation steps for a given bilevel image or bit plane
FIGURE
The singlelayer JBIG encoder Reproduced with permission from L Shen
and RM Rangayyan Lossless compression of continuustone images by com
bined interbitplane decorrelation and JBIG coding Journal of Electronic
Imaging
c
SPIE and IS%T
FIGURE
Two contextmodel templates used in the singlelayer JBIG encoder ' Pixel
being encoded X Pixels in the context model A Adaptive pixel Repro
duced with permission from L Shen and RM Rangayyan Lossless compres
sion of continuustone images by combined interbitplane decorrelation and
JBIG coding Journal of Electronic Imaging
c
SPIE
and IS%T
In a D form of JBIG coding runlength coding is used to encode each
line in the image by using a variablelength coding scheme see Gonzalez and
Woods for further details of this and other JBIG procedures See Sections
and for discussions on the results of application of JBIG
The JPEG Standard
JPEG is a continuoustone image compression standard that supports both
lossless and lossy coding The standard includes three sys
tems
A lossy baseline coding system based upon blockwise application of the
DCT
An extended coding system for greater compression higer precision and
progressive transmission and recovery
A lossless independent coding system for reversible compression
In the baseline coding system the pixel data are limited to a precision
of b The image is broken into blocks shifted in gray level and
transformed using the DCT The transform coecients are quantized with
variablelength code assignment to a maximum of b Due to the presence
of block artifacts and other errors in lossy compression this mode of JPEG
would not be suitable for the compression of medical images
The JPEG standard is based upon the wavelet transform
JPEG oers the option of progressive coding from lossy toward lossless
as well as the option of coding ROIs with higher quality than that for the
other regions in the given image which may be of interest in some
applications
The lossless mode of JPEG uses a form of predictive DPCM coding A
linear combination of each pixel s neighbors at the left upper and upper
left positions is employed to predict the pixel s value and then the dierence
between the true value of the pixel and its predicted value is coded through
an entropy coder such as the Human or arithmetic coder Lossless JPEG
denes seven linear combinations known as prediction selection values PSV
For a continuoustone image fmn m M n N the predic
tors used for an interior pixel fmn m M n N in lossless
JPEG coding are as follows
PSV no prediction or
"
fmn which indicates entropy encod
ing of the original image directly
PSV
"
fmn fm n
PSV
"
fmn fmn
PSV
"
fmn fm n
PSV
"
fmn fm n fmn fm n
PSV
"
fmn fm n fmn fm n
PSV
"
fmn fmn fm n fm n
PSV
"
fmn fm n fmn
where
"
fmn is the predicted value for the pixel fmn The boundary
pixels could be treated through various possible ways which will not make a
signicant dierence in the nal compression ratio due to their small popula
tion A typical method is to use PSV for the rst row and PSV for the
rst column with special treatment of the rst pixel at position using
the value of
K
where K is the number of precision bits of the pixels
Sung et al evaluated the application of JPEG for the compres
sion of mammograms and suggested that compression ratios of up to
were possible without visual loss preserving signicant medical information
at a condence level of It was also suggested that compression of up to
could be achieved without aecting clinical diagnostic performance
See Sections and for discussions on the results of applica
tion of lossless JPEG to several test images
The MPEG Standard
The MPEG standard includes several schemes for the compression of video
images for various applications based upon combinations of the DCT and
DPCM including motion compensation The techniques ex
ploit data redundancy and correlation within each frame as well as between
frames furthermore they take advantage of certain psychovisual properties of
the human visual system Recent versions of MPEG include special features
suitable for videoconferece multimedia streaming media and videogame
systems Most of such special features are not of relevance in lossless
compression of biomedical image data
The ACR NEMA and DICOM Standards
The proliferation of medical imaging technology and devices of several types
in the s and s led to a situation where due to the lack of standards
interconnection and communication of data between imaging and comput
ing devices was not possible In order to rectify this situation the ACR
and NEMA established the ACR NEMA standard on dig
ital imaging and communication specifying the desired hardware interface
a minimum set of software commands and a consistent set of data formats
to facilitate communication between imaging devices and computers across
networks This was followed by another standard on data compression the
ACR NEMA PS specifying the manner in which header data
were to be provided such that a recipient of the compressed data could iden
tify the data compression method and parameters used and reconstruct the
image data The standard permits the use of several image decorrelation and
data compression techniques including transform DCT predictive DPCM
Human and LempelZiv coding techniques The DICOM standard
includes a number of enhancements to the ACR NEMA standard including
conformance levels and applicability to a networked environment
Segmentationbased Adaptive Scanning
Shen and Rangayyan proposed a segmentationbased lossless image
coding SLIC method based on a simple but ecient regiongrowing proce
dure An embedded regiongrowing procedure was used to produce an adap
tive scanning pattern for the given image with the help of a discontinuityindex
map that required a small number of bits for encoding The JBIG method was
used for encoding both the errorimage data and the discontinuityindex map
data The details of the SLIC method and the results obtained are described
in the following paragraphs
Segmentationbased coding
Kunt et al proposed a contourtexture approach to picture coding they
called such approaches secondgeneration image coding techniques The
main idea behind such techniques is to rst segment the image into nearly
homogeneous regions surrounded by contours such that the contours corre
spond as much as possible to those of the objects in the image and then
to encode the contour and texture information separately Because contours
can be represented as D signals and the pixels within a region are highly
correlated such methods are expected to lead to high compression ratios
Although the idea appears to be promising its implementation meets with a
series of diculties A major problem exists at its very important rst step
segmentation which determines the nal performance of the segmentation
based coding method It is well recognized that there are no satisfactory
segmentation algorithms for application to a wide variety of general images
Most of the available segmentation algorithms are sophisticated and give good
performance only for specic types of images
In order to overcome the problem mentioned above in relation to segmen
tation Shen and Rangayyan proposed a simple regiongrowing
method In this procedure instead of generating a contour set a discontinu
ity map is obtained during the regiongrowing procedure Concurrently the
method also produces a corresponding error image based upon the dierence
between each pixel and its corresponding center pixel The discontinuity
map and the error image are then encoded separately
Regiongrowing criteria
The aim of segmentation in image compression is not the identication of
objects or the analysis of features instead the aim is to group spatially
connected pixels lying within a small graylevel dynamic range The region
growing procedure in SLIC starts with a single pixel called the seed pixel (
in Figure Each of the seed s connected neighboring pixels from (
to ( in the order as shown in Figure is checked with a regiongrowing
or inclusion condition If the condition is satised the neighboring pixel is
included in the region The four neighbors of the newly added neighboring
pixel are then checked for inclusion in the region This recursive procedure is
continued until no spatially connected pixel meets the growing condition A
new regiongrowing procedure is then started with the next pixel in the image
that is not already a member of a region the procedure ends when every pixel
in the image has been included in one of the regions grown
FIGURE
Demonstration of a seed pixel ( and its connected neighbors ( to
( in region growing Reproduced with permission from L Shen and RM
Rangayyan A segmentationbased lossless image coding method for high
resolution medical image compression IEEE Transactions on Medical Imag
ing
c
IEEE
The regiongrowing conditions used in SLIC are the following
The neighboring pixel is not a member of any of the regions already
grown
The absolute dierence between the neighboring pixel and the corre
sponding center pixel is less than the limit errorlevel to be dened
later
Figure demonstrates the relationship between a neighboring pixel
and its center pixel at the specic stage illustrated in the gure pixel A
has become a member of the region being grown and its four neighbors
namely B C D and E are being checked for inclusion E is already a
member of the region Under this circumstance pixel A is the center pixel
of the neighboring pixels B C D and E When a new neighboring pixel
is included in the region being grown its errorlevelshiftup dierence with
respect to its center pixel is stored as the pixel s error value If only the
rst of the two regiongrowing conditions is met the discontinuity index of the
pixel is incremented By this process after region growing a discontinuity
index image data part and an errorimage data part will be obtained The
maximum value of the discontinuity index is Most of the segmentation
based coding algorithms reported in the literature include contour coding and
region coding instead of these steps the SLIC method uses a discontinuity
index map and an errorimage data part
The errorlevel in SLIC is determined by a preselected parameter error
bits as
errorlevel
errorbits
For instance if the error image is allowed to take up to bpixel errorbits
the corresponding errorlevel is the allowed dierence range is then
The error value of the seed pixel of each region is dened as the
value of its lower errorbits bits the value of the higher N errorbits bits
of the pixel is stored in a highbits seeddata part where N is the number
of bits per pixel in the original image data
The three data parts described above are used to fully recover the original
image during the decoding process The regiongrowing conditions during
decoding are that the neighboring pixel under consideration for inclusion be
not in any of the previously grown regions and that its discontinuity index
equal When the conditions are met for a pixel its value is restored as the
sum of its errorlevelshiftdown error value and its center pixel value except
for the seed pixels of every region If only the rst of the two conditions is
satised the discontinuity index of that pixel is decremented Thus the
discontinuity index generated during segmentation is used to guide region
growing during decoding The highbits seeddata part is combined with
the errorimage data part to recover the seed pixel value of each region
Figure provides a simple example for illustration of the regiongrowing
procedure and its result The image in the example shown in Figure
FIGURE
Demonstration of a neighboring pixel B C D or E being checked for inclu
sion against the current center pixel A during region growing Reproduced
with permission from L Shen and RM Rangayyan A segmentationbased
lossless image coding method for highresolution medical image compression
IEEE Transactions on Medical Imaging
c
IEEE
a is a section of an eye of the Lenna image shown in Figure
a The value of errorbits was set to for this example Figure
b shows the result of region growing The corresponding three data
parts namely the discontinuityindex image data errorimage data and high
bits seed data are shown in Figure c d and e respectively The
full Lenna image and the corresponding discontinuityindex image
data scaled and errorimage data scaled are shown in Figure
The SLIC procedure
The complete SLIC procedure is illustrated in Figure At the encoding
end the original image is transformed into three parts discontinuityindex
image data errorimage data and highbits seed data by the regiongrowing
procedure The rst two data parts are encoded using the Gray code see
Table broken down into bit planes and nally encoded using JBIG
The last data part is stored or transmitted as is it needs only N errorbits
bits per region
At the decoding end the JBIGcoded data les are JBIGdecoded rst and
then the Graycoded bit planes are composed back to binary code Finally
the three parts are combined together by the same regiongrowing procedures
as before to recover the original image
A lossless compression method basically includes two major stages one is
image transformation with the purpose of data decorrelation the other is
encoding of the transformed data However in the SLIC procedure image
transformation is achieved in both the regiongrowing procedure and later
in the JBIG procedure whereas encoding is accomplished within the JBIG
procedure
Results of image data compression with SLIC
Five mammograms and ve chest radiographs were used to evaluate the per
formance of SLIC The procedure was tested using bpixel and
bpixel versions of the images obtained by direct mapping and by discard
ing the two leastsignicant bits respectively from the original b images
The method was also tested with commonly used b nonmedical test im
ages The performance of SLIC was compared with that of JBIG
JPEG adaptive LempelZiv ALZ coding HINT and D
LP coding In using JBIG for direct encoding of the image or for
encoding the errorimage data part and the discontinuityindex data part pa
rameters were selected so as to use three lines of the image NLPS in
the underlying model D and no progressive spatial resolution buildup
The lossless JPEG package used in the study includes features of automatic
determination of the best prediction pattern and optimal Human table gen
eration The UNIX utility compress was used for ALZ compression For the
b test images the D Burg LP algorithm see Section followed by
FIGURE
A simple example of the regiongrowing procedure and its result with error
bits set to be a Original image with the size of rows by columns a
section of the Lenna image shown in Figure b The result
of region growing The seed pixel of every region grown is identied a few
regions include only the corresponding seed pixels c The discontinuityindex
image data part d The errorimage data part e The highbits seeddata
part from left to right Reproduced with permission from L Shen and
RM Rangayyan A segmentationbased lossless image coding method for
highresolution medical image compression IEEE Transactions on Medical
Imaging
c
IEEE
a
b c
FIGURE
The Lenna image and the results with SLIC with errorbits
a Original image b The discontinuityindex image data part scaled
c The errorimage data part scaled See also Figure Reproduced
with permission from L Shen and RM Rangayyan A segmentationbased
lossless image coding method for highresolution medical image compression
IEEE Transactions on Medical Imaging
c
IEEE
FIGURE
Illustration of the segmentationbased lossless image coding SLIC proce
dure Reproduced with permission from L Shen and RM Rangayyan A
segmentationbased lossless image coding method for highresolution medical
image compression IEEE Transactions on Medical Imaging
c
IEEE
Human error coding DBurgHuman was utilized instead of the JPEG
package The DBurgHuman program was designed specically for b
images the JPEG program permits only b images
The SLIC method has a tunable parameter which is errorbits The data
compression achieved was found to be almost the same for most of the b
medical images by using the value or for errorbits Subsequently error
bits was used in the remaining experiments with the b versions of the
medical images The performance of the SLIC method is summarized in Ta
ble along with the results of ALZ JBIG JPEG and HINT compression
It is seen that SLIC has outperformed all of the other methods studied with
the testimage set used except in the case of one mammogram and one chest
radiograph for which JBIG gave negligibly better results On the average
SLIC improved the bit rate by about and as compared with
JBIG JPEG and HINT respectively
TABLE
Average Bits Per Pixel with SLIC errorbits ALZ JBIG JPEG
Best Mode and HINT Using Five b Mammograms m to m and
m and Five b Chest Radiographs c to c and c by BitsPixel
Image Entropy
bpixel H
ALZ JBIG JPEG HINT SLIC
m
m
m
m
c
c
c
c
m
c
Average
The lowest bit rate in each case is highlighted Reproduced with permission
from L Shen and RM Rangayyan A segmentationbased lossless image
coding method for highresolution medical image compression IEEE Trans
actions on Medical Imaging
c
IEEE
Whereas the SLIC procedure performed well with highresolution medical
images its performance with lowresolution general images was comparable
to the performance of JBIG and JPEG as shown in Table When
SLIC used an optimized JBIG algorithm with the maximum number of lines
that is NLPS row the number of rows of the image the last column of
Table in the inherent model instead of only three lines NLPS in
Table its performance was better than that of JBIG and comparable
to that of JPEG
TABLE
Average Bits Per Pixel with SLIC errorbits ALZ JBIG and
JPEG Best Mode Using Eight Commonly Used b Images
SLIC
Image Entropy NLPS
bpixel order ALZ JBIG JPEG row
Airplane
Baboon
Cameraman
Lenna
Lenna
Peppers
Sailboat
Tiany
Average
Reproduced with permission from L Shen and RM Rangayyan A
segmentationbased lossless image coding method for highresolution medi
cal image compression IEEE Transactions on Medical Imaging
c
IEEE
In the studies of Shen and Rangayyan setting errorbits was
observed to be a good choice for the compression of the b versions of the
medical images Table lists the results of compression with the ALZ
JBIG DBurgHuman HINT and SLIC methods The SLIC technique
has provided lower bit rates than the other methods The average bit rate
is bpixel with SLIC whereas HINT JBIG and DBurgHuman have
average bit rates of and bpixel respectively
TABLE
Average Bits Per Pixel with SLIC errorbits ALZ JBIG
DBurgHuman DBH and HINT Using Five b Mammograms
m to m and m and Five b Chest Radiographs c to c and c
Image Entropy
bpixel order ALZ JBIG DBH HINT SLIC
m
m
m
m
c
c
c
c
m
c
Average
The lowest bit rate in each case is highlighted Reproduced with permission
from L Shen and RM Rangayyan A segmentationbased lossless image
coding method for highresolution medical image compression IEEE Trans
actions on Medical Imaging
c
IEEE
Most of the previously reported segmentationbased coding techniques
involve three procedures segmentation contour coding and tex
ture coding For segmentation a sophisticated procedure is generally em
ployed for the extraction of closed contours In the SLIC method a simple
singlescan neighborpixel checking algorithm is used The major problem af
ter image segmentation with the other methods lies in the variance of pixel
intensity within the regions which has resulted in the application of such
methods to lossy coding instead of lossless coding In order to overcome
this problem the SLIC method uses the mostcorrelated neighboring pixels
to generate a lowdynamicrange error image during segmentation Contour
coding is replaced by a coding step applied to a discontinuityindex map with
a maximum of levels and texture coding is turned into the encoding of a
lowdynamicrange error image Further improvement of the performance of
SLIC may be possible by the application of more ecient coding methods
to the error image and discontinuityindex data parts and by modifying the
regiongrowing procedure
The SLIC method was extended to the compression of D biomedical im
ages by Lopes and Rangayyan The performance of the method varied
depending upon the nature of the image However the advantage of im
plementation of SLIC in D versus D on a slicebyslice basis was not sig
nicant For example both D and D SLIC with the compression utility
bzip as the encoding scheme after decomposition of the image by the
segmentationbased procedure reduced the data in a D CT head examina
tion from the original bvoxel size to bvoxel the zerothorder entropy
of the image was bvoxel The inclusion of the SLIC procedure improved
the performance of the compression utilities gzip and compress by
Acha et al extended the SLIC method to the compression of color
images in the application of diagnosis of burn injuries Color images of size
pixels at bpixel in the RGB format were compressed to the
rate of bpixel application of the JPEG lossless method resulted in a
rate of bpixel In a further study Serrano et al converted the
same images as in the study of Acha et al from the RGB format to
the YIQ luminance inphase and quadrature system in a lossless manner
and showed that the bit rate could be further reduced by the use of SLIC to
bpixel
Enhanced JBIG Coding
The SLIC procedure demonstrates a successful example of the application of
the JBIG method for coding the discontinuityindex map and errordata parts
We may therefore expect the incorporation of decorrelation procedures such
as predictive algorithms into JBIG coding to provide better performance
Shen and Rangayyan proposed a combination of multiple decor
relation procedures including a lossless JPEGbased predictor a transform
based interbitplane decorrelator and a JBIGbased intrabitplane decorre
lator the details of their procedure are described in the following paragraphs
Although JBIG includes an ecient intrabitplane decorrelation procedure
it needs an ecient preprocessing algorithm for interbitplane decorrelation
in order to achieve good compression eciency with continuoustone images
There are several ways in which bit planes may be created from a continuous
tone image One common choice is to use the bits of a foldedbinary or Gray
representation of intensity The Gray code is the most common alternative
to the binary code for representing digital numbers see Table The
major advantage of the Gray code is that only one bit changes between each
pair of successive code words which is useful to provide a good bitplane
representation for original image pixels most neighboring pixels have highly
correlated and close values and thus most neighboring bits within each Gray
coded bit plane may be expected to have the same value It has been shown
that by using the Gray representation JBIG can obtain compression ratios
at least comparable to those of lossless JPEG
Coding schemes other than the Gray code may be derived based on specic
requirements For instance for a Kbit image fmn the prediction error
emn could be represented as
emn fmn
"
fmn
where
"
fmn is the predicted value of the original pixel fmn In general
up to K bits could be required to represent the dierence between two
Kbit numbers However because the major concern in compression is to
retrieve fmn from emn and
"
fmn we could make use of the following
binary arithmetic operation
"emn emn % f
z
K bits
g
b
where % is the bitwise AND operation and the subscript b indicates binary
representation Then only K bits are necessary to represent the prediction
error "emn The original pixel value may be retrieved as
fmn "emn
"
fmn % f
z
K bits
g
b
With the transformation as above the value of
K
v for "emn denotes
a prediction error emn of either
K
v or v In general the lower and
the higher ends of the value of the error "emn appear more frequently than
midrange values
The F
transformation could make bitplane coding more ecient by
increasing the run length in the mostsignicant bits
v
F
v
v
v v
K
K
v v
K
with the inverse transform given by
v F
v
v
K
p v
p even
q v
q odd
The F
transformation has a similar function but through the
reversal of the higherhalf of the value range
v
F
v
v v
K
v f
z
K bits
g
b
v
K
where is the bitwise exclusive OR operation its inverse transform is
v F
v
F
v
Shen and Rangayyan investigated the combined use of the PSV sys
tem in JPEG see Section and bitplane coding using JBIG along
with one of the Gray F
and F
transforms In using JBIG for direct coding
of the original image or for coding the prediction error image the method
was parameterized to use the threeline model template and stripe size equal
to the number of rows of the image with no progressive spatial resolution
buildup see Section The lossless JPEG scheme was set to generate
the optimal Human table for each PSV value
The methods were tested with a commonly used set of eight images see
Table A comparison of the compression eciencies of lossless JPEG
JBIG and PSVincorporated JBIG with one of the Gray F
or F
trans
forms is shown in Figure for one of the test images Each group of
vertical bars in each case consists of ve bars corresponding to the entropy
of the prediction error image followed by the actual bit rates with lossless
JPEG coding and PSVincorporated JBIG bitplane coding with the Gray
F
and F
transformations respectively In the gure there are three hor
izontal lines representing the zerothorder entropy of the original image the
bit rate by direct JBIG coding with the Gray transformation and the best
bit rate among all of the methods tested The best performance for each test
image was achieved with PSVincorporated JBIG bitplane coding with the
F
transformation and PSV except for the Lenna image for
which the best rate was given with PSV and JBIG bitplane coding of the
prediction error after F
transformation The F
and F
transforms provided
similar performance and performed better than the Gray transform with the
F
transform giving slightly lower bit rates in most of the cases
In the results obtained by Shen and Rangayyan it was observed that
the zerothorder entropies of the prediction error images were much lower than
those of the original images that the bit rates by lossless JPEG compression
were always higher than the zerothorder entropies of the prediction error
images and that PSVincorporated JBIG bitplane coding provided bit rates
lower than the zerothorder entropies of the prediction error images with
the exceptions being the Baboon and Sailboat images These observations
show that a simple prediction procedure such as the PSV scheme employed
in lossless JPEG is useful for decorrelation In particular prediction with
PSV followed by the F
transform and JBIG bitplane coding achieves an
average bit rate that is about bpixel lower than those achieved with direct
Graycoded JBIG compression and the optimal mode of lossless JPEG This
also indicates that the F
transform is a better interbitplane decorrelator
than the Gray code for prediction error images and that the intrabitplane
decorrelation steps within the JBIG algorithm are not redundant with prior
decorrelation by the PSV system in JPEG
FIGURE
Comparison of image compression eciency with the enhanced JBIG scheme
The ve bars in each case represent the zerothorder entropy of the prediction
error image and actual bit rates with lossless JPEG and PSVincorporated
JBIG with the Gray F
and F
transformation respectively from left to
right Reproduced with permission from L Shen and RM Rangayyan
Lossless compression of continuustone images by combined interbitplane
decorrelation and JBIG coding Journal of Electronic Imaging
c
SPIE and IS%T
Shen and Rangayyan applied the best mode of the PSVincorporated
JBIG bitplane coding scheme the PSV predictor followed by JBIG bit
plane coding of F
transformed prediction error denoted as PSVFJBIG
to the JPEG standard set of continuoustone test images the results are
shown in Table The table also shows the zerothorder entropies of the
prediction error images the results of the best mode of lossless JPEG coding
for b component images only the results of direct Graytransformed JBIG
coding and the results of the best mode of the CREW technique Compression
with Reversible Embedded Wavelets one of the methods proposed to the
Committee of Next Generation Lossless Compression of Continuoustone Still
Pictures
The results in Table demonstrate that the enhanced JBIG bitplane
coding scheme for continuoustone images performs the best among the four
algorithms tested In terms of the average bit rate the scheme outperforms
direct JBIG coding and the best mode of CREW by and bpixel
respectively and achieves much lower bit rates than the zerothorder entropies
of the prediction error images by an average of bpixel for the entire test
set of images If only the b component images are considered the enhanced
JBIG technique provides better compression performance by
and bpixel in terms of average bit rates when compared with lossless
JPEG coding direct Graytransform JBIG coding the zerothorder entropy of
the prediction error image and the best mode of lossless CREW respectively
For comparison with SLIC in the context of radiographic images the en
hanced JBIG PSVFJBIG procedure was applied for compression of the
ve mammograms and ve chest radiographs listed in Tables and
The bit rates achieved for the b and b versions of the images are listed in
Tables and respectively The results indicate that the enhanced
JBIG procedure provides an additional improvement over SLIC on the b
image set and that the method lowers the average bit rates to bpixel
from the bpixel rate of SLIC for the b images
Lowerlimit Analysis of Lossless Data Compres
sion
There is as yet no practical technique available for the determination of the
lowest limit of bit rate in reversible compression of a given image although
such a number should exist based upon information theory It is there
fore dicult to judge how good a compression algorithm is other than by
comparing its performance with those of other published methods or with the
zerothorder entropy of the decorrelated data if available the latter approach
is commonly used by researchers when dierent testimage sets are involved
and when other compression programs are not available Either of the two
approaches mentioned above can only analyze relatively how well a compres
sion algorithm performs in comparson with the other techniques available or
the zerothorder entropy It should be apparent from the results presented in
the preceding sections that the usefulness of comparative analysis is limited
For example SLIC provided the best compression results among the methods
TABLE
Compression of the JPEG Test Image Set Using Enhanced JBIG
EJBIG Lossless JPEG b Component Images Only Direct
Graytransformed JBIG Coding and CREW Coding BitsComponent
with the Best Bit Rate Highlighted in Each Case
Image Best Direct PSV Best
colsrowscompbits JPEG JBIG H
e
EJBIG CREW
hotel
gold
bike
woman
cafe
tools
bike
water
cats
aerial
aerial
cmpnd
cmpnd
nger
x ray
cr
ct
us
mri
faxball
graphic
chart
chart s
Average all
Average b comp only
Note H
e
zerothorder entropy of the PSV prediction error comp num
ber of components cols number of columns Reproduced with permission
from L Shen and RM Rangayyan Lossless compression of continuustone
images by combined interbitplane decorrelation and JBIG coding Journal
of Electronic Imaging
c
SPIE and IS%T
TABLE
Comparison of Enhanced JBIG PSVFJBIG or EJBIG with
JBIG JPEG Best Lossless Mode HINT and SLIC
errorbits Using Five b Mammograms m to m and m
and Five b Chest Radiographs c to c and c by
BitsPixel
Image bpixel H
JBIG JPEG HINT SLIC EJBIG
m
m
m
m
c
c
c
c
m
c
Average
The lowest bit rate is highlighted in each case See also Table Note
H
zerothorder entropy
TABLE
Comparison of Enhanced JBIG PSVFJBIG or EJBIG with
JBIG DBurgHuman DBH HINT and SLIC errorbits
Using Five b Mammograms m to m and m and Five b
Chest Radiographs c to c and c by BitsPixel
Image bpixel H
JBIG DBH HINT SLIC EJBIG
m
m
m
m
c
c
c
c
m
c
Average
The lowest bit rate is highlighted in each case See also Table Note
H
zerothorder entropy
tested in Section however it is seen in Section that the enhanced
JBIG scheme performs better than SLIC
Even if it were to be not practical to achieve the lowest possible bit rate
in the lossless compression of a given image it is of interest to estimate an
achievable lowerbound bit rate for an image Information theory indicates
that the lossless compression eciency is bounded by highorder entropy val
ues However the accuracy of estimating highorder statistics is limited by
the length of the data the number of samples available the number of in
tensity levels and the order The highest possible order of entropy that can
be estimated with high accuracy is limited due to the nite length of the
data available In spite of this limitation highorder entropy values can pro
vide a better estimate of the lowerbound bit rate than the commonly used
zerothorder entropy Shen and Rangayyan proposed methods for
the estimation of highorder entropy and the lowerbound bit rate which are
described in the following paragraphs
Memoryless entropy
A memoryless source is the simplest form of an information source in which
successive source symbols are statistically independent Such a source
is completely specied by its source alphabet A fa
a
a
a
n
g and
the associated probabilities of occurrence fpa
pa
pa
pa
n
g The
memoryless entropy HA which is also known as the average amount of
information per source symbol is dened as
HA
n
X
i
pa
i
log
pa
i
The m
th
order extension entropy is dened as
H
m
A H
m
a
i
a
i
a
i
a
i
m
X
A
m
p a
i
a
i
a
i
a
i
m
log
p a
i
a
i
a
i
a
i
m
where p a
i
a
i
a
i
a
i
m
is the probability of a symbol string from the
m
th
order extension of the memoryless source and A
m
represents the set of
all possible strings with m symbols following a
i
For a memoryless source
using the property
p a
i
a
i
a
i
a
i
m
pa
i
pa
i
pa
i
pa
i
m
we get
HA
H
m
A
m
Markov entropy
A memoryless source model could be restrictive in many applications due
to the fact that successive source symbols can be signicantly interdepen
dent which means the source has memory Image sources are such exam
ples in which there always exists some statistical dependence among neigh
boring pixels even after the source symbol stream has been decorrelated
A source possessing dependence or memory as above may be modeled as
a Markov source in which the probability of occurrence of a source sym
bol a
i
depends upon the probabilities of occurrence of a nite number m of
the preceding symbols The corresponding m
th
order Markov entropy
Ha
i
ja
i
a
i
a
i
m
may be computed as
H a
i
ja
i
a
i
a
i
m
X
A
m
p a
i
a
i
a
i
a
i
m
log
p a
i
ja
i
a
i
a
i
m
where p a
i
a
i
a
i
a
i
m
is the probability of a particular state and
p a
i
ja
i
a
i
a
i
m
is the PDF of a
i
conditioned upon the occurrence
of the string fa
i
a
i
a
i
m
g It may be shown that
H a
i
ja
i
a
i
a
i
m
H
a
i
ja
i
a
i
a
i
m
H a
i
ja
i
Ha
i
HA
where the equality is satised if and only if the source symbols are statistically
independent
Estimation of the true source entropy
Although a given image or its prediction error image could be modeled with
a Markov source the order of the model and the conditional PDFs will be
usually unknown However it is seen from Equation that the higher the
order of the conditional probability the lower is the resulting entropy which
is closer to the true source entropy In order to maintain a reasonable level
of accuracy in the estimation of the conditional probability larger numbers
of data strings are needed for higherorder functions the estimation error is
given by
mK
N ln
for an m
th
order Markov source model with
K
intensity levels and N data
samples Thus for a specic image with known K and N the highest order
of conditional Markov entropy that could be calculated within a practical
estimation error of conditional probability is bounded by
maxm b
log
N ln
K
c
where bxc is the !oor function that returns the largest integer less than or
equal to the argument x Therefore given an error limit the only way
to derive a higherorder parameter is to decrease K by splitting data bytes
because the data length N cannot be extended
For example the highest orders of Markov entropy that may be calculated
with a probability estimation error of less than for a bpixel
image are and for the original data one b data part K
split data with two b data parts K split data with four b data
parts K and split data with eight b data parts K respectively
Figure shows the Markov entropy values up to the maximum possible
order of maxm with for four forms of splitting for the test image
Airplane It is obvious that the entropy values become larger with splitting
into more data parts due to the high correlation present among the data parts
although the maximum order could go higher after splitting and the entropy
values are decreasing with increasing order for each form of splitting see also
Table This indicates that decorrelation of the data bits is needed before
splitting in order to get a good estimate of the entropy because the source
entropy is not changed by any reversible transformation
From the results of the enhanced JBIG algorithm it is seen that the Gray
F
and F
transformations provide good decorrelation among bit planes af
ter PSVbased prediction This is demonstrated with four plots of the bi
nary Gray F
and F
representations of PSV prediction error data of
the Airplane image in Figure as well as in Table The binary
representation is seen to lead to poor decorrelation among the bits of the
prediction error data it actually makes the maximumorder entropy increase
to bpixel while splitting the error data into eight b data parts the
zerothorder entropies of the original error data and the original image are
and bpixel respectively It is also seen that the highestorder
entropy values that could be estimated within the error limit specied are
increasing with increasing number of data parts when using the binary rep
resentation On the other hand with the Gray F
or F
transformation
the situation is dierent It is seen that the highestorder entropy values with
splitting become smaller than the entropy of the original data when one of the
three transformations is used which shows their ecient decorrelation eect
among the bit planes of the prediction error image Finally it is seen that
the F
transform is the best among the four representation schemes with the
lowest estimated entropy value of bpixel achieved when the prediction
error is split into four b data parts
The lowest estimated entropy value is not guaranteed to occur when the
prediction error is split into four b data parts The F
transform was ob
served to always provide the best or nearbest performance Table lists
the lowest estimated Markov entropy values with PSV prediction for the
testimage set used together with their bit rates obtained with the enhanced
JBIG scheme PSVFJBIG and the zerothorder entropies of the PSV
prediction error images It is seen that the higherorder entropy values pro
FIGURE
Markov entropy values up to the maximum order possible with error limit
for four forms of splitting for the bpixel Airplane image see
also Table Note Order indicates memoryless entropy Reproduced
with permission from L Shen and RM Rangayyan Lossless compression
of continuustone images by combined interbitplane decorrelation and JBIG
coding Journal of Electronic Imaging
c
SPIE and
IS%T
TABLE
Estimated Values of the Highest Possible Order of Markov Entropy
bpixel for the Airplane Image
Data Code One part Two parts Four parts Eight parts
transform b bpart bpart bpart
Original Binary
Airplane Gray
image
F
F
PSV Binary
prediction Gray
error of F
Airplane F
Values are shown with and without prediction combined with four dierent
code representation transformation schemes and with error limit
Reproduced with permission from L Shen and RM Rangayyan Lossless
compression of continuustone images by combined interbitplane decorrela
tion and JBIG coding Journal of Electronic Imaging
c
SPIE and IS%T
Figure a
Figure b
Figure c
vide lower estimates of the bitrate limit than the zerothorder entropies The
average entropy value decreases to from bpixel with higherorder
entropy estimation while using binary representation by using the F
trans
form instead of binary representation of the error data the average Markov
entropy value further reduces to bpixel
The disadvantage of using the zerothorder entropy to measure the perfor
mance of a data compression algorithm is clearly shown in Table the
enhanced JBIG PSVFJBIG coding scheme achieves an average bit rate
of bpixel compared with the average zerothorder entropy of bpixel
for the prediction error images Considering the higherorder entropy values
shown it appears that the compression eciency of the enhanced JBIG tech
nique could be further improved An important application of highorder
entropy estimation could be to provide a potentially achievable lower bound
of bit rate for an original or decorrelated image if the highorder entropy is
estimated with adequate accuracy
d
FIGURE
Plots of the Markov entropy values up to the maximum order possible with
error limit with four forms of splitting for PSV prediction error
of the bpixel Airplane image with a Binary representation
b Gray representation c F
transformation and d F
transformation
Note Order indicates memoryless entropy Reproduced with permission
from L Shen and RM Rangayyan Lossless compression of continuustone
images by combined interbitplane decorrelation and JBIG coding Journal
of Electronic Imaging
c
SPIE and IS%T
TABLE
Lowest Estimated Markov Entropy Values with PSV Prediction
for Eight b Test Images
JPEG with PSV
b image Bit rate Lowest entropy
columns rows H
e
EJBIG F
Binary
Airplane
Baboon
Cameraman
Lenna
Lenna
Peppers
Sailboat
Tiany
Average
Also shown are bit rates via enhanced JBIG bitplane coding of F
transformed PSV prediction error PSVFJBIG or EJBIG and the
zerothorder entropies H
e
of the PSV prediction error images in
bpixel Reproduced with permission from L Shen and RM Rangayyan
Lossless compression of continuustone images by combined interbitplane
decorrelation and JBIG coding Journal of Electronic Imaging
c
SPIE and IS%T
Application Teleradiology
Teleradiology is commonly dened as the practice of radiology at a distance
Teleradiology oers a technological
approach to the problem of eliminating the delay in securing the consultation
of a radiologist to patients in rural and remote areas The timely availabil
ity of radiological diagnosis via telecommunication could potentially reduce
the morbidity mortality and costs of transportation to tertiary healthcare
centers of patients in remotely situated areas and in certain situations in
developing countries as well
In the military environment a study in indicated that over
of the medical facilities with radiographic equipment had no radiologists as
signed to them and an additional had only one radiologist In such cases
teleradiology could be a vehicle for redistributing the imagereading workload
from understaed sites to more adequately staed central locations
According to a study conducted in the province of Alberta Canada
had a total of healthcare centers with radiological imaging facilities out
of which only had resident radiologists Sixtyone of the other cen
ters depended upon visiting radiologists The remaining centers used to
send their radiographs to other centers for interpretation with a delay of
days in receiving the results The situation was comparable in
the neighboring provinces of Saskatchewan and Manitoba and it was observed
that the three provinces could benet signicantly from teleradiology Even
in the case of areas served by contract radiologists teleradiology can per
mit evaluation and consultation by other radiologists at tertiary healthcare
centers in emergency situations as well as in complicated cases
Early attempts at teleradiology systems consisted of analog
transmission of slowscan TV signals over existing telephone lines ultrahigh
frequency UHF radiolinks and other such analog channels Analog
transmission and the concomitant slow transmission rates were satisfactory for
lowresolution images such as nuclear medicine images However the trans
mission times were prohibitively high for highresolution images such as chest
radiographs Furthermore the quality of the images received via analog trans
mission is a function of the distance which could result in an unpredictable
performance of radiologists with the received images Thus the natural pro
gression of teleradiology systems was toward digital transmission The initial
choice of the transmission medium was the ordinary telephone line operating
at bps bits per second Several commercial teleradiology systems
were based upon the use of telephone lines for data transmission Improve
ments in modem technology allowing transmission speeds of up to Kbps
over standard telephone lines and the establishment of a number of Kbps
lines for commercial use by telephone companies made this medium viable for
lowresolution images
The major reason for users reluctance in accepting early teleradiology sys
tems was the inability to meet the resolution of the original lm Spatial reso
lution of even pixels was found to be inadequate to capture the
submillimeter features found in chest radiographs and mammograms
It was recommended that spatial resolution of the order of
pixels with at least shades of gray would be required to capture ac
curately the diagnostic information on radiographic images of the chest and
breast This demand led to the development of highresolution laser digitiz
ers capable of digitizing Xray lms to images of the order of
pixels with bpixel by the mid s Imaging equipment capable of di
rect digital acquisition of radiographic images to the same resolution as above
were also developed in the late s Teleradiology system designers were
then faced with the problem of dealing with the immense amount of data in
volved in such highresolution digital images The transmission of such large
amounts of data over ordinary telephone lines involved large delays which
could be overcome to some extent by using parallel lines for increased data
transfer rate The use of satellite channels was also an option to speed
up image data transmission but problems associated with image data
management and archival hindered the anticipated widespread acceptance of
highresolution teleradiology systems Such diculties motivated advanced
research into image data compression and encoding techniques
The development of PACS and teleradiology systems share some historical
common ground Although delayed beyond initial predictions both PACS
and teleradiology established their presence and value in clinical practice by
the late s The following paragraphs provide a historical review of telera
diology
Analog teleradiology
The rst instance of transmitting pictorial information for medical diagnosis
dates back to when GershonCohen and Cooley used telephone lines
and a facsimile system adapted to convert medical images into video signals
for transmitting images between two hospitals km apart in Philadelphia
PA In a pioneering project in Jutras conducted what is
perhaps the rst teleradiology trial by interlinking two hospitals km apart
in Montreal Qu$ebec Canada using a coaxial cable to transmit tele!uoroscopy
examinations The potential of teleradiology in the provision of the services
of a radiologist to remotely situated areas and in the redistribution of radiol
ogists workload from understaed centers to more adequately staed centers
was immediately recognized and a number of clinical evaluation projects were
conducted Most of the
early attempts consisted of analog transmission of medical images via standard
telephone lines dedicated coaxial cables UHF radio microwave and satellite
channels with display on TV monitors at the receiving terminal James et
al give a review of the results of the early experiments Andrus and
Bird describe the concept of a teleradiology system in which the radi
ologist stationed at a medical center controls a video camera to zoom in on
selected areas of interest of an image at another site located far away and
observes the results in real time on a TV screen Steckel conducted
experiments with a system using an line closedcircuit TV system for
transmitting radiographic images within a hospital for educational purposes
and concluded that the system s utility far outweighed disadvantages such as
the inability to view a sequence of images belonging to a single study
In Webber and Corbus used existing telephone lines and slow
scan TV for transmitting radiographs and nuclear medicine images The
resolution achieved was satisfactory for nuclear medicine images but both the
spatial resolution and the grayscale dynamic range radiometric resolution
were found to be inadequate for radiographs A similar experiment using
telephone lines and slowscan TV by Jelasco et al resulted in
correct interpretation of radiographs Other experiments with slowscan TV
over telephone lines demonstrated the inadequacy of this medium and
also that the diagnostic accuracy with such systems varied with the nature of
the images
Webber et al used UHF radio transmission in for transmitting
nuclear medicine images and radiographs While the system worked satisfac
torily for nuclear medicine images evaluation of chest Xray images needed
zoom and contrast manipulation of the TV monitor Murphy et al used
a microwave link for the transmission of images of chest radiographs acquired
with a remotely controlled video camera over a distance of about km and
indicated that it would be an acceptable method for providing health care to
people in remote areas
Andrus et al transmitted Xray images of the abdomen chest bone
and skull over a km round loop using a MHz line TV channel in
cluding three repeater stations The TV camera was remotely controlled using
push buttons and a joystick to vary the zoom aperture focus and direction
of the camera It was concluded that the TV interpretations were of accept
able accuracy Such realtime operation calls for special skills on the part
of the radiologist requires coordination between the operator at the image
acquisition site and the radiologist at the receiving center and could take up
a considerable amount of the radiologist s valuable time Moreover practical
microwave links exist only between and within major cities and cannot serve
the communication needs of teleradiology terminals in rural and remote areas
In addition the operating costs over the duration of interactive manipulations
could reach high levels and render such a scheme uneconomical
In Lester et al used a satellite ATS for analog transmis
sion of videotaped radiologic information and concluded that satisfactory
radiographic transmission is possible if a satisfactory sensor of radiographic
images were constructed In Carey et al reported on the results
of an analog teleradiology experiment using the Hermes spacecraft They re
ported the eectiveness of TV !uoroscopy to be of that with conventional
procedures Page et al used a twoway analog TV network with the
Canadian satellite ANIKB to transmit radiographic images from northern
Qu$ebec to Montr$eal and reported an initial accuracy in TV interpretation
of with respect to lm reading The accuracy rose to after a
month training of the participant radiologists in the use of the TV system
The noise associated with analog transmission the low resolution of the TV
monitors used and the requirement on the part of the radiologists to partici
pate in realtime control of the imageacquisition cameras made the concept
of TV transmission of radiographic images unacceptable Furthermore the
noise associated with analog transmission is dependent upon the distance
Not surprisingly James et al reported that their teleradiology system
transmitting emergency department radiographs via a satellite channel from a
local TV studio was unacceptable due to a decrease in the accuracy of image
interpretation to about with respect to that with standard protocols
Digital teleradiology
Given the advantages of digital communication over analog methods
the natural progression of teleradiology was toward the use of digital data
transmission techniques The advent of a number of digital medical imaging
modalities facilitated this trend Digital imaging also allowed
for image processing enhancement contrast scaling and !exible manipula
tion of images on the display monitors after acquisition Many of the initial
attempts at digital teleradiology were
based upon microcomputers and used lowresolution digitization display and
printers The resolution was of the order of to pixels with
shades of gray mostly because of the unavailability of highresolution
equipment Gayler et al described a laboratory evaluation of such a
microcomputerbased teleradiology system based upon a b for
mat for image acquisition and display and evaluated radiologists performance
with routine radiographs They found the diagnostic performance to be signif
icantly worse than that using the original lm radiographs Nevertheless they
concluded that microcomputerbased teleradiology systems warrant further
evaluation in a clinical environment
In Rasmussen et al compared the performance of radiologists
with images transmitted by analog and digital means and lightbox viewing
of the original lms The resolution of digitization used was pixels
with bpixel The digital images were converted to analog signals for analog
transmission It was concluded that the resolution used would provide satis
factory radiographic images for gross pathological disorders but that subtle
features would require higher resolution
Gitlin Curtis et al and Skinner et al followed the
laboratory evaluation of Gayler et al with eld trials using standard
telephone lines at bps for the transmission of b images
from ve medicalcare facilities to a central hospital in Maryland A relative
accuracy of with videoimage readings was reported as compared
to standard lm interpretation This was a substantially higher accuracy than
that obtained in a preceding laboratory study the improvement was
attributed to the larger percentage of normal images used in the eld trial
and to the higher experience of the analysts in clinical radiology
In a eld trial in Gitlin used a matrix of pixels
bps telephone lines and lossy data compression to bring down the
transmission times A relative accuracy with video interpretation of with
respect to standard lm interpretation was observed The relative accuracy
was observed to be dependent upon the type of data compression used among
other factors
Gordon et al presented an analysis of a number of scenarios and
tradeos for practical implementation of digital teleradiology In related pa
pers Rangayyan and Gordon and Rangaraj and Gordon dis
cussed the potential for providing advanced imaging services such as CT
through teleradiology
In DiSantis et al digitized excretory urographs to
matrices and transmitted the images over standard telephone lines after data
compression to a receiving unit approximately km away A panel of three
radiologists interpreted the images on video monitors and the results were
compared with the original lm readings performed about a week earlier
An agreement of was found between the lm and video readings in the
diagnosis of obstructions However only of urethral calculi detected
with the original radiographs were also detected with the video images This
result demonstrated clearly that whereas a resolution of pixels
could be adequate for certain types of diagnosis higher resolution is required
for capturing all of the diagnostic information that could be present on the
original lm
In Kagetsu et al reported on the performance of a commercially
available teleradiology system using b images and transmission
over bps standard telephone lines after data compression Exper
iments were conducted with a wide variety of radiographs over a fourmonth
period An overall relative accuracy of was reported between the re
ceived images on video display and the original lms Based on these results
Kagetsu et al recommended a review of the original lms at some later date
because of the superior spatial and contrast resolution of lm
Several commerical systems were released for digital teleradiology in the
late s Although such systems were adequate for handling lowresolution
images in CT MRI and nuclear medicine they were not suitable for han
dling largeformat images such as chest radiographs and mammograms Ex
periments with such systems demonstrated the inadequacy of lowresolution
digital teleradiology systems as an alternative to the physical transportation
of the lms or the patients to centers with radiological diagnostic facilities Al
though higher resolution was required in the digitized images the substantial
increase in the related volume of data and the associated diculties remained
a serious concern Furthermore the use of lossy data compression schemes to
remain within the datarate limitation of telephone lines and other lowspeed
communication channels was observed to be unacceptable
Highresolution digital teleradiology
The development of highresolution image digitizers and display equipment
and the routine utilization of highdatarate communication media paved the
way for highresolution digital teleradiology In Carey et al
reported on the performance of the DTR teleradiology system from
DuPont consisting of a pixel laser digitizer with quan
tization levels a T satellite transmission channel at Mbps and a
DuPont laser lm recorder with possible shades of gray A nonlinear
mapping was performed from the original quantization levels to
levels on the lm copy to make use of the fact that the eye is more sensitive to
contrast variations at lower density With this mapping at the lower end of
the gray scale small dierences in gray values correspond to larger dierences
in optical densities than at the higher end of the gray scale Thus the over
all optical density range of the lm is much larger than can be obtained by
linear mapping Carey et al transmitted radiographic and ultrasono
graphic images over the system from Seaforth to London in Ontario Canada
and reported a relative accuracy of in reading the lasersensitive lm as
compared to the original lm It was concluded that the lasersensitive lm
clearly duplicated the original lm ndings However they also reported
contouring on the lasersensitive lm which might have been due to the
nonlinear mapping of the original gray levels to levels on the lm
Certain portions of the original gray scale with rapidly changing gray levels
could have been mapped into the same optical density on the lm giving rise
to contouring artifacts
Barnes et al suggested that the challenge of integrating the increas
ing number of medical imaging technologies could be met by networked mul
timodality imaging workstations Cox et al compared images digitized
to b and displayed on monitors with b pix
els digital laser lm prints and conventional lm They reported signicant
dierences in the performance of the three display formats digital hard copy
performed as well as or better than conventional lm whereas the interactive
display failed to match the performance of the other two They suggested
that although the dierences could be eliminated by training the personnel
in reading from displays and by using image enhancement techniques it was
premature to conclude either way
In Batnitzky et al conducted an assessment of the then
available technologies for lm digitization display generation of hard copy
and data communication for application in teleradiology systems They con
cluded that b laser digitizers displays with scan lines of
bpixel hard copiers that interpolate
matrices to matrices and the merger of computer and com
munication technologies resulting in !exible widearea networks had paved
the way for the acceptance of nal interpretation teleradiology completely
eliminating the need to go back to the original lms Gillespy et al
described the installation of a DuPont Clinical Review System consisting of
a laser lm digitizer with b pixels and a b
display unit and reported that clinicians were generally satised with the
unit Several studies on the contrast and resolution of highresolution digi
tizers demonstrated that the resolution of the original lm
was maintained in the digitized images
Several systems are now available for digital teleradiology including high
resolution laser digitizers that can provide images of the order of
b pixels with a spatial resolution of m or better highluminance
monitors that can display up to pixels at bpixel and non
interlaced refresh rates of fps and laserlm recorders with a spatial
resolution of m that can print images of size b pixels
Satellite cable or beroptic transmission equipment and channels may be
leased with transmission rates of several Mbps However the large amount
of data related to highresolution images can create huge demands in data
transmission and archival capacity Lossless data compression techniques can
bring down the amount of data and have a signicant impact on the practical
implementation of teleradiology and related technologies
The introduction of data compression encoding and decoding in digital
teleradiology systems raises questions on the overall throughput of the sys
tem in transmission and reception storage and retrieval of image data patient
condentiality and information security The compression of image data re
moves the inherent redundancy in images and makes the data more sensitive
to errors In dedicated communication links appropriate error con
trol should be provided for detecting and correcting such errors In the case
of packetswitched communication links the removal of redundancy by data
compression could result in increased retransmission overhead However with
sophisticated digital communication links operating at typical biterror rates
of in
and channel utilization throughput eciency of about using
highlevel packetswitched protocols the advantages of data compres
sion far outweigh the overheads due to the reasons mentioned above
Highresolution digital teleradiology is now feasible without any sacrice in
image quality and can serve as an alternative to transporting patients or lms
Distance should no longer be a limitation in providing reliable diagnostic
service by citybased expert radiologists to patients in remote or rural areas
Remarks
Lossless data compression is desirable in medical image archival and transmis
sion In this chapter we studied several lossless data compression techniques
A lossless compression scheme generally consists of two steps decorrelation
and encoding The success of a lossless compression method is mainly based
upon the eciency of the decorrelation procedure used In practice a decor
relation procedure could include several cascaded decorrelation blocks each of
which could accomplish a dierent type of decorrelation and facilitate further
decorrelation by the subsequent blocks Some of the methods described in
this chapter illustrate creative ways of combining multiple decorrelators with
dierent characteristics for achieving better compression eciency
Several informationtheoretic concepts and criteria as applicable to data
compression were also discussed in this chapter A practical method for the
estimation of highorder entropy was presented which could aid in the lower
limit analysis of lossless data compression Highorder entropy estimation
could aid in the design analysis and evaluation of cascaded decorrelators
A historical review of selected works in the development of PACS and telera
diology systems was presented in the concluding section in order to demon
strate the need for image compression and data transmission in a practical
medical application PACS teleradiology and telemedicine are now well es
tablished areas that are providing advanced technology for improved health
care
Study Questions and Problems
Selected data les related to some of the problems and exercises are available at the
site
wwwenelucalgarycaPeopleRangaenel
The probabilities of occurrence of eight symbols are given as
Derive the Human code for the source
For a symbol source derive all possible sets of the Human code The
exact PDF of the source is not relevant
Create a few strings of symbols and generate the corresponding Human
codes Verify that the result satises the basic requirements of a code in
cluding unique and instantaneous decodability
For the image shown below design a Human coding scheme Show all
steps of your design
Compute the entropy of the image and the average bit rates using direct
binary coding and Human coding
Discuss the similarities and dierences between the KarhunenLoeve discrete
Fourier and the WalshHadamard transforms Discuss the particular aspects
of each transform that are of importance in transformdomain coding for image
data compression
For the image shown below compute the bit rate using
a direct binary coding
b horizontal runlength coding and
c predictive coding or DPCM using the model
fmn fmn
Show and explain all steps State your assumptions if any and explain your
procedures
For the image shown below prepare the bit planes using the direct
binary and Gray codes Examine the bit planes for the application of run
length coding
Which code can provide better compression Explain your observations and
results
Laboratory Exercises and Projects
Write a program in C C or MATLAB to compute the histogram and
the zerothorder entropy of a given image Apply the program to a few images
in your collection Study the nature of the histograms and relate their charac
teristics as well as entropy to the visual features present in the corresponding
images
Write a program to compute the zerothorder and rstorder entropy of an
image considering pairs of graylevel values Apply the program to a few
images in your collection and analyze the trends in the zerothorder and rst
order entropy values
What are the considerations complexities and limitations involved in com
puting entropy of higher orders
For the image in the le RajREyedat with bpixel create bit planes using
a the binary code and b the Gray code Compute the entropy of each bit
plane Compute the average entropy over all of the bit planes for each code
What is the expected trend
Do your results meet your expectations Explain
Develop a program to compute the DFT of an image Write steps to compute
the energy contained in concentric circles or squares of several sizes spanning
the full spectrum of the image and to plot the results
Apply the program to a few images in your collection Relate the spectral
energy distribution to the visual characteristics of the corresponding images
Discuss the relevance of your ndings in data compression
Develop a program to compute the error of prediction based upon a few simple
predictors such as
a
fmn fm n
b
fmn fmn
c
fmn fm n fmn fm n
Derive the histograms and the entropies of the original image and the error of
prediction for a few images with each of the predictors listed above Evaluate
the results and comment upon the relevance of your ndings in image coding
and data compression