Physics in Medicine & Biology

Paper • The following article is Open access

Bilateral filtering using the full noise covariance matrix applied to x-ray phase-contrast computed tomography

S Allner¹, T Koehler^2,3, A Fehringer¹, L Birnbacher¹, M Willner¹, F Pfeiffer^1,4 and P B Noël^1,4

Published 22 April 2016 • © 2016 Institute of Physics and Engineering in Medicine
Physics in Medicine & Biology, Volume 61, Number 10 Citation S Allner et al 2016 Phys. Med. Biol. 61 3867 DOI 10.1088/0031-9155/61/10/3867

Download Article PDF

Article metrics

2461 Total downloads
Video abstract views

Author e-mails

sebastian.allner@tum.de

Author affiliations

¹ Department of Physics and Institute of Medical Engineering, Technische Universität München, Garching, Germany

² Philips Research, Hamburg, Germany

³ Institute for Advanced Study, Technische Universität München, Munich, Germany

⁴ Department of Radiology, Klinikum rechts der Isar, Technische Universität München, Munich, Germany

Dates

Received 26 August 2015
Accepted 29 February 2016
Published 22 April 2016

Buy this article in print

Journal RSS

Abstract

The purpose of this work is to develop an image-based de-noising algorithm that exploits complementary information and noise statistics from multi-modal images, as they emerge in x-ray tomography techniques, for instance grating-based phase-contrast CT and spectral CT.

Among the noise reduction methods, image-based de-noising is one popular approach and the so-called bilateral filter is a well known algorithm for edge-preserving filtering. We developed a generalization of the bilateral filter for the case where the imaging system provides two or more perfectly aligned images. The proposed generalization is statistically motivated and takes the full second order noise statistics of these images into account. In particular, it includes a noise correlation between the images and spatial noise correlation within the same image.

The novel generalized three-dimensional bilateral filter is applied to the attenuation and phase images created with filtered backprojection reconstructions from grating-based phase-contrast tomography. In comparison to established bilateral filters, we obtain improved noise reduction and at the same time a better preservation of edges in the images on the examples of a simulated soft-tissue phantom, a human cerebellum and a human artery sample. The applied full noise covariance is determined via cross-correlation of the image noise.

The filter results yield an improved feature recovery based on enhanced noise suppression and edge preservation as shown here on the example of attenuation and phase images captured with grating-based phase-contrast computed tomography. This is supported by quantitative image analysis. Without being bound to phase-contrast imaging, this generalized filter is applicable to any kind of noise-afflicted image data with or without noise correlation. Therefore, it can be utilized in various imaging applications and fields.

Export citation and abstract BibTeX RIS

Previous article in issue

Next article in issue

Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

1. Introduction

Image noise is one of the major factors which can significantly reduce the diagnostic quality in computed tomography (CT). As it is bound to the radiation dose, it can be decreased by increasing the radiation exposure. Consequently, higher-quality CT images normally imply a higher radiation dose applied to the patient, which can be even higher for multi-modality CT like dual-energy CT. As extended radiation exposure increases the chance of a radiation induced cancer (Brenner and Hall 2007), there is a need for tools to reduce the noise level in images while preserving diagnostic information.

One way to improve image quality and to reduce noise is to use iterative reconstruction techniques (Koehler et al 2011a, Gaass et al 2012, Noël et al 2013, Hahn et al 2015). Those techniques have proven high dose reduction potentials in the clinical arena (Noël et al 2013). However, those iterative approaches require increased computation time compared to analytical reconstructions. A computationally more attractive solution is to use an image de-nosing algorithm. For classical absorption-based x-ray CT several methods have been proposed to de-noise in projection (Manduca et al 2009) or in image space (Kalra et al 2003, Bruder et al 2009). The simplest way to achieve reduced noise in images is low-pass filtering. This however leads to a loss in spatial resolution and thus often to a loss of diagnostic quality. More advanced image-based methods fall into the class of so-called edge-preserving de-noising (Kervrann and Boulanger 2006, Chatterjee and Milanfar 2012). The general idea behind this type of algorithm is to detect (implicitly or explicitly) edges and to avoid low-pass filtering across these edges. The bilateral filter (Tomasi and Manduchi 1998) is one popular member of this class. Standard bilateral filtering and its extension for two aligned images exploiting their complementary edge information (Koehler and Roessl 2012) assume white noise in the image to be filtered.

One limitation in absorption imaging is the availability of only a single signal. However, new technologies, such as spectral photon-counting (Roessl and Proksa 2006, 2007) or grating-based phase-contrast CT (Pfeiffer et al 2007b, Bravin et al 2013), are on the horizon and offer multiple signal bands from a single data acquisition. Especially phase-contrast imaging, which provides absorption, phase, and scattering information, has undergone a remarkable development (Schleede et al 2012, Grandl et al 2014, Herzen et al 2014). The enhanced soft-tissue contrast is particularly useful in computed tomography and yields complementary data to standard absorption imaging (Pfeiffer et al 2007a, Herzen et al 2009). Like in all imaging systems, the reconstructed phase-contrast images suffer from noise.

In this work, we explore the possibility to leverage all available information to generate fast and reliable data with reduced noise. We apply a more general approach where a possible noise correlation is taken into account. For phase-contrast CT, the attenuation image has predominantly high-frequency noise, whereas the phase image suffers from lower-frequency noise exhibiting a high spatial correlation. The goal of this work is to include this additional information into the bilateral filter. We also investigate the advantages of three-dimensional (3D) over two-dimensional (2D) filtering. Furthermore, we show the influence of the noise covariance on the weighting of neighboring image voxels. Finally, the capability of the proposed algorithm to reduce noise and to preserve diagnostic quality is illustrated on simulated data for a multi-inlay soft-tissue phantom and on experimental phase-contrast reconstructions of a human cerebellum as well as a human heart artery sample.

2. Materials and methods

2.1. Bilateral filtering

The concept of bilateral filtering (Tomasi and Manduchi 1998) is based on a selective weighting scheme for averaging neighboring pixels to achieve a de-noised image. The general form of a bilateral filter contains a distance-dependent domain filter part $d\left(x,{{x}^{\prime}}\right)$ and a gray-value-dependent range filter part $r\left(\,f(x),f\left({{x}^{\prime}}\right)\right)$ :

with x being the position of the central pixel, ${{x}^{\prime}}$ the positions of neighboring pixels, and N(x) a normalization factor. While the domain filter accounts for local weighting of neighboring pixels the range filter part enforces the value-dependent component to prevent filtering across edges. For domain and range filter part normally a Gaussian function is used. It depends either on the Euclidean pixel distance or on the pixel value difference:

with ${{\sigma}_{d}}$ being a width parameter of the filter kernel size and ${{\sigma}_{f}}$ the noise standard deviation of the considered reconstruction values (e.g. attenuation noise standard deviation ${{\sigma}_{a}}$ ). The Gaussian weighting of gray value differences implies the model of white or uncorrelated noise.

In image processing of colored images three color bands have to be treated. In most cases they are saved as red, green and blue band (rgb). When they are filtered separately, the colors at the edges of features can be corrupted because the contrast levels in the different bands can be very different. Therefore, the transition of colors can vary drastically between bands which creates color seams at edges. This can be avoided by transferring the images into the CIE-lab color space which is related to the human perception (Tomasi and Manduchi 1998). There, combined filtering of the image bands produces best results because small changes in color values corresponds to small changes of color in visual perception.

In contrast to de-noising colored images, where noise is assumed to be white and uncorrelated in each band, different modalities x-ray imaging show different noise structures. Grating-based phase-contrast computed tomography shows high-frequency noise in the attenuation band and low-frequency noise in the phase band (Koehler et al 2011b). In this work, we include the spatial noise distribution in the extended filter approach into the weighting process. The different aligned signals of grating-based phase-contrast imaging are from this point forward referred to as bands.

The class of bilateral filters applies a selective weighting scheme according to which a voxel is averaged with its neighbors to obtain the output image. The relative weight for each neighboring voxel depends on the spatial distance and gray-value distance to the central voxel. Since a voxel has two bands we consider in the following two voxel positions and two bands to derive the relative weight for the averaging process.

This filter approach is developed on the example of grating-based phase-contrast imaging and therefore we refer to the volumetric image gray values of the bands by a for attenuation and p for phase. However, it is not limited to this particular imaging technique and can be applied to any registered multi-band data with or without intra-band (spatial) and/or inter-band noise correlation.

The general filter formula is derived by investigating the probability of two pixels having a certain gray value distance under the influences of noise. Given a vector of gray values $\left({{a}_{1}},{{a}_{2}},{{p}_{1}},{{p}_{2}}\right)$ with their respective noise realizations $\left({{n}_{1}},{{n}_{2}},{{n}_{3}},{{n}_{4}}\right)$ , the following Gaussian model describes the noise distribution in two voxels in two bands:

where n₁ and n₂ are the noise values of the two considered voxels in one band, n₃ and n₄ the noise of the same voxels in the other complementary band. The matrix S is the inverse of the noise covariance matrix of the two two-band voxels. In the case of two image bands the covariance matrix C between two voxels is a $4\times 4$ matrix. Symmetry arguments and the assumption of stationary noise lead to:

with ${{s}_{kl}}={{(S)}_{kl}}$ being the elements of the inverted covariance matrix. Since the variables are in the order $\left({{a}_{i}},{{a}_{j}},{{p}_{i}},{{p}_{j}}\right)$ with i being the treated central voxel and j the corresponding voxels of the weighting array the upper left $2\times 2$ sub-matrix represents the spatial noise correlation of the first band and the bottom right $2\times 2$ sub-matrix the spatial noise correlation of the other band. The other matrix entries relate to the linkage between the two bands.

The idea for the proposed range filter is to use the marginal likelihood that a certain difference is observed in the two bands, given the assumption that the true values of the two voxels are the same in each band, respectively. Therefore, we introduce the difference of the gray values in one band as $\Delta a$ and $\Delta p$ . Given the assumption that the true values of the two voxels are the same, these differences are completely due to noise, i.e. ${{n}_{2}}=\Delta a-{{n}_{1}}$ and ${{n}_{4}}=\Delta p-{{n}_{3}}$ . After inserting these into (4), the marginal likelihood is calculated by integrating over all possible noise realizations:

Exploiting the symmetry of S leads to:

This equation represents the generalized range filter which takes the noise correlation into account by utilizing the inverse covariance matrix. For the case of completely uncorrelated noise, the covariance matrix becomes diagonal C = c_ii and therefore, also all off-diagonal elements of S vanish. As the diagonal elements of the inverse covariance matrix of uncorrelated noise are just the inverse variances of the noise ( ${{s}_{ii}}=\frac{1}{\sigma _{i}^{2}}$ ), the range filter boils down to:

Thus, this cuts down to the multi-band range filter (Koehler and Roessl 2012) except for the additional factor of two in the denominator of the exponent. It can be explained by the statistically motivated derivation of the filter presented here. The variance of the differences $\Delta a$ and $\Delta p$ needs to be used, which is twice the variance of the noise in the band. It also can be explained from a different point of view. In a normal Gaussian distribution the statistically varying quantity would be subtracted from the mean value of the distribution. Therefore, the statistical properties of the gray values apply. Here, we do not subtract a mean value, but a different voxel value being also afflicted by statistical fluctuations resulting in a doubled variance of the distribution.

2.2. Investigated range filters

We applied the filters to attenuation and phase reconstructions after filtered backprojection (FBP). The following equations show the filtered attenuation $\tilde{a}$ band. The formulation of the corresponding filtered phase $\tilde{p}$ is straightforward. The investigated filter based on equation (7) is referred to as covariance-based filter which takes additional information from the complementary band and noise correlations of the datasets into account:

where x_j is the considered voxel in the neighborhood ${{\mathcal{N}}_{i}}$ of the central voxel x_i. For the example of grating-based phase-contrast imaging, the mixed term vanishes in the filter formula as attenuation and phase noise are not correlated (Weber et al 2011). However, for other imaging modalities with inter-band correlation this might not be negligible. We compared this filter to the standard Gaussian bilateral filter, where the image bands are treated separately. The single-band filter is defined as:

The multi-band filter (Koehler and Roessl 2012) uses a weighting procedure deduced from a common range factor derived simultaneously in both bands. It can be written as:

It uses the edge information from both bands to recover more reliable weights but does not consider noise correlation.

2.3. Determination of the noise covariance

To derive the full noise covariance matrix, a common approach is evaluating noise in a region that is known to contain no signal information. The fundamental assumptions for this approach is that the noise in the images is at least approximately stationary (Kak and Slaney 2001). In other words the noise properties do not vary across the image. This hypothesis can be verified for images generated with FBP, keeping in mind that we do not suffer from photon starvation. Additionally, the noise variance in the projections of our measurements, and therewith in the image volume, is almost constant across the projections as all objects were submerged in a water bath during measurement and do not contain highly absorbing materials.

The covariance analysis in Fourier domain is based on the Wiener–Khinchin theorem (Cohen 1998) for stationary processes as considered here as cross-correlation of the sample-free region of the different bands. For noise analysis, the mean value was subtracted for each band and zero padding to the doubled size was performed. The 3D spatial correlation was then calculated according to:

with A being the 3D spatial correlation of the background volume B ( ${{B}^{\prime}}$ zero padded) between band k and l. It is normalized by the number of voxels in the background region N_B and $\mathcal{F}$ represents the 3D Fourier transform. For k = l the spatial correlation of the noise region within each band is called the autocorrelation. The main diagonal elements of the covariance matrices (see equation (5)) correspond to the central order of the respective autocorrelation (variance). The spatial correlation corresponds to the remaining elements of the covariance matrix. For grating-based phase-contrast imaging the inter-band correlation between attenuation and phase A_kl vanishes for $k\ne l$ as the attenuation and phase noise do not influence each other (Weber et al 2011). This way the covariance matrix is filled with the respective spatial autocorrelation depending on the distance to the central voxel for each voxel in the weighting array. The axial correlation is expected to be isotropic with a short-range correlation for the attenuation band. The phase image provides long-range linkage due to the low-frequency noise originating from the phase integration during reconstruction. Figure 1 exhibits the spatial correlation of noise in both investigated signal bands.

Figure 1. Refer to the following caption and surrounding text. — **Figure 1.** The axial and transversal spatial correlation of the attenuation (left side) and the phase (right side) noise. They were acquired from a sample-free region of the human artery sample. The correlation is shown in a correlation range of [−0.1, 1].
Download figure:
Standard image High-resolution image

As theoretically expected the axial correlation of attenuation is of short range whereas the the corresponding phase investigations show long range correlation. Additionally, the FBP algorithm would suggest an uncorrelated noise structure in other slices because of its slice-independent nature. We observed a significant coronal correlation of short range comparable in attenuation and phase which could not be explained by the projection geometry. It stems from a correlation between adjacent detector pixels. This originates from charge sharing of the photon-counting detector. The effect is caused by a photon with at least the doubled detector threshold energy triggering two pixels instead of one. This charge sharing correlation of 1 px range is isotropic in projection space which pursues into the reconstruction via FBP and causes the observed vertical linkage. Note, that this effect is consequently also considered in the weighting of the covariance-filter.

2.4. Tools for quantitative analysis

2.4.1. Contrast-to-noise ratio.

A very basic attempt to quantify the quality of noisy images is the contrast-to-noise ratio of homogeneous regions. This well known technique compares the value difference of two regions of interest to the corresponding standard deviation of the noise:

where ${{\bar{a}}_{i}}$ are the mean gray values within regions 1 and 2. The $\sigma _{i}^{2}$ correspond to the respective noise variance. It yields reliable information about noise suppression within flat regions of the image. Therefore, a high $\text{CNR}$ is desirable. However, the behavior of de-noising filters at edges or structures cannot be characterized by CNR.

2.4.2. Structural similarity index.

A popular way to characterize image quality, if there is a reference image it can be compared to, is the structural similarity index. It attempts to quantify the image impression for human perception with respect to the reference by means of luminance, contrast, and structure (Wang et al 2004, Wang and Bovik 2009). It is used in this work as:

The image is characterized via the local mean value μ and the local standard deviation σ (within a certain neighborhood of the image). In this formulation f refers to the evaluated and r to the reference image. Here, ${{\sigma}_{fr}}$ corresponds to the local cross-correlation of these two images. The constants ${{c}_{\text{L}}},{{c}_{\text{C}}},{{c}_{\text{S}}}$ ensure numerical stability. They are small and set according to Wang et al (2004). For perfectly matching image and reference the structural similarity is equal to one. Any difference between the two images results in a decreased similarity. Therefore, the image quality is superior for higher the structural similarity. However, for this evaluation a reference image is mandatory which is often not available.

2.4.3. Histogram entropy.

The histogram or image entropy is a quantitative analysis tool which utilizes the gray value distribution of a given image for quality assessment. It originates from the basis of information theory (Shannon and Weaver 1949):

where p_k is the probability density function of the image and K being a normalization constant. Here, we apply an estimated probability density function (Parzen window) and normalize it to the maximum entropy for the numerical evaluation (Donath et al 2006):

with $\Delta g$ being the bar width and n_j being the number of gray values contributing to the estimated probability density function. The width of the Parzen window is h and the total number of contributing pixels is N. An image without any feature information is just white noise which results in the maximum entropy H_max. Therefore, the image contains more feature information from an information theory point of view for lower histogram entropy. The histogram entropy can also be used for optimization purposes (Schüller et al 2015). Therefore it can be regarded as image quality metric.

2.5. Simulation of a soft-tissue phantom

For characterizing qualitative and quantitative performance of the proposed bilateral filter, a soft-tissue phantom with different material inlays is used (Koehler et al 2011b). It is a mathematical phantom using inlays with the characteristic values of human tissue in attenuation and phase.

2.6. Scanning parameters

The filters were applied to image data of a human artery and a cerebellum sample which were acquired at a symmetric grating-based Talbot-Lau interferometer lab setup with a rotation stage, a fixed rotating-anode source and a PILATUS II photon-counting detector. All gratings have a period of 5.4 μm and an inter-grating distance of 80 cm. The distance between source and detector is 137 cm and the distance between source and rotation axis is 99 cm which translates into almost parallel beam geometry. To avoid phase wrapping the formalin-fixated samples were put into plastic cylinders (3 cm in diameter) and imaged in a water bath. The x-ray tube was operated at 40 kVp resulting in an effective energy of ≈27 keV. Both samples were measured with 11 steps for each of the 1200 projections. The heart artery was imaged at different exposure times per step (long exposure: 3.6 s, medium exposure: 0.4 s, and short exposure: 0.144 s). The cerebellum sample was measured with only a long exposure of 5.0 s per step. From those stepping images absorption and differential phase projection were acquired (Weitkamp et al 2005).

3. Results

Here, we want to assess the performance of the novel covariance-based bilateral filter. To achieve optimal filter results, we first investigate 2D slice-based and 3D volumetric filtering in section 3.1. The influence of the extended weighting scheme is also exhibited here (see section 3.2). It depicts how the additional information changes the voxel weights. In the following sections 3.3–3.5 the filter performance is evaluated on the examples of a simulated soft-tissue phantom, a cerebellum, and heart artery sample.

Throughout all filters and samples the domain filter has a width of ${{\sigma}_{d}}=4$ px. The uncertainties of the attenuation and phase signal ${{\sigma}_{a}}$ and ${{\sigma}_{p}}$ for the single-band and multi-band filter were measured by calculating the standard deviation of gray values in a region known to contain no change in attenuation values, keeping in mind, that the noise in the image is approximately stationary. In this example, the regions contain only formalin.

3.1. Influence of 2D or 3D filtering

For a given range filter, the amount of noise reduction that can be achieved with bilateral filtering is influenced by the width of the Gaussian domain filter part given by the width parameter ${{\sigma}_{d}}$ . Of course, the wider the domain filter is, the more voxels contribute to the filtered output image and the better the noise reduction works. Another option for the implementation of the filter is whether to operate on two-dimensional (2D) images independently or to operate on the full three-dimensional (3D) volumes. For a given width parameter, many more pixels are used for averaging and better noise suppression is expected for 3D filtering. Another important aspect is the spatial correlation of the noise in the phase image. It is dominated by low-frequency noise, which cannot be removed efficiently by any low-pass filter of reasonable size, especially in 2D de-noising of axial slices. Since noise is almost uncorrelated across slices, the remaining low-frequency noise in axial images appears as horizontal streaks in coronal or sagittal views as shown below.

In figure 2, we illustrate the difference between 2D (axial) and 3D (volumetric) filtering. These single-band filtered images have a reduced noise level compared to the unfiltered FBP. The 3D filtered slice provides better results due to the larger number of voxels used in the domain filter. Additionally the differential phase projections cause low-frequency noise in axial slices and introduce strong noise streaks in the coronal slices as expected. These streaks are contained in 2D filtering but reduced in 3D de-noising. In all subsequent investigations, 3D filtering is preferred over 2D filtering because of its improved results.

Figure 2. Refer to the following caption and surrounding text. — **Figure 2.** The unfiltered, 2D and 3D single-band filter results of an axial and a coronal slice from phase FBP of the human artery sample. The noise level is reduced further in 3D filtering than in the 2D case. Additionally, the noise structure in coronal slices is partially overcome by 3D de-noising.
Download figure:
Standard image High-resolution image

3.2. Integrating complementary information

Pure domain filters like a simple Gaussian or mean filter smooth over edges and features of the approximate size of their kernel. Bilateral filtering prevents this by taking the gray value distance and the noise standard deviation into account. However, if an edge or feature is in the range of noise it gets washed out. This unwanted behavior is partially overcome by the multi-band or covariance-based filter because they include the combined information of both bands. If this edge or feature (being faint in one band) has a distinct value difference in the complementary image band, the two materials can be well separated there which translates into an improved overall weighting. To show the influence of the weighting scheme beyond the contribution of the Gaussian domain filter part, the range weights of the attenuation image obtained by the different bilateral filters are displayed in figure 3. For completeness the weights of the single-band covariance filter are also listed there but not pursued further in this work.

Figure 3. Refer to the following caption and surrounding text. — **Figure 3.** The weighting contribution of the range filter part for the different filter types in an axial image. The images to the left show unfiltered attenuation and phase of a human cerebellum sample with a detail zoom. A red box indicates the central voxel of the weighting array. The images to the right display the weight of the respective filter type for the attenuation image in the value range of [0, 1]. In general, the weights including covariances are higher than for the normal case resulting in a stronger smoothing. The red arrows highlight the voxels closest to the central voxel, where the weight is changed most due to incorporating the covariance information. This can be interpreted as weighting correlated voxels less.
Download figure:
Standard image High-resolution image

The attenuation detail does not show a distinct feature whereas there is an edge in the phase image. The weights in figure 3 follow the respective case in equation (19).

The region of interest shown in the second column in figure 3 contains an edge between gray and white matter of the brain sample as clearly visible in the phase image. However, it is almost completely hidden by noise in the attenuation image. Consequently, both range filter values of the single band filters (illustrated in the two images top right) would smooth significantly across the edge. The multi-band and covariance-based filter incorporate the additional image band and therefore apply more suiting weights. In addition, the weights of contributing voxels using the full noise covariance matrix are higher than the corresponding ones in the single- and multi-band filter. Finally, the voxels around the central voxel of the weighting array are less reliable because their noise exhibits a non-negligible correlation. The covariance-based filters account for that by decreasing the weight depending on the noise correlation as indicated by the red arrows in figure 3.

Using single-band, multi-band and covariance-based range weights, figure 4 shows a detail image view of the artery wall's attenuation. The feature border, namely the inner wall of the artery, is partially obstructed by noise in the unfiltered attenuation image in particular in the lower part of the image. It cannot be distinguished from the formalin since the attenuation of the formalin inside the vessel is very similar to the attenuation of the vessel wall. The single- and multi-band filter reduce the noise level by approximately the same amount, and the covariance-based filter decreases noise even further. However, on the inner edges of the artery wall the multi-band and covariance-based versions recover a more prominent edge than the single-band bilateral filter because the range filter part is influenced by the higher phase signal contrast. The best result is achieved by the covariance-based filter combining the smallest noise level and sharp edges.

Figure 4. Refer to the following caption and surrounding text. — **Figure 4.** A detail from an axial slice of the different filter results of the absorption signal of an artery wall. The image to the right shows the same detail in the phase image to give an idea about the information in the other band.
Download figure:
Standard image High-resolution image

3.3. Filtering results of the simulated soft-tissue phantom

Here, the results of the different filters by means of the simulated soft-tissue phantom with multiple inlays are characterized. In figure 5, the filtering results for this phantom are displayed for the established bilateral filters and the new covariance-based filter for the absorption and phase signal, respectively and compared to the unfiltered input volume. Visually, the covariance-based de-noising outperforms both established filters for both imaging modalities. The single- and multi-band filter produce comparable results. Only for one of the inlays in the attenuation images the multi-band filter yields superior outcome over the single-band method, shown in the detail feature of figure 5. There, complementary edge information introduced from the phase prevents from filtering over the inlay edge which is partially hidden in the noise.

Figure 5. Refer to the following caption and surrounding text. — **Figure 5.** The attenuation and phase reconstructions of the simulated phantom with results of bilateral filtering. The highlighted box shows the position of the detail feature.
Download figure:
Standard image High-resolution image

This visual image impression is supported by the quantitative analysis shown in table 1. Here, the contrast-to-noise ratio, the structural similiarity index with respect to a given reference image, and the histogram entropy are listed. As this is a simulated dataset the noiseless phantom is used as reference for the structural similarity index.

Table 1. Quantitative bilateral filter performance: contrast-to-noise ratio, structural similarity index and histogram entropy for simulated data.

		Unprocessed	Single-band	Multi-band	Covariance-based	Noiseless phantom
	CNR	10.88	21.33	21.19	31.71	—
Attenuation	SSIM	0.461	0.756	0.751	0.869	1.
	HE	0.346	0.283	0.284	0.247	0.073

	CNR	6.92	13.44	13.35	19.48	—
Phase	SSIM	0.844	0.948	0.947	0.972	1.
	HE	0.254	0.195	0.195	0.162	0.081

Note: The contrast-to-noise ratio is measured between the material in which the inlays are inserted and a background region, as the inlays are too small to give a reliable result. For the calculation of the structural similarity as well as histogram entropy the background regions outside the main cylinder are excluded.

The contrast-to-noise ratio increases from the unprocessed input to the single- and multi-band results to the covariance-based filter. However, the CNR is not calculated for the noiseless phantom as the ideal phantom does not have any noise. The structural similarity is calculated with respect to the noiseless phantom reference for all filter results. The increasing structural similarity index entailing improved de-noising performance can be seen for the attenuation as well as for the phase image. However, the SSIM cannot be interpreted as an absolute figure of merit as it is influenced by constants which enforce numerical stability. Therefore, only a relative comparison within one imaging modality is possible. The structural similarity for the noiseless phantom is perfect ( $\text{SSIM}=1.$ ) meaning total similarity because it is compared with itself. A low histogram entropy suggests reduced noise and sharp edges which result in sharper histogram peaks. Therefore, low entropy indicates better image quality. All numbers identify the covariance-based filter as the best of all investigated filters. For this phantom the single-band filter slightly outperforms the multi-band filter for each of the investigated quantitative values. This is due to the fact, that during the weighting process of the multi-band filter an additional weight of $\leqslant 1$ is multiplied to all contributing gray values. For the central voxel this weight is equal to one because $\Delta {{a}_{ii}}=\Delta {{p}_{ii}}=0$ . This increases the relative weight of the central voxel which slightly decreases the smoothing of the filter. However, this effect is negligible for reasonable sizes of the weighting array and the multi-band filter gains improved results at faint edges from incorporating complementary information as shown in figure 4.

3.4. Filtering results of the cerebellum sample

A comprehensive comparison of the single-band, multi-band, and covariance-based filter is conducted on the cerebellum sample from which axial and coronal slices are shown in figure 6. Obviously, the noise level in both attenuation and phase is reduced by all bilateral filters. Single-band and multi-band filters have approximately the same noise level while the multi-band filter provides better edge preservation because it considers the additional image band. The best results are achieved with the covariance-based filter since the noise level is reduced even further without smoothing over the edges. This is especially noticeable in the attenuation images where the feature borders often vanish in the noise of the unprocessed reconstruction because of the low soft tissue contrast of this signal. In the coronal slices of the phase the 3D filtering helps to partially overcome the horizontally oriented noise structure. The covariance-based filter yields again the best image quality and feature recovery through suppressing the noise.

Figure 6. Refer to the following caption and surrounding text. — **Figure 6.** The attenuation and phase reconstructions of the cerebellum sample compared with results of bilateral filtering.
Download figure:
Standard image High-resolution image

3.5. Filtering results for different image qualities of the heart artery sample

To perform quantitative analysis alongside qualitative investigations, a human heart artery sample was measured at different exposure times which is shown in figure 7 as the second test case for the different bilateral filter types. Since the sample was measured at different exposure times, the unfiltered long exposures can serve as reference for evaluating the de-noising. In addition, the filtered low-quality reconstructions can be visually compared to the unfiltered high-exposure references.

Figure 7. Refer to the following caption and surrounding text. — **Figure 7.** The attenuation and phase reconstructions of the artery sample of good, medium, and poor image quality compared with results of bilateral filtering.
Download figure:
Standard image High-resolution image

All filters reduce the noise level in the images with the best result being achieved by the covariance-based filter. The filters which use more bands achieve better results via an improved weighting. This is most apparent at the artery walls in the attenuation images where the complementary phase band saves the filter from smoothing over faint edges. Additionally, the edges appear less jagged in the attenuation with the multi-band and even a little less with the covariance-based filter image. This effect can be observed best in the medium exposure attenuation or the short exposure phase image at the top border of the specimen. Furthermore, the attenuation benefits more from the phase image in this example than the other way around. This can be explained by the improved signal contrast of the phase images and by the high phase sensitivity of the setup.

The quantitative analysis for the experimental results confirms the qualitative and descriptive investigations. Here, the performance of the proposed algorithm is characterized via quantitative figures of merit which are given in table 2. The contrast-to-noise ratio, the structural similarity index, and the histogram entropy of the covariance-based bilateral filter results and are compared to the outcome of established bilateral filters and unprocessed reconstructions.

Table 2. Quantitative bilateral filter performance: contrast-to-noise ratio, structural similarity index and histogram entropy for human heart artery sample.

		Unprocessed	Single-band	Multi-band	Covariance-based	Long exposure FBP
	CNR	1.07	2.06	2.07	3.08	5.39
Attenuation	SSIM	0.425	0.579	0.603	0.690	1.
	HE	0.409	0.382	0.382	0.365	0.344

	CNR	9.97	18.16	18.13	23.99	48.22
Phase	SSIM	0.623	0.748	0.752	0.797	1.
	HE	0.403	0.383	0.383	0.375	0.358

Note: The short exposure scan is filtered and the long exposure FBP serves as reference for all structural similarity calculations. The CNR is calculated from a formalin-containing region in the background and one region in the PMMA rod (bright circular region in phase FBP) in parts of the volumes not covered by the images. The histogram entropy and structural similarity are determined in the circular region inside the main plastic tube to exclude possible FBP artifacts in the outer areas.

In this example, the single-band filter has a slightly better CNR than the multi-band method only for the phase image. For the attenuation band the CNR of the multi-band filter probably is improved, because the selected region (PMMA rod) has a very little CNR (1.07) in the unprocessed image. Therefore, the edge is partially smoothed for the single-band case making the region not entirely flat. The higher structural similarity for the multi-band methods suggests superior de-noising because of the increased weighting due to incorporating the complementary band. However, the SSIM evaluation can be slightly biased as the reference image is the unfiltered long exposure FBP and not a noiseless image which entails a different noise realization. In addition, also the measurements might have been different due to minor instabilities between both exposures in the experimental setup. The covariance-based filter yields the best contrast-to-noise ratio, structural similarity, and histogram entropy for both image bands. This supports the previous findings from the simulated soft-tissue phantom.

4. Conclusion

We showed that the bilateral filter derived from a statistically motivated model reduces the noise level significantly in comparison to the classical formulation. Incorporating more bands in the weighting process and thus using complementary data decreases the possibility to smooth over faint features, especially when dealing with shorter exposure times. The covariance-based filter accounts for the noise correlation in the image bands and weights voxels with positively correlated noise less. For 3D filtering the effect of including off-diagonal elements of the inverse covariance matrix is smaller in our application, since only a few voxels in the same slice differ considerably in weighting compared to the large amount of uncorrelated weights as the noise is not correlated in different axial slices. For the combined filtering of attenuation and phase reconstructions the inter-band correlation is negligible because the two signals are independent from each other. For other applications however, the impact of using the covariance-based filter is expected to be even larger. This could be the case for instance in decomposed dual-energy images, where a strong noise correlation is present between the photo-electric and the Compton image (Morgan et al 1987, Kalender et al 1988, Gauntt and Barnes 1994).

Finally, the proposed extension to the bilateral filter reduces noise significantly and therefore has the potential to improve automated segmentation and help in medical diagnostics.

Acknowledgments

We want to thank A Hipp and K Achterhold for their support during setup construction. We also want to acknowledge Dr med. H Hetterich of the department of clinical radiology and Dr med. Christian Braun of the institute of forensic medicine, both part of the Ludwig-Maximilian-Universität München, which provided the samples for this work.

Please wait… references are loading.