Projection space denoising with bilateralfiltering and CT noise modeling for dose reduction in CT
Armando Manduca a͒
Department of Physiology and Biomedical Engineering,Mayo Clinic College of Medicine,
Rochester,Minnesota55905and Department of Radiology,Mayo Clinic College of Medicine,
Rochester,Minnesota55905
Lifeng Yu
Department of Radiology,Mayo Clinic College of Medicine,Rochester,Minnesota55905
Joshua D.Trzasko and Natalia Khaylova
Department of Physiology and Biomedical Engineering,Mayo Clinic College of Medicine,
Rochester,Minnesota55905
James M.Kofler,Cynthia M.McCollough,and Joel G.Fletcher
Department of Radiology,Mayo Clinic College of Medicine,Rochester,Minnesota55905
͑Received15May2009;revised8August2009;accepted for publication29August2009;
published2October2009͒
Purpose:To investigate a novel locally adaptive projection space denoising algorithm for low-dose
CT data.
Methods:The denoising algorithm is based on bilateralfiltering,which smooths values using a
weighted average in a local neighborhood,with weights determined according to both spatial
proximity and intensity similarity between the center pixel and the neighboring pixels.Thisfiltering
is locally adaptive and can preserve important edge information in the sinogram,thus maintaining
high spatial resolution.A CT noise model that takes into account the bowtiefilter and patient-
specific automatic exposure control effects is also incorporated into the denoising process.The
authors evaluated the noise-resolution properties of bilateralfiltering incorporating such a CT noise
model in phantom studies and preliminary patient studies with contrast-enhanced abdominal CT
exams.
Results:On a thin wire phantom,the noise-resolution properties were significantly improved with
the denoising algorithm compared to commercial reconstruction kernels.The noise-resolution prop-
erties on low-dose͑40mA s͒data after denoising approximated those of conventional reconstruc-
tions at twice the dose level.A separate contrast plate phantom showed improved depiction of
low-contrast plates with the denoising algorithm over conventional reconstructions when noise
levels were matched.Similar improvement in noise-resolution properties was found on CT colonog-
raphy data and onfive abdominal low-energy͑80kV͒CT exams.In each abdominal case,a
board-certified subspecialized radiologist rated the denoised80kV images markedly superior in
image quality compared to the commercially available reconstructions,and denoising improved the
image quality to the point where the80kV images alone were considered to be of diagnostic
quality.
Conclusions:The results demonstrate that bilateralfiltering incorporating a CT noise model can
achieve a significantly better noise-resolution trade-off than a series of commercial reconstruction
kernels.This improvement in noise-resolution properties can be used for improving image quality
in CT and can be translated into substantial dose reduction.©2009American Association of
Physicists in Medicine.͓DOI:10.1118/1.3232004͔
Key words:CT,radiation dose,noise reduction,bilateralfiltering,denoising
I.INTRODUCTION
Radiation exposure and the associated risk of cancer for pa-tients receiving CT examinations have b
een an increasing concern in recent years.1,2It is critically important to reduce the radiation dose level in CT examinations.However,dose reduction generally leads to an increased level of noise in the measured projection data and the subsequent reconstructed images,thus degrading the diagnostic value of the CT exami-nation.
Many techniques have been developed for controlling noise in CT,operating on either the raw projection measure-ment,the log-transformed sinogram,during image recon-struction,or on images after reconstruction.3–10In conven-tional shift-invariantfiltration applied during image reconstruction,the suppression of the high-frequency com-ponent in the sinogram is performed with a simple assump-tion that all the measurements are equally reliable,which may result in severe degradation of spatial resolution.3–8 More sophisticated methods have been developed to adap-
49114911 Med.Phys.36…11…,November20090094-2405/2009/36…11…/4911/9/$25.00©2009Am.Assoc.Phys.Med.
tively smooth the data by taking into account the local sta-tistics in the measurements.Hsieh used an adaptive trimmed meanfilter to smooth the data with thefilter strength chang-ing with the detected signal.4The method of Kachelriess et al.was similar except that it was generalized to include both the
2D detector array and the angular dimension and used a direct threshold on projection data as opposed to a more accurate noise model.5Some of these methods are currently implemented on clinical scanners mainly to suppress the streaking artifacts caused by x-ray photon starvation.Many other approaches have also been proposed to incorporate more explicit statistical models and to iteratively restore the log-transformed data by optimizing a penalized weighted least-squares or likelihood objective function,in some cases by relating the sinogram value to its variance,thus empiri-cally characterizing the noise of the sinogram value after processing͑beam-hardening correction,calibration,etc.͒.6,7 Other approaches model CT noise in a more complete way, including compound Poisson,off-focal,cross-talk,and other effects.8,9Iterative reconstruction methods can achieve sig-nificant denoising but at the expense of very long computa-tion times.10Techniques based entirely on image space have also been described,taking advantage of the image structure to smooth noise while preserving edges but suffering from the complicated properties of noise in image space in CT.11–13
In this work,we investigated a locally adaptive method for noise control in CT.This method is based on bilateral filtering,14which smooths the sinogram by using a weighted average in a local neighborhood,with the weights deter-mined according to both the spatial proximity and intensity simil
arity between the center pixel and the neighboring pix-els.Thisfiltering is locally adaptive and can preserve impor-tant edge information in the sinogram,thus maintaining high spatial resolution.It is closely related to anisotropic diffusion15but is much faster.Furthermore,since it origi-nated from the statistical framework of maximum a poste-riori͑MAP͒estimation,16a CT noise model can be incorpo-rated,which is critical for effective noise reduction in low-dose CT.Note that Demirkaya17proposed anisotropic diffusionfiltering in projection space,but with a crude noise model and only showing results on a synthetically corrupted Shepp–Logan phantom.
While some denoising or sinogram restoration methods have been developed to smooth the projection data by taking explicit statistical models into account,6–9the proposed ap-proach is noniterative and is easier to implement in practice. Compared to adaptivefilters developed to suppress the streaking artifacts caused by the photon starvation along some directions,4,5the proposed approach aims to perform noise reduction in the full projection dataset while preserving structural detail.
This paper is organized as follows.In Sec.II,we intro-duce bilateralfiltering and its application to CT data denois-ing,describe the effects of the bowtiefilter and automatic exposure control and their inclusion in the CT noise model,and present the evaluation methodology.In Sec.III,the evaluation res
ults are presented and discussed.Finally,we present conclusions in Sec.IV.
II.METHODS
II.A.Bilateralfiltering
Bilateralfiltering was originally proposed by Tomasi and Manduchi as a noniterative and locally adaptive method for removing additive noise from images while preserving edge information.14In addition to its intuitive appeal,the bilateral filter also has strong origins in MAP estimation.16 Suppose Q is a stationary Gaussian process with a mean of Q
Q=Q+x,͑1͒where x is a white Gaussian noise added to the original im-age Q.We use a bold letter and the corresponding normal letter to denote a stochastic process and its mean,respec-tively.To restore the original image Q from the noise con-taminated image Q,we consider minimizing the energy func-tional
E͑Q͒=͚iI͚j⍀
i
P1͑i,j͒P2͑i,j͒,͑2͒
where P1calculates a weight according to the spatial dis-tance between the center pixel and a neighboring pixel and the second penalty factor P2calculates a weight according to the pixel value difference between the center pixel and a neighboring pixel.Additionally,j is the index of a neighbor-hood pixel inside a region of⍀i centered on pixel i,and the outer sum is taken over all pixels i in image I.Minimizing this functional encourages local groups of pixels to be of similar intensity,resulting in an image with overall piecewise homogeneity.Moreover,in the restricted case where only immediate neighbors are investigated and the spatial penalty function P1is taken as the identity operator,Eq.͑1͒is the same functional that is minimized by the classical aniso-tropic diffusion process for image denoising.15
While there are many possibilities for choosing the pen-alty functions P1and P2,the most common assignments are the Gaussian penalties given by
P1͑i,j͒=expͩ−͑i−j͒22d2ͪ,͑3͒P2͑i,j͒=1−expͩ−͑Q i−Q j͒222ͪ,͑4͒
where parameters d andcan be used for controlling the spatial and intensity range of the weighting.Following Elad,16the bilateralfilter arises as thefirst iteration of a Jacobi process minimizing Eq.͑2͒and is formally defined by the weighted averaging of a given image pixel Q i,namely,
Medical Physics,Vol.36,No.11,November2009
Qˆi=͚j⍀
i
W1͑i,j͒W2͑i,j͒Q j
͚j⍀
i
W1͑i,j͒W2͑i,j͒
,͑5͒
where W1and W2are the differentials of P1and P2and,after certain factors cancel,
W1͑i,j͒=P1͑i,j͒=expͩ−͑i−j͒22d2ͪ,
W2͑i,j͒=expͩ−͑Q i−Q j͒222ͪ,͑6͒
which are again Gaussian.Under the specification in Eqs.͑3͒and͑4͒with Gaussian weighting functions,the bilateralfilter is identical to the Nadaraya–Watson estimator commonly employed in nonparametric kernel regression.While re-peated iterations are possible to further minimize the energy
functional͑2͒,in practice,only a single pass of the bilateral filter is generally needed to obtain substantial denoising re-sults.Furthermore,a separable approximation of the bilateral filter,as proposed by Pham and van Vliet,18is available and allows for a dramatic increase in computational performance with little to no degradation in efficacy.This approximation consists simply of running a one-dimensionalfilter in one directionfirst,then applying thefilter to the results of this step in the other direction,analogous to the way FFTs are usually calculated.For all applications in this work,the sepa-rable approximation is employed.
II.B.CT noise model
Denoising with bilateralfiltering can be directly applied in image space after reconstruction,and indeed the sharper edges there are well suited to techniques like bilateralfilter-ing,but the noise model in image space is very complex due to the data preprocessing and image reconstruction.In co
n-trast,noise in projection space is relatively easier to model under certain simplified assumptions.In this study,we will focus only on the application of bilateralfiltering to projec-tion space denoising.
Noise in CT images originates from data noise in the pro-jection measurement,which has two principal sources: Quantum noise and electronic noise.The electronic noise is the result of electronicfluctuation in the detector photodiode and other electronic components.The quantum noise is due to the limited number of photons collected by the detector. Although a current CT detector is not a photon-counting el-ement but an energy integrator that generates a signal pro-portional to the total energy deposited in the detector,a photon-counting model is still a good approximation and is widely used for characterizing noise properties of the CT data.10,19,20More accurate noise models have been investi-gated,such as the compound Poisson model that takes into account the polychromatic x-ray beam and energy integration.21As explained therein,the actual residual error introduced by a photon-counting model is only a few percent for typical photonflux level in clinical CT protocols.As also indicated in that work,the impact from the bowtiefilter and tube current modulation on noise characteristics of CT data
is even more significant.Therefore,for simplicity,we use a
photon-counting model and will consider the effect of the
bowtiefilter and tube current modulation in this work.
For a given attenuating path in the imaged subject,denote
the incident and the penetrated photon numbers as
N0͑␣,,͒and N͑␣,,͒,respectively,where␣and v de-note the index of detector bins along the radial and longitu-
dinal directions,andthe index of projection angle.In the
presence of noise,the measured data should be considered as
a stochastic process.Ideally,the line integral along the at-
tenuating path is given by P=−ln͑N/N0͒.Herein,we neglect the detector index͑␣,,͒in all the variables unless explic-itly indicated.
We assume that N0is a deterministic constant and N is
Poisson distributed with mean N.We also neglect the elec-
tronic noise and assume that the data collected on each de-
tector bin are uncorrelated.The number of input photons
N0͑␣,,͒must now be estimated as a function of␣,,and .
II.C.Incorporation of the effect of x-ray beam bowtie filter
The incident number of photons varies for each projection angle due to the use of the automatic exposure control͑AEC͒technique22and is also nonuniform across the radiationfield of the x-ray beam mainly due to the use of the beam-shaping bowtiefilter.23To accurately quantify the noise properties in the projection data and preserve the noise pattern in the de-noised image,these effects must be taken into consideration. This task can simply be accomplished by expressing the in-cident number of photons N0as a function of detector bin index␣and v based on the estimation of the nonuniformity across the x-ray radiationfield and a function of projection anglebased on the estimation of the tube current modula-tion.
As shown in Fig.1͑left͒,a bowtiefilter is usually used in the x-ray beam.Because the cross section of most patients is ovally or circularly shaped,the attenuation in the peripheral region is much less than that in central region.The
purpose F IG.1.Effect of x-ray beam bowtiefilter.Left:Illustration of the bowtie filter in the x-ray beam to reduce the incident x-ray intensity in the periph-eral region of the x-ray fan-beam.Right:An example of noise-equivalent number of x-ray quanta curves for140,120,100,and80kV that can be used for characterizing the effect of the nonuniform photon intensity caused by the bowtiefilter.
Medical Physics,Vol.36,No.11,November2009
of the bowtie filter is to reduce the incident x-ray intensity in the peripheral region so that the radiation dose to the patient,especially the skin dose,can be minimized.As a conse-quence,the x-ray intensity incident to the patient is highly nonuniform across the fan beam,which affects the noise properties in the measured CT data.We empirically deter-mine the effect of the bowtie filter by measuring the variance of the transmission from an air scan.The inverse of the vari-ance is a number that can be used to estimate the incident x-ray intensity across the x-ray beam.21Figure 1͑right ͒dis-plays an example of the noise-equivalent incident quanta number on a single detector row across the x-ray beam ob-tained for different tube potentials ͑kV ͒.By fitting the vari-ance with a third-or fourth-order Gaussian equation,four calibration curves of the incident x-ray intensity for the kVs at 80,100,120,and 140are obtained.Such calibration curves may be different,depending on the configuration of the detector collimation.
II.D.Incorporation of AEC
AEC is widely used for dose reduction in abdominal CT.Figure 2shows one example of tube current modulation from a CT scanner.The curve represents the reference signals of the tube as a function of the table position,which corre-sponds to a certain tube angle.The reference signal is pro-portional to the tube current used for that table position and projection angle.As can be seen,the tube current
oscillates during the gantry rotation in order to adapt to the attenuation level of the patient along different orientations.This auto-matic tube current modulation leads to a continuous change in the incident x-ray intensity,which will also affect the noise characteristics of the CT data.We incorporate this ef-fect by extracting the reference signal from each projection
frame and then estimating the corresponding reference x-ray intensity.The calibration curves determined from the bowtie filter are used for this estimation.
II.E.Sinogram smoothing with bilateral filtering
To incorporate the noise model in bilateral filtering,we first convert the measured sinogram P to a dataset represent-ing a map of detected number of photons,which is simply expressed as N =N 0exp ͑−P ͒.
͑7͒
Practically,the x-ray photons emitting from the tube are polychromatic.The total number of incident photons before the attenuation N 0can be estimated as the noise-equivalent photon number,as detailed above.The Anscombe transform 24is a standard statistical tool to convert Poisson distribute
d data to data with an approximately normal distri-bution with a constant variance,Q =2ͱ͑N +3/8͒.
͑8͒
This transformation is considered valid when the mean value of the Poisson data is greater than 20.The minimum flux levels in CT data are at least several hundreds in all the datasets we have tested thus far.Here we also ignore the small constant additive factor of 3/8͑negligible at these count numbers,but easily included if desired ͒and the mul-tiplicative factor of 2,and simply take the square root of the photon number to generate the dataset which is subsequently used for denoising,with a constant standard deviation of 1/2.The square-root-transformed data are given by Q =ͱN =ͱN 0exp ͑−P ͒.
͑9͒
Equation ͑5͒can be applied on Q to obtain the denoised
dataset Q
ˆ.The steps are then reversed to convert Q ˆback to the log-transformed sinogram,which can be used for regular image reconstruction.
II.F.Evaluation of noise and resolution properties
We performed phantom studies to evaluate the noise-resolution properties of bilateral filtering with CT noise mod-eling.These studies were performed on a dual-source CT scanner ͑Somatom Definition,Siemens Medical Solutions,Forchheim,Germany ͒.A phantom with a small acrylic cyl-inder and a thin wire that is typically used for quality control on Siemens scanners was scanned with the following param-eters:120kVp,0.5s rotation time,32ϫ0.6mm 2detector collimation,672ϫ32detector matrix size,helical pitch of 1.0,and mA s values of 40and 80.The images were recon-structed using a series of kernels available on the scanner with a slice thickness of 1mm and an FOV size of 50mm.The raw data were downloaded from the scanner and de-noised with bilateral filtering.The Gaussian filter with a given spatial standard deviation d in Eq.͑3͒was approxi-mated to a given filter length w ,and we found it useful empirically to lock the ratio of d to w as 1/6.The denoising parameters were filter length w =5͑corresponding to
d
F I
G .2.Effect of AEC.Displayed is an example of tube current modulation from an abdominal CT exam.The curve represents the reference signals ͑I 0͒on the detector as a function of table position.The reference signal is pro-portional to the tube current used for that table position and projection angle.The tube current oscillates as the table translates to adapt to the different attenuation levels of the patient along different projection angles.
Medical Physics,Vol.36,No.11,November 2009
=5/6and evaluated fromϪ2to2͒andϭ0.7,1,1.4,1.8, 2.2,and2.8.We also explored different values of w and found,as expected,that largerfilters gave lower noise but lower resolution.For space reasons,we only discuss results with wfixed at5,again determined empirically.The de-noised data were uploaded to the scanner and reconstructed with both the B70f and B40f kernels͑the B70f kernel is an extremely sharp but noisy kernel,while the B40f is slightly sharper and noisier than is used at some institutions,but it is employed in our clinical practice for routine body CT͒.The modulation transfer function͑MTF͒on each image,either directly reconstructed on the scanner or after bilateralfilter denoising,was calculated using a method similar to that in Ref.25.The noise level at a small region of interest͑ROI͒15 mm away from the wire was also measured.
A second phantom consisted of a stadium-shaped water tank͑lateral width of30cm,height of22.5cm͒.Eight thin plates with different contrast levels͑120HU:2;70HU:3; 40HU:3͒were placed in the central region of the water tank to allow visualization of these lower contrast materials.The scanning parameters were detector collimation64ϫ0.6with a z-flying focal spot,120quality reference mA s,rotation time of0.5s͑the term“quality reference mA s”is used in the automatic exposure control software in Siemens scanners to prescribe mA s.Its value equals the effective mA s used for a standard patient size͒.The images were reconstructed with commercial kernels at B45,B40,B30,B20,and B10 with a slice thickness of0.6mm.The raw data were de-noised with w=5andϭ1and reconstructed with various kernels.
II.G.Evaluation in CT colonography
We performed a preliminary study to evaluate the noise-
resolution properties of bilateralfiltering.The tested CT raw
data were obtained from a patient study undergoing a nonca-
thartic CT colonography.The patient was scanned on the
adaptive
dual-source CT system with the following scanning param-
eters:120kVp,100quality reference mA s,0.5s rotation
time,32ϫ0.6mm2detector collimation,and a helical pitch
of1.4.The volume CT dose index͑CTDI vol͒is7mGy.The images were reconstructed on the scanner with a slice thick-
ness of1mm and several reconstruction kernels͑B10f,B20f,
B30f,and B40f͒that represent a trade-off between noise and
resolution.The raw data were downloaded from the CT scan-
ner and processed with bilateralfiltering with different pa-
rameter settings.The smoothed data were subsequently re-
loaded on the scanner and reconstructed with a B40f kernel,
and the noise and resolution properties were analyzed.
II.H.Evaluation in CT enterography
We applied bilateralfiltering tofive abdominal CT exams
performed on the dual-source scanner.Each exam was
scanned in a dual-energy mode with one tube operated at140
kV and the other at80kV,with one-half or less of the power
distributed to the80kV tube.The80kV data is acquired at
425quality reference mA s.Two sets of images͑low and
high energy͒were generated from each dual-energy scan.Besides the dual-energy processing that is used to provide material-specific information,26–28the dual-energy images are mixed together in a linear fashion to form a single set of images,which serves a similar purpose as the images from a conventional single-energy scan at120kV.
The80kV images obtained from a dual-energy scan have enhanced iodine signal compared to140kV͑t
he iodine sig-nal is approximately doubled͒.However,they are often sub-ject to increased noise contamination due to the higher at-tenuation of the low-energy x rays.Low-energy images alone are often insufficient for diagnostic purposes despite the im-proved contrast enhancement.We applied bilateralfiltering ͑w=5,ϭ1͒,followed by B40f reconstruction to the80kV images and compared these images to B40f and B20f recon-struction alone.
III.RESULTS
III.A.Noise-resolution properties
Figure3compares the noise-resolution properties of bilat-eralfiltering with the body kernels available on the scanner on the thin wire phantom.The noise was represented by the standard deviation͑HU͒in a small ROI15mm away from the thin wire in the phantom.The spatial resolution was rep-resented by the spatial frequency at10%of the maximum value on the MTF curve.The solid curve linking the
solid F IG.3.Noise-resolution properties of bilateralfiltering and body kernels. Noise is expressed as the standard deviation͑HU͒in a small ROI close to the thin wire.The spatial resolution was quantified a
s the spatial frequency at10%of the maximum value on the MTF curve.The solid curve linking the solid triangles͑᭡͒was obtained from the image scanned with40mA s and reconstructed with kernels of B10f,B20f,B30f,B40f,and B50f͑from left to right͒.The solid diamond symbol͑ࡗ͒was from the same scan and represents a special body kernel B25f.The dashed curve linking the open triangles͑᭝͒represents the noise-resolution results obtained from images after applying bilateralfiltering to the40mA s scan data with ten different smoothing parameters.From left to right,thefirstfive points were recon-structed with the B40f kernel,with wfixed at5,andϭ2.2,1.8,1.4,1.0, 0.7,respectively;the secondfive points were reconstructed with a sharper B70f kernel,with wfixed at5,andϭ2.8,2.2,1.8,1.4,and1.0,respec-tively.The solid curve linking the solid circles͑b͒was obtained from the image scanned with80mA s and reconstructed with kernels of B10f,B20f, B30f,B40f,and B50f͑from left to right͒.The dashed curve linking the open circles͑᭺͒represents the noise-resolution results obtained from images after applying bilateralfiltering to the80mA s scan data with the same ten smoothing parameters as above.The noise-resolution results on the40mA s data after bilateralfiltering approach or exceed the noise-resolution proper-ties of the80mA s data,andfiltering is effective on the80mA s data as well as on the40mA s data.
Medical Physics,Vol.36,No.11,November2009
版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。
发表评论