IEEE
TRANSACTIONS
ON
BIOMEDICAL
ENGINEERING,
VOL.
BME28,
NO.
2,
FEBRUARY
1981
tomography,
in
Statistical
Mechanics
and
Statistical
Methods
in
Theory
and
Application,U.
Landman,
Ed.
New
York:
Plenum,
1977.David
Nahamoo
was
born
in
Hamedan,
Iran,
in
1953.
HeDreceived
the
B.Sc.
degree
in
electrical
engineering
from
Tehran
University,
Iran,
and
the
M.Sc.
degree
in
electrical
engineering
from
Imperial
College
of
Science
and
Technology,
University
of
London,
England,
in
1975and
1976,
respectively.
He
is
currently
working
toward
a
Ph.D.
degree
in
electrical
engineering
at
Purdue
University,
West
Lafayette,
IN.
He
has
been
a
Research
Assistant
at
Purdue
University
since
August,
1977.
Mr.
Nahamoo
is
a
member
ofPhi
Kappa
Phi.
Carl
R.
Crawford
was
born
in
Milwaukee,
WI,
in
1956.
He
received
the
B.S.
and
M.S.
degrees
in
electrical
engineering
fromPurdue
University,
West
Lafayette,
IN,
in
1976
and
1977,
respec
14
A
dz
flyl
ively.
He
is
currently
pursuing
a
Ph.D.
degree
in
electrical
engineering
at
Purdue
University.
Mr.
Crawford
is
a
member
of
Eta
Kappa
Nu
and
Tau
Beta
Pi.
Avinash
C.
Kak
(M'7
1),
for
a
photograph
andbiography,
see
this
issue,
p.49.
RapidExecution
of
Fan
Beam
Image
Reconstruction
Algorithms
Using
Efficient
Computational
Techniques
and
SpecialPurpose
Processors
BARRY
K.
GILBERT,MEMBER,
IEEE,
SURENDER
K.
KENUE,
RICHARD
A.
ROBB,
ALOYSIUS
CHU,
ARNOLD
H.
LENT,
AND
EARL
E.
SWARTZLANDER,
JR.,
SENIOR
MEMBER,
IEEE
AbstractRapid
advances
during
the
pastten
years
of
several
forms
of
computerassisted
tomography(CT)
have
resulted
in
the
develop
ment
of
numerous
algorithms
to
convert
raw
projection
data
into
cross
sectional
images.
These
reconstruction
algorithms
are
either
iterative,
in
which
a
large
matrix
algebraic
equation
is
solved
by
successive
ap
proximation
techniques;
or
closed
form.
Continuing
evolution
of
the
closed
form
algorithms
hasallowed
thenewest
versions
to
produce
excellent
reconstructed
images
in
most
applications.
This
paper
will
review
several
computer
software
and
specialpurpose
digital
hardwareimplementations
of
closed
form
algorithms,
either
proposed
during
the
past
several
years
by
a
number
of
workers
or
actually
implemented
in
commercial
or
research
CT
scanners.
The
discussion
willalso
cover
a
number
of
recentlyinvestigated
algorithmic
modifications
which
reduce
Manuscript
received
March
2,
1980;
revised
September
1,
1980.
This
work
was
supported
in
part
by
the
U.S.
Public
Health
Service
under
Grant
HL04664,
the
National
Institutes
of
Health
underGrant
RR
00007,
the
U.S.
Air
Force
under
Grant
F3361579C1875,
andthge
Fannie
E.
Rippel
Foundationunder
a
research
grant.
B.
K.
Gilbert,
R.
A.
Robb,
and
A.
Chu
are
with
the
Biodynamics
Re
search
Unit,
Department
of
Physiology
and
Biophysics,
Mayo
Founda
tion,
Rochester,
MN
55901.
S.
K.
Kenue
was
with
the
Biodynamics
Research
Unit,
Department
of
Physiology
and
Biophysics,
Mayo
Foundation,
Rochester,
MN
55901.
He
is
now
with
the
Picker
Corp.,
Cleveland,
OH.
A.H.
Lent
was
with
the
Biodynamics
Research
Unit,
Department
of
Physiology
and
Biophysics,
Mayo
Foundation,
Rochester,
MN
55901.
He
is
now
with
the
New
Ventures
Division,
Technicare
Corp.,Salem,
OH.
E.E.
Swartzlander,
Jr.
is
with
TRW
Defense
andSpace
Systems
Group,
Huntsville,
AL.
the
amount
of
computation
required
to
execute
the
reconstruction
process,as
well
as
several
new
specialpurpose
digital
hardware
imple
mentations
under
development
in
our
laboratories
at
the
Mayo
Clinic.
INTRODUCTION
D
APID
advancesduring
thepast
ten
years
of
several
forms
^
of
computer
assisted
tomography
(CT)
have
resulted
in
the
development
of
numerous
algorithms
to
convert
raw
projection
data
into
crosssectional
images.
These
reconstruc
tion
algorithms
areeither
iterative,
in
which
alarge
matrix
algebraic
equation
is
solved
by
successive
approximation
tech
niques;or
closed
form.
Closed
form
reconstruction
algo
rithms
which
operate
in
the
frequency
domain
are
generally
referred
to
as
Fourieralgorithms;
those
which
operate
in
the
spatial
domain
are
referred
to
as
filtered
back
projection
algorithms.
Though
the
iterative
algorithms
are
frequently
employed
when
the
signal/noise
ratio
of
the
inputdata
is
low,
or
when
several
of
theprojections
are
unavailable,
they
are
in
less
frequent
use
for
high
quality
data,
since
iterative
solutions
consume
large
amounts
of
computer
time.
Continuing
evolu
tion
of
the
closed
form
algorithms
has
allowed
the
newest
versionsto
produce
excellent
reconstructed
images
in
most
applications.
This
paper
will
review
several
computer
software
and
special
purpose
digital
hardware
implementations
of
closed
form
00189294/81/02000098$00.75
©
1981
IEEE
GILBERT
et
al.:
FAN
BEAM
IMAGE
RECONSTRUCTIONALGORITHMS
algorithms,
either
proposed
during
the
past
several
years
by
a
number
ofworkers
oractually
implemented
in
commercial
orresearch
CT
scanners.
The
discussion
will
also
cover
several
new
specialpurpose
digital
hardware
implementations
under
development
in
our
laboratories
at
the
Mayo
Clinic,
Rochester,
MN.
Alternate
implementations,
using
analog
electronics
or
optical
methods
[1],
though
of
considerable
interest,
have
been
omitted
from
considerationhere
in
the
interests
of
brevity.
For
excellentrecent
overviews
of
the
numerous
types
of
iterative
and
closed
form
reconstruction
algorithms,
see
[2][4].
Early
Xray
CT
scanners
used
a
pencil
beam
of
Xrays,
thereby
allowing
anassumption
that
all
Xrays
from
a
given
angle
of
view
were
parallel
to
one
another,
in
turn
allowing
so
called
parallel
beam
reconstruction
algorithms
to
be
imple
mented.
In
order
to
permit
more
rapid
scanning
and
datacollection
than
was
feasible
for
pencil
beam
systems,
the
newer
Xray
CT
scanners
now
employ
a
fan
beam
of
Xrays
with
an
edge
to
edge
beam
divergence
of
approximately
30°450.
Though
the
fan
beam
algorithms
were
more
diffi
cult
to
derive
and
frequently
are
more
costly
to
execute
than
the
parallel
beam
algorithms,
they
are
today
in
common
use.
Implementation
of
any
of
the
closed
form
reconstruction
methods
requires
a
careful
scrutiny
ofevery
stage
of
the
process,
including
an
analysis
of
thealgorithm's
numerical
char
acteristics,
and
selection
of
the
appropriateprocessor
on
which
to
execute
it.
Section
I
of
this
paper
will
concentrate
on
the
algorithmic
aspects
of
such
an
implementation.
Although
a
brief
description
will
be
presented
of
various
fan
beam
algo
rithmscurrently
extant,
primary
emphasis
will
be
given
to
the
filtered
back
projection
approach.
Methods
of
executing
the
filtration
portion
of
this
algorithm,
as
well
as
simplifications
thereof
to
decreasethe
required
amount
ofcomputation,
will
be
discussed.
The
requirement
for,
and
mechanization
of,
an
interpolation
of
extravalues
into
the
filtered
projections
will
then
be
described.Several
optimizations
of
the
back
projec
tion
computations
using
approximation
methods
will
be
covered,
followed
by
an
example
of
improvementswhich
can
be
identified
in
the
fundamental
numerical
operations
employed
in
reconstruction,
e.g.,
the
inner
product.Section
II
of
the
paper
concentrates
on
physical
implemen
tation
of
the
algorithms,
whether
in
software
or
on
specially
designed
processors.
A
discussion
of
several
implementations
of
the
filtration
step
is
followed
by
descriptions
both
of
serial
processing
( projectiondriven )
and
of
parallel
processing
( pixeldriven )
approaches.
Finally,the
impactof
advanced
technology
on
processorefficiency
will
beexamined.
I.
ALGORITHMIC
CONSIDERATIONS
A.
Fan
Beam
Reconstruction
byRebinningThe
first
attempt
to
generalize
the
parallel
beam
filtered
back
projection
reconstruction
algorithm
of
Ramachandran
and
Lakshminaryanan
[5]
to
a
fan
beam
energy
source
is
based
upon
therecognition
that
projections
generated
by
a
fan
beam
source
and
collected
from
around
3600
of
the
object
to
be
reconstructed
contain
data
samples,
or
ray
sums,
which
could
equally
well
have
been
formed
by
energy
sourcesgenerat
ing
parallel
beams
of
penetrating
radiation.
Judicious
rearrange
ment
of
ray
sums
from
all
fan
beam
generated
projections
could
thus
create
a
new
set
of
projections
with
the
same
characteristics
as
thosegenerated
by
parallelsets
of
penetrat
ing
rays.
The
rearranged
projections
couldthen
be
reconstructed
into
crosssectional
images
by
application
of
parallel
beam
reconstruction
algorithms
such
as
thatderived
in
[5].
These
rebinning
algorithms
[6]
[9]
have
been
implemented
in
software
by
several
workers
and
also
in
hardware
for
at
least
two
specialpurpose
computers
dedicated
to
thereconstructiontask
[10].
One
of
the
major
deficiencies
of
rebinningalgorithms
is
the
requirement
to
execute
a
twodimensional
interpolation
of
ray
sums
in
both
the
angular
dimension
andbetween
measured
ray
sums.
If
no
interpolation
is
carried
out,
rebinning
algo
rithms
generate
parallel
projections
in
which
the
number
of
ray
sums
is
usually
very
much
less
than
the
number
of
srcinal
ray
sums
recorded
in
thefan
beam
projections.
Thisreduction
in
sample
density
results
in
a
decrease
in
bandwidth
of
the
information
contentof
the
parallel
projections,
and
is
thusnever
used
in
actual
practice.
In
some
schemes,
additional
projections
are
generated
by
linear
interpolation
of
projections
between
two
consecutivefanangles;thedesired
parallel
ray
sums
are
thenobtained
by
interpolation
from
the
two
closest
fan
beam
ray
sums.
Linear
interpolation
between
measured
ray
sums
usually
suffices,
butdoes
degrade
the
resolution
of
the parallelized projection
data.
Higher
orderinterpolation
(e.g.,
cubic
splines)
canbeused
at
the
expenseof
morecompu
tation.
In
addition,
the
filtration
and
back
projection
(see
SectionIC
for
a
discussion
of
back
projection)
of
the paral
lelized
projections
cannotproceed
until
all
projections
have
been
collected.
However,
the
rebinningalgorithms
are
more
computationally
efficient
than
the
direct
fan
beam
algorithmsduring
the
back
projection
portion
of
theprocedure,
since
they
omit
a
weighting
operation
which
must
beexecuted
in
the
direct
fan
beam
algorithms
(as
described
later).
B.
Fan
Beam
Reconstruction
by
Fourier
Techniques
Several
workers
have
presented
parallel
beam
and
fan
beam
reconstruction
algorithms
which
operate
partially
or
totally
in
the
frequency
(Fourier)
domain.
In
onesuch
algorithm,
the
raw
projections
are
first
back
projected; the
back
projected
image
is
then
filtered
in
the
twodimensional
(2D)
Fourier
domain
by
a
suitably
chosen
filter
function,
e.g.,
a
ramp,
ora
ramp
multiplied
by
a
window
function.
Fourier
filtration
is
accomplished
by
taking
the
forward2D
fast
Fourier
trans
form
(FFT)
of
the
back
projected
image,
then
multiplying
each
frequency
component
of
the
transformedimage
by
the
filter
function,
and
finally
by
forming
the
inverse
2D
FFT
of
the
transformed
image
[11],[12].
Although
the
use
of
the
FFT
technique
tends
to
minimize
the
amountof
computation
required
to
execute
the
filtration
process(see
below),
this
method
requires
four
times
morecomputer
memory
than
for
the
spatial
domain
techniques.
C.
Truly
Fan
Beam
Reconstruction
Algorithms
The
development
in
1975
of
a
truly
fan
beam
filtered
back
projection
reconstruction
algorithm
[13],
followed
by
a
rigor
ous
derivation
thereof
[14]
[16],
resulted
in
a
conversion
bymany
workers
from
parallel
beam
or
rebinning
approaches
to
the
direct
reconstruction
of
thefan
beam
projection
data.
99
IEEE
TRANSACTIONS
ON
BIOMEDICAL
ENGINEERING,
VOL.
BME28,
NO.
2,
FEBRUARY
1981
This
algorithm
is
composed
of
three
stages.First,
the
unpro
cessed
raw
projection
vectors
are
multiplied
by
a
precalculated
weighting
vector
of
the
same
length.
Each
Xray
projection
is
then
filtered
with
a
carefully
formulated
finite
impulse
response
(FIR)
filter
kernel
[17],
either
by
direct
onedimensional
con
volution
in
the
spatial
domain,
or
by
means
of
the
FFT
or
numbertheoretic
transforms
[18],
[191
and
direct
vector
multiplication
of
the
data
vector
and
filter
function
vector
in
the
frequency
domain
[2],
followed
by
inverse
transformation
back
to
the
spatial
domain.
Afterpreweighting
and
filtration,
the
projections
are
back
projected
into
the
image
space
to
generatethegray
scale
values
of
the
image
picture
elements
(pixels).
To
generate
an
image
pixel
by
back
projection,
one
ray
sum
from
each
of
the
filtered
projections
is
selected
according
toa
prespecified
selection
algorithm;
all
ray
sum
values
are
added
together
tocreate
thegray
scale
value
of
the
pixel
[2]
.
Stated
alternately,
since
each
ray
sum
representsthe
sum
of
the
energy
absorp
tion
coefficients
of
the
pixels
through
which
its
associated
ray
traversed,
includingappropriate
corrections
for
theactual
ray
path
length
through
each
pixel,
back
projection
is
thepro
cedure
by
which
the
magnitude
ofeach
ray
sum
is
partitioned
among
those
pixels.
In
addition,
if
the
individualrays
in
the
planar
beam
diverge
in
their
passage
through
the
structure,
as
is
the
case
when
the
energy
source
radiates
in
a
fanshaped
pat
tern,
the
intensity
of
the
unattenuated
beam
decreases
with
increasing
distance
from
the
source.
Hence,
a
weighting
term
must
be
incorporated
into
the
calculation
of
the
contribution
of
each
ray
sum
to
its
associated
image
pixels
as
a
function
of
the
relative
locations
of
the
source
and
the
pixels.
The
back
projection
equations
for
datacollected
with
such
a
fanshaped
beam
are
formulated
so
that
the
individual
ray
which
passed
through
eachimage
pixel
at
a
given
angular
position
of
thesource
is
identified,
its
corresponding
ray
sum
is
indexed
from
the
array
of
all
measured
ray
sums,
and
interpolation
betweenmeasured
ray
sums
is
executed
if
a
measured
ray
did
not
traverse
the
pixel
in
question
(see
below).
The
appropriate
weighting
factor
to
correct
for
the
divergence
of
the
beam
is
then
calculated
and
applied
to
the
ray
sum
to
produce
thecontribution
to
the
pixel
from
that
source
position.
The
procedure
is
repeated
for
all
angular
positions
of
the
energy
source
until
the
ray
sums
from
all
angles
of
viewhave
been
back
projected
to
form
thegray
level
of
that
pixel;
theproce
dure
is
then
repeated
for
all
the
other
image
pixels.
Several
versions
of
direct
fan
beam
reconstruction
algorithms
have
appeared
during
the
past
few
years
[3],
[13],[14],
[16],
[20],[21].
The
equations
for
fan
beam
filtered
back
projection,
assum
ing
a
curved
Xray
detector
strip
[14][16]
,
are
Zji
=
cJ i
XAi
N
yj,,
=
E3
ali
Zjj
i=N
A~r,
0)

E,
Wj,
ro
YK(,rO)
j=i
preweighting
(1)
filtration
(2)
back
projection
(3)
{r
cos
(j3
4
1
K(j,
r,
))
=
tan'O
a
D
+
rsin
pi
4))
})
ray
sum
index
(4)
131+l

131
Wjr4=
2U2
back
U
=
[(r
cos
i3

0))2
}
projection
+
(D
+
r
sin
(pi
0
4))2
1
1/2
(5)
in
which
the
Xi,
are
samples
of
the
unprocessed
projection;
the
letter
c
is
a
constant;J(i)
is
a
preweighting
vector;
and
the
a1,I
are
the
elementsof
a
convolution
kernel
which
is
used
to
filter
all
projections
[each
projection
contains
(2N+
1)
ele
ments].
The
Z1,i
are
elements
of
a
preweighted
projection
prior
tofiltration,
and
Ym,1
are
the
elementsof
aprojection
after
filtering.
The
valuef(r,
4)
is
the
magnitude
of
a
back
projected
pixel
at
location
(r,
4)
(the
sign
of
0
is
positive
for
counterclockwise
angular
displacements
measured
from
the
positivexaxis),
and
the
Yj,K
represents
the
Kth
ray
sum
from
the
jth
filtered
projection;the
subscripts
on
Y,
i.e.,
[j,
K(j,
r,
4)],
representthe
indexof
ray
sum
Y,
or
theray
sum
index,
referred
to
the
image
pixel
at
location
(r,
4),
from
the
jth
projection.
Although
thesubscript
i
is
an
independent
variable,
K
must
be
computed
from
(4).
The
Wjir,,
represents
amultiplicative
weight,
applied
to
the
Kth
ray
sumfrom
the
jth
filtered
projection
undergoing
back
projection
into
the
image
pixel
at
location
(r,
0).
D
represents
the
energy
source
tocenterofobject
distance;
Piy
is
theangle
between
the
yaxis
and
a
line
connecting
the
jth
source
position
and
the
center
of
theobject;
and
ca
is
the
angular
increment
between
individualrays
in
each
fan
beam.
D.
Overview
of
Mechanizations
of
the
Filtration
Step
The
linearfiltration
of
the
individual
projections
can
be
carried
out
by
any
viable
filtration
method.
For
computed
tomography
applications,
direct
convolution
filtration
has
beenmost
frequently
employed,
followed
in
popularity
by
an
application
of
the
classical
FFT
to
implement
a
socalled
FFTfast
convolution
linear
filtration
algorithm
[22].
The
number
of
multiplication
and
addition
operationsrequired
for
Fourier
filtration
of
vectors
containing
2561024
or
more
ele
ments
is
considerably
smaller
than
for
direct
convolution
approaches.
The
relative
cost
of
arithmetic
executed
on
generalpurpose
computers
has
been
sogreat
that
thereduction
in
absolute
numbers
of
computations
achieved
by
the
FFTfast
convolution
method
has
resulted
in
a
significant
reduction
in
processing
cost.
In
addition,
although
the
FFT
is
capable
of
transforming
vectors
composed
of
complex
elements,
the
elements
of
the
projection
data
vectors
to
be
filtered
are
pure
real
(i.e.,
not
complex
elements).
In
such
a
case,
the
effective
throughputof
the
FFT
algorithm
canbe
improved
considerably
by
packing
two
real
vectors
together
and
transforming
them
simultaneously;
such
approaches
are
referred
to
as
decima
tion
in
time
or
decimation
in
frequency
[23],
[24].
None
theless,
several
characteristics
of
the
classical
FFT
algorithms
present
their
own
unique
sets
of
complications.
FFT
convolu
100
GILBERT
et
al.:
FAN
BEAM
IMAGERECONSTRUCTION
ALGORITHMS
tion,
even
of
real
vectors,
requires
executionof
numerical
operations
in
the
domain
of
complex
numbers;
further,
because
these
computations
require
use
of
trigonometric
func
tions
and
thereby
of
irrational
numbers,
the
availability
of
floating
point
arithmetic
software
and
hardware
capability
is
desirable
(though
not
absolutely
mandatory)
to
minimize
com
putational
roundoff
error.
Inaddition,to
ensurethat
aliasing
errors
[22]
do
not
appear
in
the
data
vector
during
the
trans
formation
or
inverse
transformation
processes,
a
set
of
zero
elements
must
be
included
at
each
end
of
the
vector
to
be
transformed,
thereby
effectively
increasing
the
bandwidth
of
the
transformoperation
itself.
A
variety
of
FFT
algorithms
are
documented
in
[24]
[27].
The
theory
of
linear
transformationshas
recently
been
expanded
by
the
development
of
numbertheoretic
and
finite
field
transforms
[18],
[191,
[28],
[29]
,
a
class
of
linear
trans
formations
defined
to
exist
only
on
carefully
selected
finite
fields
of
numbers,
e.g.,
the
subspace
of
real
rational
(integer)
values.
Although
these
new
transforms
require
a
much
lower
absolute
count
of
arithmeticoperations
even
in
comparison
to
the
classical
FFT,
they
appear
torequire
significantly
greater
numbers
of
logical
operations
and
memory
references,
as
well
as
very
complicated
control
sequences
to
assure
correct
genera
tion
of
intermediate
results.
Though
the
finite
field
transforms
appear
slightly
faster
than
the
FFT
when
executed
on
general
purpose
computers,
it
has
yet
to
be
demonstrated
conclusively
that
they
are
superior
to
the
classical
FFT
or
direct
convolu
tion
methods
in
computedtomography
applications.
Direct
convolution
of
real
vector
operands
of
fixed
precision
can
be
carried
out
entirely
within
the
domain
of
real
integers,
thereby
eliminatingthe
requirement
forfloating
point
pro
cessors.
The
conceptual
simplicity
of
direct
convolution,
and
the
straightforward
design,
preparation,
and
testing
of
appro
priate
software
or
hardware
makes
implementation
of
the
direct
convolution
approach
very
attractive.
In
addition,
direct
convolution
can
be
simplified
by
exploiting
specific
numerical
characteristics
of
the
filter
function
itself.
For
example,
a
considerable
computational
saving
can
be
achieved
if
the
filter
function
is
even,
i.e.,
symmetric
about
its
ordinate,
in
the
spatial
domain.
Inthat
case
halfthe
kernel
may
be
reflected
about
theordinate;
since
the
filter
weights
a1
and
ai
have
the
same
numerical
value,
elements
of
theprojectionvector
which
multiply
an1
and
ai
are
first
summed,
and
there
after
multiplied
by
ai.
Only
the
order
of
execution,
not
the
number,
of
additions
is
altered,
but
the
requisite
number
of
multiplications
is
nearly
halved.
That
this
approach
can
always
be
invoked
is
assured
by
the
requirement
that
a
linear
phase
response
be
maintained
by
the
filter
function
to
prevent
phase
distortion
of
the
filtered
projections,
and
thereby
of
the
complete
reconstructions.
Linear
phase
response
is
guaranteed
only
if
the
filter
function
is
pure
real
in
the
spatial
frequency
domain,
and
therefore
aneven
function,
i.e.,
symmetric
about
its
ordinate,
in
the
spatial
domain
[30].
Another
simplification
of
direct
convolution
filtration
is
based
on
the
extension
of
a
wellknown
signal
processingcon
cept,
that
of
the halfband
filter
kernel.
In
such
an
imple
mentation,
all
elements
in
the
filter
vector
with
even
indices,
except
the
central
(zeroth)element,
are
defined
to
be
zero
[31].
The
resulting
restriction
in
equivalent
filtration
band
width
of
thehalfband
function,
in
comparison
tothefullbandkernel,
does
not
preclude
satisfactory
image
quality
provided
thatthe
bandwidth
of
the
projection
data
is
not
too
great.
As
an
extension
of
this
concept,
it
is
alsopossible
to
conceptualize
a
filter
kernel
in
which
several
consecutive
elements
are
identicallyzero,
followed
by
one
or
more
consecutive
non
zero
elements,
and
so
on
[32][34],
yieldinga
fractional
band
filter
kernel
[35],with
aneven
more
restricted
filtration
bandwidth
(e.g.,
a
lowpass
filter
with
an
even
lower
cutoff
frequency).
The
tendency
ofsuch
fractionalband
filter
kernels
to
smooth
the
reconstructed
images
excessively
can
be
partially
overcome
by
judicious
selection
of
thevalues
of
the
remaining
nonzero
elements.
Algorithms
have
been
developed
recently
[32]
[34]
for
generating
fractionalbandkernels
which
appear
to
yield
satisfactory
reconstructed
images
from
high
resolution
phantoms
scanned
in
research
Xray
and
ultra
sound
CT
scanners.Direct
convolution
filtration
can
be
further
simplified
by
exploitingrecently
developed
techniques
[32],
[33],
[36]
which
can
convert
agiven
convolution
kernel
into
one
in
which
the
central
value
is
afixed
precision
irrational
number,
while
all
othervalues
of
the
kernel
are
simple
fixedprecision
binary
multiples
of
one
another,
i.e.,
values
such
as
2,
16,
64,
or

128.
The
conversion
ofnonbinary
filter
weights
to
fixed
precision
binary
numbersmakes
it
possible
to
execute
the
filtration
primarily
by
binary
shift
and
add
operations,
nearly
obviatingthe
requirement
for
execution
of
multiplications.
This
technique
is
so
powerful
when
implemented
on
binary
arithmetic
computers
that
it
will
be
discussed
in
greater
detail.
It
is
well
known
that
the
integral
of
a
continuous
convolu
tion
kernel
is
zero;
however,
when
such
a
kernel
is
converted
intoa
sampled
versionfor
digital
processing,
and
truncated
to
an
acceptable
(i.e.,
noninfmite)
length
to
allow
reasonable
execution
times,the
sum,
E,
of
the
resulting
kernel
can
no
longer
be
zero.
Such
nonzerosum
convolution
kernels
are
widely
used
in
computedtomography
applications
with
accept
able
results.
Satisfactory
image
qualityhas
been
achieved
by
emulating
the
characteristics
of
these
nonzerosum
kernels
during
the
formation
of
a
binary
valued
kernel
from
an
arbi
trary
convolution
kernel.
The
binary
kernel
is
generated
in
the
following
manner.
Let
S
be
a
binary
approximation
to
the
central
element
of
thetruncated
nonbinary
kernel.
The
truncated
nonbinary
kernel
is
divided
by
the
value
of
its
own
central
element
to
yield
a
normalized
kernel
K,,
(i).
The
binarykernel
KB
(i)
=2
(except
for
the
central
element
i
=
0)
is
obtained
from
K,,
(i)
such
that
lKn(i)

2ml
<
W~n
(i)
2ni,
m
=A
n
Iml

Inl
<
1
(6)
assuming
m
and
n
areintegers.
The
central
element
KB
(0)
is
formed
by
summing
the
other
binaryelements.
Every
element
is
then
multiplied
by
S,
which
is
merely
a
number
of
the
form
101
IEEE
TRANSACTIONS
ON
BIOMEDICAL
ENGINEERING,
VOL.
BME28,
NO.
2,
FEBRUARY
1981
2 .
Lastly,
the
central
element
is
adjusted
(incremented
or
decremented)
by
thevalue
E.
Thus,
except
for
thecentral
element,
all
other
elements
are
binary
numbers
of
the
form

2 .
Lastly,
acorrection
must
be
made
to
all
pixels
of
the
final
reconstructedimage,
consisting
of
a
constant
offset
to
correct
its
average
gray
level,
followed
by
a
multi
plicativescaling
of
the
entire
image.
It
is
of
interest
to
note
that
the
generation
of
a
binarykernelfrequently
results
in
several
adjacent
elements
with
the
same
numerical
value.
Thus,
a
furthersimplification
in
computation
can
be
achievedduring
the
execution
of
each
inner
product
phase
of
the
convolution
by
first
summing
the
data
values
corresponding
to
the
identical
kernel
elements,
and
then
by
executing
a
single
binary
shift
operation,
in
contradistinction
to
theusual
case
in
which
a
binary
shift
operation
is
executed
for
each
data
value.
Forexample,
theusual
convolution
formulation
for
the
projection
P(i)
and
kernel
K(i)
can
be
rewritten
in
a
somewhat
simplified
format
from
that
of
(2)
as
N
(P*K)=Z
K(nj)P(n).
(7)
n=o
When
executed
using
binary
kernels,
(7)
canbe
rewritten
as
L
(P
*
KB),
P(j)KB(O)
+
g;n
2n
(8)
n=l
Sn
g7
=
E
P(t)
(9)
t=Rn
where
Rn
and
Sn
are
the
lower
and
upper
indices
in
the
binarykernel
of
identical
terms
of
power
2n,
and
L
is
the
total
num
ber
of
such
groups.
Furthermore,
a
considerableincremental
savingsin
addition
instructions
can
be
realized
by
a
recognition
that
with
increas
ing
distance
from
the
central
elementof
the
kernel,
i.e.,
at
its
symmetric
tail
ends,
the
number
of
consecutive
kernel
ele
ments
with
identical
numerical
valuesincreasesrapidly.
Rather
than
summing
all
projectiondata
values
to
be
multi
plied
by
a
givenkernel
element
every
time
the
index
j
is
incre
mented
(i.e.,
every
time
the
projection
is
indexed
past
the
kernel
by
one
element
position),
it
is
less
costly
tosubtractthe
oldest
projection
sample
from
the
running
sum
andadd
a
new
projection
sample
to
the
running
sum.
This
procedure
is
really
efficient
only
if
employed
for
the
elements
in
the
tails
of
the
kernel,
where
many
consecutive
elements
are
iden
tical.
This
procedure
is
expressed
by
the
relationship
g7.
=
g7

P(Rn)
+
P(Rna+)i
(1
0)
The
algorithm
described
abovewasused
to
generate
binary
filter
kernels
whose
performancewas
then
tested
on
simulated
and
real
Xray
and
ultrasound
projection
data.
The
fractional
band
kernels
and
their
binary
versions
were
testedfor
accuracy
of
reconstructed
density
values
for
Xray
projectiondata
as
well
as
for
ultrasound
timeofflight
data
[32],[33].
Recon
struction
of
the
same
projectiondata
sets,
both
with
binary
kernels
and
with
the
srcinalfullvalued
kernels,
followed
by
0.50
Z
0.25
0L
0
0.25
0.50
FREQUENCY
Fig.
1.
Modulus
of
Fourier
transform
of
convolution
kernels.
The
kernels
shown
are
RamLak
[14]
and
SheppLogan
[37]
and
their
binary
versions.
Since
kernels
are
symmetric
about
ordinate,
only
half
of
Fourier
domain
is
shown.
Reproduced
with
permission
from
[33].
comparison
of
the
resultingimages,
demonstrated
that
the
binary
kernel
reconstructions
were
both
stable
and
accurate.
However,
several
errors
in
the
reconstructeddensity
values
were
observed
near
the
skull
boundaries
of
very
high
contrastskull
phantoms
[37]
T
he
frequency
response
of
the
binary
kernels
is
given
in
Fig.
1,
demonstrating
the
close
approxima
tion
of
the
binary
kernels
to
thesrcinal
kernels.
If
in
some
applicationsthe
errors
in
reconstructed
image
density
values
caused
by
the
binary
kernel
approximations
are
intolerable,
the
followingrecently
developed
technique
for
improvement
of
the
results
can
be
employed.
Since
the
convo
lution
algorithm
is
both
associative
and
commutative,
the
con
volutionof
a
projection
with
an
arbitrary
fullvalued
kernel
can
be
subdivided
into
severalsuccessive
convolutions
of
the
projection
with
a
sequence
of
different
binary
kernels.
The
resulting
intermediate
convolved
projections
can
then
be
summed
to
yield
the
final
filtered
projection.
Forsuch
a
technique,the
binary
approximation
to
thesrcinal
kernel
is
obtained,
and
then
convolved
with
the
preweighted
projection
(e.g.,
using
binary
shift
and
add
instructions).
Next,
a
residual
convolutionkernel
is
obtained,
which
is
the
elementby
element
difference
between
thesrcinalkernel
and
the
binary
kernel.
A
binary
approximation
to
the
residual
kernel
is
then
generated,the
preweighted
projection
is
again
convolved,
this
time
with
theresidual
binary
kernel,
and
the
result
is
added
to
the
initial
convolved
projection
generated
by
the
first
con
volution
described.
If
desired,
this
successive
approximation
process
can
be
repeated.
In
effect,the
elements
of
the
full
valued
kernel
are
rewritten
as
sums
and
differences
of
binary
numbers.Note
that
the
central
element
of
the
kernel
is
always
a
nonbinary
value
and
hence
requires
the
execution
of
a
fixedprecision
multiplicationoperation
for
the
evaluation
of
its
contribution
at
each
stage
of
the
convolution.
This
technique
102