Journal of Medical and Biological Engineering, 30(5): 283-288
283
Development of a Multi-scale Noise Reduction Algorithm for
Cone-beam Micro-CT Using Fuzzy Logic-based
Anisotropic Filter
Zhiwei Tang
Guangshu Hu*
Department of Biomedical Engineering, Tsinghua University, Beijing 100084, China
Received 16 Sep 2009; Accepted 18 Nov 2009; doi: 10.5405/jmbe.30.5.03
Abstract
In micro-CT, smaller voxels require much higher dose to achieve adequate signal-to-noise ratio compared to
conventional CT, thus, necessitating a method to reduce statistical noise and then to improve dose efficiency. We
developed a method for multi-scale filtration of micro-CT data before reconstruction. The study aims are to increase the
signal-to-noise ratio (SNR) of image data while at the same time maintaining image quality, in particular spatial
sharpness of the images. In our approach, each projection is decomposed into successive sub-bands representing
different levels of resolution from coarse to fine. Within each sub-band, pixels are distinguished into high-contrast
regions (edges) and smooth regions using a fuzzy logic operation. Anisotropic filters are then used to adaptively
smooth the image data. The degree of smoothing is controlled by the factor obtained by fuzzy logic operation. Both the
fuzzy membership function and fuzzy rules are defined in such a way that the performance is optimized in terms of its
robustness to noise. The potential of our algorithm was evaluated using phantom and clinical data. It was shown that
SNR could be improved by a factor of 4 at the same radiation dose without loss of sharpness or introduction of new
image artifacts. Alternatively, on the basis of clinical performance, our method shows that the radiation dose can be
significantly reduced while at the same time maintaining the noise level and the image sharpness. Additionally, our
approach can be extended to other imaging methods as an efficient tool for edge-preserving noise reduction.
Keywords: Micro-CT, Noise reduction, Multi-scale, Fuzzy logic, Anisotropic filter
1. Introduction
Small-animal imaging has played a critical role in
phenotyping and drug discovery, and provided a basic
understanding of mechanism of disease. Among the newly
developed small-animal imaging modalities, micro-CT, with its
spatial resolution on the order of several microns, is believed to
have great potential for accurate 3D imaging of anatomy in a
minimally invasive way [1,2]. The signal-to-noise ratio (SNR)
is a special concern for micro-CT, especially if longitudinal
studies are required. The major source of measurement noise is
the x-ray photon flux. The smaller voxels of micro-CT (usually
less than 100 um) require much higher dose, since any voxel,
independent of its size, must interact with a certain number of
x-ray photons for adequate image quality. Therefore, in order to
compensate for the reduction of the voxel size, the photon flow
must be much larger in micro-CT to achieve a SNR comparable
to that of clinical CT. If the photon noise (measured by its
* Corresponding author: Guangshu Hu
Tel: +86-10-62784568; Fax: +86-10-62784568
E-mail: [email protected]
variance) is kept constant and the linear dimension of the
voxel is reduced by 2 (i.e. a reduction of voxel volume by 8)
the dose must be increased by up to 16 times [3,4]. Thus,
although desirable, a high enough SNR for micro-CT is
difficult to achieve, due to the associated increase in radiation
dose risk.
To improve noise characteristics or dose efficiency in
micro-CT, the aim of this study is to develop an
edge-preserving multi-scale filtering approach which treats the
projection data used for reconstruction. In contrast to
conventional filtering, our methods have the potential to reduce
random noise while trying to keep the structure information in
images [5]. By using fuzzy logic [6,7], we aim to define a
confidence measure with respect to the existence of edge by
adaptive judgment, where the judgment is adapted to the local
gray-scale statistics and to the context spatial of the pixel. The
anisotropic filters are subsequently used to filtering image on a
multi-scale structure. The filter is numerically robust (choice of
the parameters is not critical) and faster than robust statistics
estimation smoothing. Additionally, it can be straightly
integrated with fuzzy logic manner and multi-scale concept.
J. Med. Biol. Eng., Vol. 30. No. 5 2010
284
0
xa

BIG =  B ( x )  
xa
1
2. Methods
2.1 Fuzzy logic edge detector
In previous studies, binary edge detector was usually
introduced, which could distinguish only between edge and
no edge. This edge detector was used to control the spatial
filter by switching the filter on or off [8]. A disadvantage of
this hard switching was that in the case of higher noise level,
due to edge ambiguities, the filter introduced blurring and
artifacts.
Fuzzy logic is an extension of classical theory called
binary logic. It handles the concept of partial membership
(partial truth). In fuzzy logic theory, the membership degree
can be taken as a value intermediate between zero and one.
Fuzzy sets are generally represented by membership function,
i.e., a mapping from the underlying universe into the unit
interval [0, 1]. The membership function of the intersection of
two fuzzy sets A and B (A∩B) is modeled by triangular norm,
while the membership function of the union of two fuzzy sets A
or B (A∪B) is modeled by triangular co-norm (see Table 1 for
some examples) [7]. Additionally, we denote γA(ε) and γB(ε) as
the corresponding membership in the fuzzy sets A and B,
respectively, where ε is an element of universe U.
Table 1. Common fuzzy operations.
Fuzzy Intersection (A∩B) (γA(ε)∧γB(ε))
Algebraic product
γA(ε)γB(ε)
Fuzzy Union (A∪B) (γA(ε)∨γB(ε))
Algebraic sum
γA(ε)+γB(ε)- γA(ε)γB(ε)
In this study, we propose a fuzzy logic edge detector for
denoising. The determined degree of the edge confidence is
expressed as a real number between zero (no edge for sure) and
one (edge for sure). The uncertainties at lower processing
levels are combined in a specific way (using fuzzy rule) to
obtain a meaningful result at a higher processing level. This
provides more freedom for further processing tasks.
We assume that when an edge-like structure extends the
pixel (x, y) in a certain direction, it causes large derivative
values perpendicular to the direction. Thus, we use the
parameter θ(x, y) to measure whether an edge extends pixel
(x, y). θ(x, y) is calculated as follows:
θ(x, y) = σ - k∣Ip(x, y)- μ∣
if x  a
if a  x  b
(2)
if x  b
that is, the possibility where edge exists increases for higher γB
value. Specifically, between x = a and x = b, the membership
function increases linearly.
Subsequently, γB of the specific neighbors from the 2-D
window (3 × 3) are analyzed and combined in a meaningful
way (defined by the fuzzy logic rule) to determine the edge
confidence measure α(x, y) for each pixel (x, y). The degree of
α(x, y) is expressed as a real number measuring the existence of
edge between two extremes: zero (no edge for sure) and one
(edge for sure). Considering the fact that when an edge extends
a pixel(x, y), it will occupy a group of connected pixels, we
design the fuzzy rule R (see Table 2) as there is an edge if γB of
the central pixel is BIG and at least three neighboring pixels
(x+i, y+j) (i, j = -1, 0, +1) in the 2-D window are also
considered as BIG (see Figure 1 for some examples). The
intersection (∧) of all combinations of γB(x, y) and three
different neighboring BIG membership functions γB(x+i, y+j)
(i, j = -1, 0, +1) are computed, where the algebraic product is
used as the intersection operator. The results of these outcomes
are summed, using the∨ operator, which is the algebraic sum
(see Table 1) of all possible instances(see Table 2) to obtain the
edge confidence α(x, y). The activation degree of R indicates to
what degree α(x, y) can be considered as true. The edge
confidence α(x, y) will be used to control the strength of
smoothing in adaptive filtering as explained in next section.
Table 2. Fuzzy rule description (∧: And; ∨: Or;  : there exists).
Subfact
SF1
SF2
SF3
SF4
Rule
R
Fact
γB(x,y) = BIG
γB(x+i1,y+j1) = BIG
γB(x+i2,y+j2) = BIG
γB(x+i3,y+j3) = BIG
Antecedent
i1 , j2 ,i2 , j3 ,i4 , j4 ( SF1  SF2  SF3  SF4 )
Consequent
There exists
an edge at
[(i1 , j1 )  (i2 , j2 )  (i3 , j3 )]  (i1,2,3 , j1,2,3 )  (0, 0) (x,y).
(1)
2
where Ip(x,y) is the pixel value at (x,y), σ and μ are the local
variance and mean values in the neighborhood, and k is a
constant coefficient usually between 1.5~2. If θ is big, it means
that Ip(x, y) is close to its mean value of neighborhood and the
local variances are probably generated by the edge-like
structure. Contrariwise, the local variances may be probably
produced by noise. Actually, there would be gradual transition
from small θ value caused by noise to large θ caused by edge.
To differentiate between noise and edge, we define the
membership function BIG for the signal θ(x, y), i.e., the
membership γB[θ(x, y)], as follows:
Figure 1. Examples of fuzzy rule in a sliding window for fuzzy logic
edge detector: B, big; X, unsure.
2.2 Fuzzy logic bilateral filter
The anisotropic filter is a classic example of an
edge-preserving noise reduction method. Bilateral filter is a
special class of anisotropic filter, which obtains the pixel value
Noise Reduction for Cone-Beam Micro-CT
at location x with an average of nearby and similar pixel values
by a combined domain (space) and range (similarity) filter
[9,10]. A commonly used bilateral filter is the shift-invariant
Gaussian filtering, in which both the domain (space) filter and
the range (similarity) filter function are Gaussian functions
[11]. For an image with value I(p) at pixel p, a bilateral
Gaussian filtering operation is defined by:
Gg ,h
1
I ( p) 
Np

g
( p  q ) h ( I ( p )  I ( q )) I ( q )
(3)
qD
where:
Np 

g
( p  q ) h ( I ( p )  I ( q )),  m (t )  exp(
qD
t
)
2m 2
2
(5)
The actual value of the statement “▽D(x, y) is SMALL” is
FS[▽D , h].
Combining with domain filter, we defined the proposed
fuzzy logic bilateral filtering in terms of the calculated fuzzy
derivatives (fuzzy rule R1), for each of the eight neighbors, as
follows:
1
I out ( x , y ) 
1
 
g
( k , l ) wr ( x  k , y  l ) I p ( x  k , y  l )
k 1 l  1
1
(6)
1
 
k 1 l  1
g
( k , l ) wr ( x  k , y  l )
Table 3. Pixel positions involved in the calculation of derivatives.
Direction D
Neighboring derivative
Central derivative
 DD ( x, y )
 CD ( x, y )
NW
 D ( x  1, y  1)
 D ( x, y )
W
 D ( x  1, y )
 D ( x, y )
SW
 D ( x  1, y  1)
 D ( x, y )
S
 D ( x, y  1)
 D ( x, y )
SE
 D ( x  1, y  1)
 D ( x, y )
E
 D ( x  1, y )
 D ( x, y )
NE
 D ( x  1, y  1)
 D ( x, y )
N
 D ( x, y  1)
 D ( x, y )
(4)
D denotes a sliding window of the neighboring pixels centered
at the pixel p. g and h are filtering parameters measuring the
width of φg and φh, respectively.
The general idea behind our method is to integrate the
edges feature orientations into bilateral filter while filtering the
unstationary noise. Consider a 3 × 3 neighborhood of a pixel
Ip(x, y). The derivative ▽ D(x, y) is defined as the absolute
difference between the central pixel (x, y) and its neighbor in
the direction D; D∈{NW,W,SW,S,SE,E,NE,N}. For example,
▽N(x, y) = ∣Ip(x, y-1)-Ip(x, y)∣.
It is assumed that when an edge structure extends in a
certain direction D, it causes large derivative values
perpendicular to the direction D at the current pixel position (x, y)
and its neighbor pixels as well. For example, if a edge extends in
W to E direction, either ▽N(x, y) together with ▽N(x, y-1) or
▽S(x, y) with ▽S(x, y+1) are expected to be relatively large. In
such a manner, by taking into account those derivative values for
each direction and combining them in a fuzzy logic way, we
determine the averaging weight for every neighboring pixel. We
use the following logic rule R1 to determine the fuzzy derivative
value for a certain direction D at (x, y) position. The rule R1 is
C
described as, if the derivative of central pixel  D ( x, y ) is
D
SMALL and the derivative  D ( x, y ) of neighboring pixel in
direction D is also SMALL, we conclude that there is SMALL
derivative value SDD(x, y) in direction D, i.e. a small edge
magnitude along the direction perpendicular to D. In Table 3, we
give an enumeration of the pixel positions that are involved in the
calculation of central and neighboring derivatives. The
intersection (and) operator is defined using algebraic product (see
Table 1). The fuzzy membership function SMALL, decreases for
larger▽D, is defined according to the Gaussian function:
FS[▽D , h] = φh(▽D)
285
The weighting coefficient wr for neighboring pixel in D’ th
direction, which will replace the range filter coefficient in
regular Gaussian bilateral filtering, is defined as the output of
the fuzzy controller fixed by rule R1, i.e., the weight for the
D’ th direction is equal to SDD(x, y). The smaller the derivative
 CD ( x, y ) and  DD ( x, y ) , the larger will be weighting
coefficients wr, and hence the corresponding neighboring pixel
values will be taken more into account for the filtering.
Considering the edge confidence α(x, y) calculated in the
last section, we can modify the filtering formula in order to
control the amount of filtering as follows:
I out ( x , y )  I p ( x , y )  (1   ( x , y ))
1
1
 
( k 1 l 1
g
( k , l ) wr ( x  k , y  l ) I p ( x  k , y  l )
1

1
  g ( k , l ) wr ( x  k , y  l )
 I p ( x , y )) (7)
k 1 l 1
Currently, the smoothing parameter g and h are determined
empirically.
2.3 Multi-scale analyzing and filtering
Multi-scale methods decompose the image into multiple
frequency bands and perform processing separately for each
band. Edges and noise mainly belong to the high-frequency
portion, and contrast information is conveyed primarily in the
low-frequency. The proposed fuzzy logic-based analysis and
filtering is integrated into a multi-scale structure to process
micro-CT images. As shown in Figure 2, each image is
decomposed into a succession of sub-bands using Laplacian
Pyramid [12]. Within each sub-band, pixels are analyzed and
segmented into highly textured region (edges) and weakly
textured region using fuzzy logic edge detector. Highly
anisotropic filtering is applied to the highly textured region
using fuzzy logic-based bilateral kernels that are elongated
along the edges and that will not blur in the direction
perpendicular to these edges. In the weakly textured region,
slight anisotropic filtering is applied to reduce high-frequency
random noise. The result of the analysis phase is used to
control the amount of the filtering. As the noise is mainly
contained in the high-frequency levels, the low-frequency
levels will be left unfiltered to maintain the coarser contrast
information. Those operations are applied to each of the
individual sub-bands before recomposing the pyramid, as
indicated in the figure.
286
J. Med. Biol. Eng., Vol. 30. No. 5 2010
Low-pass and
down-sample
Fuzzy logic
edge detector
Anisotropic
filtering
Sum and
up-sample
Low-pass and
down-sample
Fuzzy logic
edge detector
Anisotropic
filtering
Sum and
up-sample
Low-pass and
down-sample
Up-sample
Figure 2. Block diagram of multi-scale analyzing and filtering.
2.4 Experiments and data processing
The algorithms were evaluated using measured projection
data acquired with our micro-CT system. The system is a
bench-based cone-beam computed tomography system with
cone angle of 9.5º. The micro-CT system consists of a
Hamamatsu flat panel detector, an Oxford micro-focal x-ray
tube and a rotational bench. The x-ray tube can be operated
with 0~50 kV, 0~1 mA anodic current, and 50-μm focal spot.
The detector is a 50-mm Ti:CsI-CMOS device with 50-μm
pixels. The radiological constants are 40 kVp, 0.2 mAs per
frame. A 0.1 mm aluminum filter is placed between the source
and the object. Projection images are reconstructed into a 3D
volume (400 × 400 × 400) by the short-scan Feldkamp
algorithms from images of 190 view angles. The pixel pitch of
reconstructed volume is 0.1 mm.
Projection image was first decomposed into a four-level
Laplacian Pyramid by recursive low-pass filtering and
subtraction. Since the random noise is mainly contained in
high-frequency band, the lowest frequency level is left
unfiltered. Each high-frequency level was first processed with
fuzzy logic edge detector to obtain two-dimensional edge
confidence data α(x, y). The parameter k in Eq. (1) was set as
1.8. The parameters a and b in Eq. (2) determine the slope of
membership function. Note that a and b are adapted to depend
on the global noise level σav, which is represented by extraction
of average local variance of all pixels. The choices of a = 0 and
b = 2.5σav were experimentally found to be best in terms of
optimal division between noise and edge. Then the band image
is smoothed with anisotropic filter. The parameter g of domain
filter is chosen as 1 pixel to avoid unnecessary blurring of
coarser structures. The parameter h is set according to the
global noise level σav as h = 2σav. The resulting image is then
recomposed from the filtered details and low-pass image.
Figure 4 illustrates the reconstructed images from original data,
linear Gaussian (5 × 5) filtered data and our multi-scale
anisotropic filtered data. The result shows that our filter has
effectively reduced noise and maintained structure information,
in contrast to the loss of edge caused by the conventional
linear Gaussian filter.
For more quantitative comparison, we calculate the
signal-noise-ratio (SNR) of a uniform area (marked with a
square in Figure 4). The SNR is defined as the local mean
divided by the local standard deviation (Stdv) of image
densities in the selected area. Table 4 lists the Means, Stdvs,
and SNRs of the uniform area. The SNR of the image
reconstructed from multi-scale anistropic filter is almost the
same as that of the images reconstructed from linear Gaussian
filter, without removing intrinsic structure.
Figure 3. PMMA-air phantom.
(a) Original image
(b) Gaussian filtering
3. Results
To investigate the efficiency of our approach, a
thorax-like phantom was scanned and tested. The phantom is a
cylinder made of poly(methyl methacrylate) (PMMA) with
two perforative air holes (see Figure 3). The high contrast
between the air and PMMA region is well representative of the
contrast between the lungs and soft tissue in the thorax. The
image quality and the influence of our filtering algorithms can
be judged best by regarding the image themselves and the
difference images of the filtered data minus the original image.
(c) Multi-scale anisotropic filtering
Figure 4. Reconstructed slice of thorax-like phantom with filtered
projection images (left), difference image between original
and filtered one (right).
Noise Reduction for Cone-Beam Micro-CT
287
Figure 5. Schematic plot of the method of on-line MTF calculation: (a) the cross section image of the air hole is used as a sharp edge; (b) the
pixel values across the edge are used as to derive the LSF; (c) the Fourier transform of the LSF is then calculated to determine the
MTF; (d) MTF.
In order to quantify the spatial resolution of the processed
image, we proposed an on-line modular transfer function (MTF)
measurement method. The cross-sectional view of the air hole
in the reconstructed image provides a sharp edge (Figure 5).
The line spread function (LSF) can be obtained by derivative of
the pixel values crossing the edge. Then, the MTF is calculated
through the Fourier transform of the LSF (Figure 5). Using this
method, the spatial resolution (quantified by MTF) of
reconstructed images by different filtering methods was
calculated. The MTF for various noise reduction algorithms are
shown in Figure 6. In order to measure the improvement of
dose efficiency, we reconstructed high-dose projections whose
dose is 5 times higher than regular projections. The results
indicate that the linear Gaussian filter produced a much lower
spatial resolution, while our method provided almost the same
spatial resolution as high-dose images.
(a) Original image
High dose
Gaussian filter
(b) Gaussian filtering
Multi-scale
anisotropic filter
(c) Multi-scale anisotropic filtering
Figure 6. MTF for various noise reduction algorithms.
We also used our filtering algorithms with low-dose true
rat data to demonstrate its applicability for clinic. The results
are shown in Figure 7. Noise has been reduced about 75%
without losing edge after fuzzy multi-scale anisotropic filtering.
It can be clearly see that the overall image quality is greatly
improved with our denoising method.
4. Discussion and conclusion
The aim of this study was to improve the
reconstructed-image quality of micro-CT, which is degraded
Figure 7. Reconstructed slice of true rat with filtered projection images
(left), different image subtracting original image (right).
by the random noise caused by insufficient radiation dose. As
described above, our filtering method effectively reduces noise
without losing structure information. Averagely, noise can be
typically reduced to 25% of its original level.
Iterative reconstruction methods with Bayesian priors can
also reduce noise without blurring noise, but need a large
number of iterations to ensure convergence [13]. The fuzzy
logic operation is entirely based on a series of multiplicative
and additive operations which are computationally inexpensive
and will therefore not increase the calculation time significantly.
Our method contains several parameters to be tuned
J. Med. Biol. Eng., Vol. 30. No. 5 2010
288
experimentally. However, for both high and low noise level, the
fuzzy logic-based multi-scale bilateral filtering gives excellent
results with the same parameters set. Therefore, the selection of
parameters is not critical and the method is robust in use. For
the speed and robustness, our filtering method combined with
FBP reconstruction is suitable for the routine clinical use.
We have demonstrated the fuzzy logic based multi-scale
bilateral filtering is a sufficient way to reduce noise without
losing structures information for micro-CT images. This new
method has great potential in future micro-CT applications as
an available tool to improve image quality and increase dose
efficiency by means of data preprocessing. This will make
micro-CT studies more feasible.
[4]
[5]
[6]
[7]
[8]
[9]
[10]
References
[1]
[2]
[3]
P. Ruegsegger, B. Koller and R. Muller, “A microtomography
system for the nondestructive evaluation of bone architecture,”
Calcif. Tissue Int., 58: 24-29, 1996.
M. J. Paulus, S. S. Gleason, S. J. Kennel, P. R. Hunsicker and D.
K. Johnson, “High resolution X-ray computed tomography: an
emerging tool for small animal cancer research,” Neoplasia, 2:
62-70, 2000.
N. L. Ford, M. M. Thornton and D. W. Holdsworth,
“Fundamental image quality limits for microcomputed
tomography in small animals,” Med. Phys., 30: 2869-2877,
2003.
[11]
[12]
[13]
N. Beckmann, R. Kneuer, H. U. Gremlich, H.
Karmouty-Quintana and M. Muller, “In vivo mouse imaging and
spectroscopy in drug discovery,” NMR Biomed., 20: 154-185,
2007.
C. K. Chu, I. K. Glad, F. Godtliebsen and J. S. Marron,
“Edge-preserving smoothers for image processing,” J. Amer.
Statistical Assoc., 93: 526-541, 1998.
D. van de Ville, M. Nachtegael, D. van der Weken and E. Kerre,
“Noise reduction by fuzzy image filtering,” IEEE Trans. Fuzzy
Syst., 11: 429-436, 2003.
E. E. Kerre, Fuzzy Sets and Approximate Reasoning, Xian,
China: Xian Jiaotong University Press, 1998.
V. Zlokolica and W. Philips, “Motion and detail adaptive
denoising of video,” Proc. SPIE, 5298: 403-412, 2004.
P. Perona and J. Malik, “Scale-space and edge detection using
anisotropic diffusion,” IEEE Trans. Pattern Anal. Mach. Intell.,
12: 629-639, 1990.
S. Weeratunga and C. Kamath, “A comparison of PDE-based
non-linear anisotropic diffusion techniques for image
denoising,” Proc. SPIE, 5014: 201-212, 2003.
F. Godtliebsen, E. Spjstvoll and J. S. Marron, “A nonlinear
Gaussian filter applied to images with discontinuities,” J.
Nonparametr. Statist., 8: 21-43, 1997.
P. J. Burt and E. H. Adelson, “The Laplacian pyramid as a
compact image code,” IEEE Trans. Commun., 31: 532-540,
1983.
J. A. Fesseler, E. P. Ficaro, N. H. Clinthorne and K. Lange,
“Grouped coordinate ascent algorithms for penalized likelihood
transmission image reconstruction,” IEEE Trans. Med. Imaging,
16: 166-175, 1997.