It has long been established that the appropriate way of reslicing volume MR images is to use the method of sinc interpolation [2, 4]. We have recently needed to implement this method ourselves and have found, like other authors before us, that large convolution kernels are needed in order to produce accurate reslice data, suitable for subtraction. This requirement has led many groups to investigate the use of specialised hardware and software in order to perform data analysis within sensible timescales. However, we have found that the major component of the error introduced from interpolation with small kernels, is actually due to a first order normalisation problem introduced by truncation. In this paper we demonstrate the characteristics of this problem on real data and show how it can be eliminated, so that accurate reslice data can be obtained with small kernels. Unlike other recent suggestions for correcting such effects [3], the required changes in computation are simple and significantly reduce the processing requirement for a given interpolation accuracy. Renormalised Sinc Interpolation. There are several techniques that one can adopt to solve the problem of image interpolation. One is to assume a particular prior functional model for a local region of the image data, estimate the function parameters from a maximum likelihood metric and then recompute intermediate sites from the functiona

A. Jackson

B. Vokurka

D. Moriarty

N. A. Thacker

CiteSeerX

Tina Memo No. 2000-011
Short Version published in: Proc. MIUA, pp 7-10, London. 10th-11th July, 2000.
Modal Division and its Application to
Medical Image Analysis.
N.A.Thacker and A.J.Reader.
Last updated
18 / 12 / 2000
Imaging Science and Biomedical Engineering Division,
Medical School, University of Manchester,
Stopford Building, Oxford Road,
Manchester, M13 9PT.Modal Division and its Application to Medical Image
Analysis.
N.A.Thacker1, A.J.Reader2.
1: Imaging Sciences and Biomedical Engineering, University of Manchester
2: Department of Instrumentation and Analytical Science, UMIST
email: nat@sv1.smb.man.ac.uk
Abstract
Statistical instability is the problem associated with the calculation of poorly dened variables from data
corrupted by noise. The problem of statistical instability of numerical algorithms is well known in the
image processing literature. Problems with instability of a small number of calculations is never worse
than when attempting Fourier deconvolution. Each parameter in the Fourier domain aects every data
point in the spatial domain so that unstable calculation of just one parameter distorts the solution for
the entire space. The standard solution to this problem is often referred to as "optimal" or "Wiener"
ltering. However, Weiner ltering is not found to yield good solutions in images and other authors have
suggested alternative (empirical) techniques for solving the problem in practical situations. This paper
explains the un-suitability of Weiner ltering for image deconvolution and denes a new technique based
on an intuitive use of maximum likelihood, which we have called \modal arithmetic", and can be shown
to generate the empirical methods used in positron emission tomography (PET) as a special case. The
method is demonstrated on the deconvolution of 3D PET data.
Introduction
Image reconstruction directly from back-projected images is still a method of choice for novel 3D PET
systems (eg: PETRRA [1]), as spatial sampling accuracy of the system is retained whilst giving signicant
data size reduction. The reconstruction is normally carried out by 3D deconvolution of the system point
response function. However, for noisy images this can yield unstable results so generally some form of
constrained deconvolution is used.
Weiner ltering is based on a maximum likelihood derivation for the optimal linear lter which needs
to be applied to data in order to give the minimum error in the spatial domain after deconvolution
[2]. The technique involves estimating the noise in the power spectrum of the data (numerator) and
setting any terms in the Fourier domain within some number of deviations from zero to that value. The
problem with this approach is that it is still possible for very small values in the Fourier spectrum of the
deconvolution kernel (denominator) to amplify residual noise (beyond the selected noise threshold) to very
large values. The probability that this can happen increases as the scale of the Fourier space increases.
Unfortunately, the number of Fourier terms in the FFT of an image is generally very large. Ultimately,
the number of deviations from zero which must be used in order to ensure stable division is so great that
the resulting deconvolution is compromised. In practical applications a technique which has been applied
widely involves adding an oset term to the denominator. This limits the smallest value of the division
and ensures a stable result but it is not clear how this technique is statistically justied compared to
Wiener ltering. However, there is a more direct way of assessing the maximum likelihood value of a
calculation which does not necessarily require the assumption of a linear lter solution. In the remainder
of this paper we will examine the direct application of maximum likelihood to noisy denominators. We
will show how the method most commonly used for the deconvolution of 3D PET data [3] follows from
this alternative derivation with particular assumptions for the behavior of the data. The implications for
design of image processing algorithms is far broader than this however, as the general approach can be
applied to any non-linear calculation in order to regain computational stability.
2Methods
For division (eg: y = 1=x), small errors on the data produce instabilities in computations involving large
quantities of data. Error propagation shows that a small change in the input quantity x will give an
error on the corresponding output of
y =
x
x2
which is clearly unstable for values of x which are comparable to its error. This problem can be understood
better by considering the distribution of computed values from the range of those available for input. We
start by assuming a Gaussian distribution for the denominator.
Pxx = A exp( (x   x0)2=22)
Where x0 is the central value of x with a standard deviation of . If we take a small area of data from
the probability distribution for x (ie: Pxx), we can associate this with an equal number of solutions in
the output space y (ie: Pyy) (gure 1 (a) (i) and (ii)) giving:
Py = A x
2exp( (x   x0)
2=2
2)
￿
￿
￿
￿
￿
￿
￿
￿
￿
￿
￿
￿
￿
￿
￿
￿
￿
￿
￿
￿
￿
￿
￿
￿
￿
￿
￿
￿
￿
￿
￿
￿
￿
￿
￿
￿
￿
￿
Py
0
area = P   D y y
y= 1/x
x
P
0
y
X
X max +
max -
x
Px
0
area = P   Dx x
(i) (ii)
(iii)
(a) Probability distributions for x and y.
P
0 mod(z)
z
(ii)
Re(z)
Im(z)
(i)
projection
axis
(b) Probability distributions for z and signed
mod(z)
Figure 1: Probability Distributions for a real and complex denominator.
This expected probability distribution for y as a function of x (gure 1(a)(iii)) can be dierentiated to
nd its maxima. Setting this to zero we can determine the modal values of this distribution:
x2   x0x   22 = 0 with xmax =
x0 
p
x2
0 + 82
2
which correspond to the positive and negative peaks due to the distribution of x spanning zero (gure
1(a)(i)). If we were to ask which value of y would be most likely to result from the division then the answer
would be 1=xmax selected with the same sign as the input value x0. Taking this value as a replacement
for the denominator provides a maximum likelihood technique of stabilising the process of division using
knowledge of measurement accuracy and could best be described as modal division. Modal division
can be used with impunity for calculations involving large quantities of noisy data without instability
problems for values around zero, with the minimum denominator limited to a value of
p
2.
A simple symmetry argument can be used to show that the technique can be easily extended to complex
numbers z, with equal independent Gaussian errors on both real and imaginary components (gure
1(b)(i)), by applying the correction to the magnitude of the complex number (gure 1(b)(ii)).
mod(zmax) =
mod(z0) +
p
mod(z0)2 + 82
2
and zmax = z0 mod(zmax)=mod(z0)
The general method of modal arithmetic for a measured value with distribution D(x) and a non-linear
function f(x) would be to nd the solution xmax of
@[
D(x)
@f(x)=@x
]=@x = 0
3with the modal solution of f(xmax). Modal arithmetic is unconditionally stable, as peaks in probability
distributions cannot occur at innity. It also has much similarity with some approaches in statistics
which advocate the use of the mode rather than the mean as the most robust indicator of a distributed
variable.
It is not immediately obvious that the method of modal division can be applied to deconvolution. The
standard approach proceeds as follows:
I = U:R implies U = I=R = IR
=(R
R)
where I and U are the complex Fourier coecients from the image data and signal source respectively, and
R is a coecient from the assumed convolution kernel. From Parsevals theorem and the orthogonality of
the Fourier domain, we know that uniform independent Gaussian errors in the spatial domain correspond
to uniform independent Gaussian errors in the spatial frequency domain. The errors on the data I are
generally assumed to come from the measurement process, either in the spatial or spatial frequency
domains. However, the observable error I can just as easily be generated from noise in each component
of the convolution kernel R , ie:

2
I = mod(U)
2
2
R
This change in the assumed image formation process makes no dierence to the observed data but makes
modal division appropriate 1. Using this relationship we can compute the error on a coecient in the
kernel R needed to generate the expected error in the measured image 2. We can then apply the method
of modal division to Fourier deconvolution.
mod(Rmax) =
mod(R) +
p
mod(R)2 + 82
I=mod(U)2
2
This method can be implemented as an iterative calculation
mod(Rt) =
mod(R0) +
p
mod(R0)2 + 82
I=mod(Ut 1)2
2
Ut = IR
t=R
tRt
with U0 = I and R0 = R. Since convolution generally involves attenuation of the value of coecients
(ie: I < U), this is expected to be computationally stable as early iterations overestimate the error on
R. Further, as the calculation involves only a slight modication to the computed values it is expected
to converge rapidly.
Results
The above theory is directly applicable to calculations in medical image processing involving a division
of voxel values. One example of this is in the calculation of perfusion ow from relative cerebral blood
volume (RCBV) and mean transit time (MTT) [4, 5]. However, the most convincing demonstration of
the method is with Fourier deconvolution. One of the standard applications of this is in the deconvolution
of 3D PET images after back projection. The expected deconvolution kernel is a symmetrical function
with 1=r2 dependency. PET data has Poisson noise characteristics and this leads to correlations in the
noise in the spatial frequency domain. However, we need only know the individual distributions for each
component in order to apply modal division. Treated separately, the co-ecients have Gaussian noise
with equal variance as modelled above.
The method suggested by Chu and Tam in [3] is a constrained deconvolution which has been used by
other groups (eg: [7]). It is derived from an assumption of smoothness in the reconstructed image and of
the form
U
0(p) =
I(p)
R(p) + p4=R(p)
1Clearly there is an innite set of possible models with dierent levels of contribution to the statistical error from both
sources but this model is the one which will result in the greatest stability.
2Interestingly, the resulting kernel error
2
R = 2
I=mod(U)2
is also indicative of the accuracy to which this coecient must be determined in order to attempt deconvolution, with the
implication that contrary to conventional teaching we do not need highly accurate estimates of all kernel co-ecients in
order to get stable results.
4where U0(p) is an estimate of the constrained solution at spatial frequency p,  is a constant proportional
to the noise level and R(p) is real. The method is not dened for asymmetric kernel functions (ie:
complex R(p)). Assuming that R(p) is small (which would therefore be potentially a cause of instability)
the correction terms in the two approaches can be directly compared, ie:
I
p
2R(p)2=mod(I(p))2  p
4=R(p) implies mod(I(p)) / R(p)
2p
 4
The scaling with noise is directly equivalent. The Chu and Tam result can be interpreted as equivalent
to modal division only on the assumption of a typical relationship between the power spectra of the data
and kernel. This result is in keeping with the assumptions of the original paper. The practical dierence
between the two approaches is that those occasions where the algorithm fails can often be recovered by
adjusting the value of . For modal division, this parameter can be replaced by a value computed using
error propagation.
Typical results from deconvolution of back-projected 3D PET images for the standard technique and
modal division are given in Figure 2 below. The image is an 18-FDG (uorodeoxy-glucose) scan of a
mouse (Fluorine 18 positron emitter). The data were Acquired on the HIDAC (II) [6] scanner at the
Paterson Institute for Cancer Research. The slice is 128x128 voxels, each of 1 mm side length. (from
a 128x128x128 matrix). The scanner has two main rotating planar detectors: each with two HIDAC
detection modules, of size 320mm x 320 mm with a detector separation of 714 mm. There is clearly very
little dierence between the nal results which implies that the required relationship between the power
spectra are indeed quite close. However, one would have to worry about the stability of the Chu and
Tam method for other deconvolution kernels.
(a) Back-projected Image (b) Modal Division
(c) Constrained Deconvolution
Figure 2: Performance of modal division and constrained deconvolution on PET data.
Conclusions
The problem of statistical stability occurs frequently in medical image analysis but is probably never
more obvious than in the deconvolution of image data, where denominators comparable to the noise level
are frequent yet every coecient must be computed stably in order to give an accurate result. Common
5methods for dealing with such data based on maximum likelihood have key failings, while more ad-hoc
techniques have often been more successful. This demonstrates that statistical methods have not yet
been correctly applied to the problem. Ultimately, appropriate statistical methods should give the best
performance. We have dened a new approach to the design of unconditionally stable algorithms based
upon maximum likelihood which we call \modal arithmetic". We have been able to show how the most
common (practical) technique for deconvolution in the 3D PET images can related to the method derived
using this technique. The results obtained, do not demonstrated the new technique to be any better than
the empirical method on this data. A more systematic evaluation is required if we wish to establish this.
In fact we have recently been able to show that the convergence point of our iterative algorithm is the
same as using Wiener Filtering. However, the results do illustrate the validity of this new methodology
for designing stable computational processes. The basic approach to the implementation of non-linear
calculations in noisy data sets has scope for general applicability. The method of modal image division
is available within the TINA software which is available as open source from our web site [8].
References
1. D.M.Duxberry et. al. Prelimiary Results from the new Large PETRRA Positron Camera, IEEE Trans.
Nucl. Sci. 46 , 1050-1054, 1999.
2. W.H.Press B.P.Flannery S.A.Teukolsky W.T.Vetterling, Numerical Recipes in C, Cambridge University
Press 1988.
3. G.Chu and K.Tam, Three Dimensional Imaging in the Positron Camera Using Fourier Techniques. Phys.
Med. Biol, 22, 2, 254-265, 1997.
4. L.Ostergaard, A.G. Sorensen, K.K. Kwong, R.M.Weissko, C.Gyldensted and B.R.Rosen. High Res-
olution Measurement of Cerebral Blood Flow Using Intravascular Tracer Bolus Passages, Parts I and
II,Experimental Comparison and Preliminary Results. MRM, 36, 715-736, 1996.
5. N.A.Thacker, X.P.Zhu, M.Nazapour, C.Moonen and A.Jackson, A New Approach to the Estination of MTT
in Bolus Passage Pefusion Techniques. This conference.
6. J.C.Hand "HIDAC camera design for positron emission tomography" PhD thesis, Univ. of Manchester UK,
1997.
7. S.Webb, R.J.Ott, J.E.Bateman, A.C.Flesher, M.A.Flower, M.O.Leach, P.Marsden, O.Khan and V.R.McCready.
Tumour localisation in oncology using positron emitting radiopharmaceuticals and a multiwire proportional
chamber . Nucl. Inst. Meth. 221 233-24, 1984.
8. URL: www.niac.man.ac.uk
6

Last updated

http://www.tina-vision.net/docs/memos/2000-011.pdf

Last updated

Abstract

Similar works

Full text

Available Versions

CiteSeerX