Evaluation of noise removal algorithms for imaging and reconstruction of vascular networks using micro-CT

Under the assumption that a contrast-enhanced vascular μCT  image has an approximate piecewise constant structure, a TV model known also as Rudin–Osher–Fatemi (ROF) (Rudin et al 1992) can be applied. TV was reported in the literature as a noise reduction and deblurring method (Vogel and Oman 1998, Chan et al 2000, Beck and Teboulle 2009a, Oliveira et al 2009) and a large number of variants of this model have been proposed. The merit of this regularization model is that the edges are well-preserved in the denoised images. Recently, an image space reconstruction algorithm based on TV regularization was proposed to overcome the quantum noise limitation associated with low-dose μCT  (Vandeghinste et al 2013) in the in vivo case.

The constrained minimization problem considered in the ROF model is given by:

$$\begin{eqnarray}&&f(u)\\ &&=\mathrm{argmin}\left\{\displaystyle \frac{1}{2}{\parallel u-z\parallel }_{2}^{2}+\alpha {| u| }_{\mathrm{TV}}:u\in {{\mathbb{R}}}^{d}\right\},\end{eqnarray} \tag{ 1 }$$

where $u\in {{\mathbb{R}}}^{d}$ is the desired true solution computed at each iteration, z is the observed or noisy data, α > 0 is the regularization parameter and ${| u| }_{\mathrm{TV}}={\parallel {\rm{\nabla }}u\parallel }_{1}$ is the TV of u. Due to the non-differentiability of the TV norm and the large dimension of the underlying images, the functional from equation (1) is difficult to minimize by conventional methods. In order to address this issue, a variety of optimization algorithms have been proposed (see Goldstein and Osher 2009 and the reference cited therein). Prior to the recent work of (Micchelli et al 2011), the leading approach was to approximate the ROF functional using an anisotropic Besov norm. In this way, an efficient method known as split Bregman algorithm (Goldstein and Osher 2009) was proposed.

More recently, the approximated solution of the ROF denoising model was computed using a specialized fixed point algorithm (Micchelli et al 2011). By using this new framework, greater gains in convergence speed were demonstrated in the test images. This alternative approach views the ROF functional as a proximity operator of $\alpha {| \cdot | }_{\mathrm{TV}}$, evaluated on the noisy image z (Micchelli et al 2011). While the proximity operator $\alpha {| \cdot | }_{\mathrm{TV}}$ is difficult to compute, it is possible to express it as the composition of a convex function with a linear transformation (Micchelli et al 2011). Thus, equation (1) can be rewritten as:

$$\begin{eqnarray}f(u) & = & \mathrm{argmin}\left\{\displaystyle \frac{1}{2}{\parallel u-z\parallel }_{2}^{2}\right.\\ & & +(\varphi \,\circ \,B)(u):u\in {{\mathbb{R}}}^{d}\},\end{eqnarray} \tag{ 2 }$$

where φ is a convex function on ${{\mathbb{R}}}^{m}$ and B is a m × d matrix defined as:

$$\begin{eqnarray}&&B:= \left|\begin{array}{l}{I}_{N}\otimes D\\ D\otimes {I}_{N}\end{array}\right|,\end{eqnarray} \tag{ 3 }$$

where $\otimes$ denotes the Kronecker product, I_N the N × N identity matrix and D the N × N matrix, defined as:

$$\begin{eqnarray}D:= \left|\begin{array}{llll}0 & & & \\ -1 & 1 & & \\ & . & . & \\ & & . & .\\ & & -1 & 1\end{array}\right|.\end{eqnarray} \tag{ 4 }$$

As it has been shown in section 4 in (Micchelli et al 2011) by choosing $\varphi \,:{{\mathbb{R}}}^{2{N}^{2}}\to {\mathbb{R}}$ the isotropic total variation (ITV) can be defined as a linear combination of the norm ${\parallel \cdot \parallel }_{2}$:

$$\begin{eqnarray}&&\varphi (x):= \alpha \displaystyle \sum _{i=1}^{{N}^{2}}{\parallel \begin{array}{l}{x}_{i}\\ {x}_{{N}^{2}+i}\end{array}\parallel }_{2},\,\,\,\,\,\,x\in {{\mathbb{R}}}^{2{N}^{2}},\end{eqnarray} \tag{ 5 }$$

while the anisotropic total variation (ATV) is defined as the composition of the norm ${\parallel \cdot \parallel }_{1}$:

$$\begin{eqnarray}&&\varphi (x):= \alpha {| x| }_{1},\,\,\,\,\,\,x\in {{\mathbb{R}}}^{2{N}^{2}}.\end{eqnarray} \tag{ 6 }$$

By viewing the TV-norm as a composition map, the explicit form of the proximity map of the l₁- or l₂-norm can be fully exploited.

In ex vivo μCT  a significant amount of data has to be post-processed and denoised prior to segmentation. The algorithmic speed enhancements outlined above represents an important consideration in the selection of an optimal denoising method in practice.

In recent years a non-iterative denoising technique for additive Gaussian noise known as non-local filtering or BM3D proposed by Dabov et al (Dabov et al 2007) has received wide attention in the image processing community. This method is based on an enhanced sparse representation by grouping 2D image fragments into 3D data arrays and applying a collaborative filtering procedure, which consists of 3D transformation, shrinkage of the transform spectrum, and inverse 3D transformation. First, a reference patch is given and by searching over the whole image, similar patches are found and grouped. Secondly for each group of patches found, the 3D wavelet coefficients are computed for denoising by thresholding. Finally the inverse 3D transform is applied and the denoised image is obtained.

This approach considers the following noisy model. For the image domain $X\subset {{\mathbb{Z}}}^{2}$ the noisy image $z:X\to {\mathbb{R}}$ is described by:

$$\begin{eqnarray}&&z(x)=u(x)+n(x),\end{eqnarray} \tag{ 7 }$$

where x is the 2D spatial coordinate that belongs to the image domain X, u is the true image and n is the zero mean Gaussian noise with variance σ². A fixed size N₁ ×  N₁ reference block denoted by ${Z}_{{x}_{\mathrm{ref}}}$ is extracted from the noisy image z within a neighborhood X and a set of blocks similar with ${Z}_{{x}_{\mathrm{ref}}}$ are found after computing the l₂-distance:

$$\begin{eqnarray}&&{d}^{\mathrm{noisy}}({Z}_{{x}_{\mathrm{ref}}},{Z}_{x})=\displaystyle \frac{{\parallel {Z}_{{x}_{\mathrm{ref}}}-{Z}_{x}\parallel }_{2}^{2}}{{({N}_{1}^{\mathrm{ht}})}^{2}},\end{eqnarray} \tag{ 8 }$$

where $\parallel \cdot \parallel$ is the l₂-norm and x is the coordinate of the top-left corner of the block. In this way by using the distance d^noisy a set of image blocks similar S^ht with the reference block ${Z}_{{x}_{\mathrm{ref}}}$ are found:

$$\begin{eqnarray}&&{S}_{{x}_{\mathrm{ref}}}^{\mathrm{ht}}=\{{x}_{\mathrm{ref}}\in X\,:{d}^{\mathrm{noisy}}({Z}_{{x}_{\mathrm{ref}}},{Z}_{x})\leqslant {\tau }_{\mathrm{match}}^{\mathrm{ht}}\},\end{eqnarray} \tag{ 9 }$$

where the ${\tau }_{\mathrm{match}}^{\mathrm{ht}}$ is a parameter defining the maximum distance d for which a selected block is similar with the reference block ${Z}_{{x}_{\mathrm{ref}}}$. This parameter is chosen in an empirical manner. Moreover, the cardinality $| {S}_{{x}_{\mathrm{ref}}}^{\mathrm{ht}}|$ will always be greater than 1 because $d({Z}_{{x}_{\mathrm{ref}}},{Z}_{x})=0$. After obtaining ${S}_{{x}_{\mathrm{ref}}}^{\mathrm{ht}}$ a 3D array of size ${N}_{1}^{\mathrm{ht}}\times {N}_{1}^{\mathrm{ht}}\,\times | {S}_{{x}_{\mathrm{ref}}}^{\mathrm{ht}}|$ denoted with ${{\bf{Z}}}_{{S}_{{x}_{\mathrm{ref}}}^{\mathrm{ht}}}$ is formed by stacking the matched noisy blocks. By using hardthresholding in the 3D transform domain the collaborative filtering of ${{\bf{Z}}}_{{S}_{{x}_{\mathrm{ref}}}^{\mathrm{ht}}}$ is realized.

In case where the variance grows asymptotically and we are dealing with large values for the standard deviation (SD) σ or small values for N₁^ht the BM3D performance will experience a sharp drop. This case is described by the authors as erroneous grouping which can be avoided by using coarse pre-filtering. First, a normalized 2D linear transform is applied on both blocks and secondly the resulting coefficients are hardthresholded:

$$\begin{eqnarray}&&d({Z}_{{x}_{\mathrm{ref}}},{Z}_{x})=\displaystyle \frac{{\parallel {\rm{\Upsilon }}^{\prime} ({T}_{2{\rm{D}}}^{\mathrm{ht}}({Z}_{{x}_{\mathrm{ref}}}))-{\rm{\Upsilon }}^{\prime} ({T}_{2{\rm{D}}}^{\mathrm{ht}}({Z}_{x}))\parallel }_{2}^{2}}{{({N}_{1}^{\mathrm{ht}})}^{2}},\end{eqnarray} \tag{ 10 }$$

where ϒ^' is the hardthresholding operator with threshold λ2Dσ and ${T}_{2{\rm{D}}}^{\mathrm{ht}}$ the normalized 2D linear transform. It has been suggested by the authors to compute ${T}_{2{\rm{D}}}^{\mathrm{ht}}$ using DCT instead of biorthogonal wavelet transform and also to increase N₁^ht from 8 to 12. Using these adjustments the problem can be solved in a certain context, since some parts of the desired signal can be erased and some artefacts introduced when the 2D-DCT transform is applied locally. BM3D has the ability to preserve edges and to avoid undesired smoothing effects by taking advantage of the similarity between the grouped fragments. In this way, a highly sparse representation of the wanted signal is achieved and by applying a shrinkage technique the noise can be well separated.

2.3.1. Shepp Logan phantom

2.3.2. Vascular phantom

In order to assess the performance of the denoising methods a synthetic vascular network was created. In a spherical volume five generations of vessels with varying lengths (14–28 voxels), diameters (1.2–2.4 voxels) and branching angles (3°–77°) were generated using a space-filling algorithm that recursively partitions uniformly distributed points through their center of mass and sprouts daughter vessel segments towards the centroids of the partitioned point clouds. This synthetic network features 3, 4 and 5-vessel junctions. The geometry of this vascular phantom is displayed in figure 2. 720 projections over 180° about the z-axis were calculated and x–y plane cross-sections of 2000 × 2000 pixels were reconstructed using the cone-beam volumetric FDK reconstruction algorithm for a nominal resolution of 1 μm. Noise was added in the sinogram domain for a peak-signal-to-noise ratio (PSNR) of 8.2 dB (see equation (12)) and synthetic images were first reconstructed using FDK, denoised using the selected algorithm in a second step and finally the network topology was extracted using the pipeline shown in figure 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

All μCT-data were obtained on Bruker SkyScan 1172 (Bruker, Kontich, Belgium).

2.4.1. Physical phantom

2.4.2. Ex vivo rat heart

Experimental investigations conformed to the UK Home Office guidance on the Operations of Animals (Scientific Procedures) Act 1986 and were approved by the University of Oxford ethical review board. The heart was excised from a female Sprague-Dawley rat (320 g) after anesthesia and cervical dislocation, and perfused with modified Krebs–Henseleit buffer in Langendorff setup at 37 °C and 80 mmHg. The heart was then cardioplegically arrested with high potassium solution (in (mM): NaCl 125.0; KCl 20; MgCl2 1.0; HEPES 5.0; Glucose 11.0; CaCl2 1.8) and fixed with low osmolality Karnovskys fixative (Solmedia, Shrewsbury, UK). The coronary vasculature was filled by injection with a silicone rubber compound microfil (Flow Tech, Carver, USA). The heart was stored in low osmolarity Karnovkys fixative at 4 °C until imaging time. It was then rinsed in iso-osmotic cacodylate buffer, embedded in 0.5% low-temperature melting agar, and mounted on the platform of the imaging system.

The following parameters were used: x-ray source voltage 55 kVp, source current 181 μA, nominal resolution (pixel size) 4 μm, Al filter, 360° rotation, 0.15° rotation step, frame averaging 5, random movement 10, and camera binning 2 × 2. The size of the heart required the scan be segmented into 4 subscans. The reconstruction parameters were: ring artifact correction 20%, beam hardening correction 50%. The acquisition time was 8 h and the size of the projection images was 3872 × 1336. The 3801 slices were reconstructed in 3 h and 45 min and the cross-section size was 2992 × 3264.

Zoom In Zoom Out Reset image size

The performances of the denoising methods were first evaluated on a 3D Shepp Logan phantom. Since the noise-free data reconstruction was available (figure 5(a)) quantitative measurements of the errors were made using normalized mean square error (NMSE) with the l₂ norm

$$\begin{eqnarray}&&\mathrm{NMSE}=100\,\times \,\displaystyle \frac{{\parallel i-r\parallel }_{2}}{{\parallel i\parallel }_{2}},\end{eqnarray} \tag{ 11 }$$

where i is the ground truth and r the denoised reconstructed image obtained using a denoising algorithm. The quality of the denoised tomograms were evaluated also using the PSNR defined as

$$\begin{eqnarray}&&\mathrm{PSNR}=20\underset{10}{\mathrm{log}}\left(\displaystyle \frac{255}{{| i-r| }_{2}}\right).\end{eqnarray} \tag{ 12 }$$

Zoom In Zoom Out Reset image size

For the experimental data, the overall quality of the reconstructed images was quantitatively assessed in terms of their signal-to-noise (SNR) and contrast-to-noise ratio (CNR), defined on the tomographic slices of the object as

$$\begin{eqnarray}&&\mathrm{SNR}=\displaystyle \frac{| { \mathcal I }| }{\sigma },\end{eqnarray} \tag{ 13 }$$

and

$$\begin{eqnarray}&&\mathrm{CNR}=\displaystyle \frac{| {{ \mathcal I }}_{A}-{{ \mathcal I }}_{B}| }{\sqrt{({\sigma }_{A}^{2}+{\sigma }_{B}^{2})}},\end{eqnarray} \tag{ 14 }$$

where ${{ \mathcal I }}_{A}$ and ${{ \mathcal I }}_{B}$ denote the mean pixel values, and σ_A and σ_B the standard variations within two selected homogeneous regions respectively.

The iterations in the TV denoising algorithms were terminated when the following condition was satisfied

$$\begin{eqnarray}&&\displaystyle \frac{{| {r}_{n}-{r}_{n-1}| }_{2}}{{| {r}_{n}| }_{2}}\leqslant \mathrm{TOL},\end{eqnarray} \tag{ 15 }$$

where r_n is the denoised image at iteration n.

For a simulated object in which the ground truth is known, additional indices can be calculated to evaluate the segmentation accuracy. The similarity between the ground truth A and a binarized denoised data set B quantified by the Jaccard index is (Jaccard 1912)

$$\begin{eqnarray}&&J(A,B)=\displaystyle \frac{| A\cap B| }{| A\cup B| }\end{eqnarray} \tag{ 16 }$$

whereas the Sørensen–Dice coefficient is defined as (Dice 1945)

$$\begin{eqnarray}&&S(A,B)=\displaystyle \frac{2| A\cap B| }{| A| +| B| }.\end{eqnarray} \tag{ 17 }$$

Both measures range between 0 (no similarity) and 1 (identical).

Another important parameter that was measured for all data sets is the noise power spectrum (NPS). The 3D NPS is defined as (Dainty and Shaw 1974):

$$\begin{eqnarray}&&\mathrm{NPS}(u,v,w)={ \mathcal F }\{A({\tau }_{x},{\tau }_{y},{\tau }_{z})\},\end{eqnarray} \tag{ 18 }$$

where ${ \mathcal F }$ is the 3D Fourier transform operator u, v and w are spatial frequency variables corresponding to spatial variables x, y and z. $A({\tau }_{x},{\tau }_{y},{\tau }_{z})$ defines the system autocovariance where τ_x, τ_y and τ_z represent the distance between voxels in x, y and z directions. In our analysis a 3D approach was used in order to provide an accurate NPS over each 3D data sets. By averaging the 2D slices of the volume and by using several VOIs and a 3D FFT the correlation between the images is explored. A similar technique was proposed using the CT ACR accreditation phantom, where two similar scans were available (Friedman et al 2013). It is known that the 2D NPS is strongly dependent on the VOI position. For this reason, 4 VOIs were chosen for each data set and for each VOI the 3D NPS was calculated. The 4 VOIs overlapped in all denoised volumes by 30%. The final NPS was obtained by averaging the results obtained for each VOI in the radial direction.

In all the studied methods, the regularization parameters were chosen carefully by trial-and-error in order to obtain the best decrease and to avoid divergence. In this sense, large and small values of the parameters leading to poor reconstruction results were identified first. Secondly, the optimal values of the parameters were then gradually refined by trial-and-error with a decreasing interval. For each parameter different VOIs were chosen in the denoised image and different metrics (SNR, CNR and SD) are measured. A plot profile between the original image and the resulting image after applying a denoising algorithm for these VOIs was visualized in order to assess the behavior of the algorithm for a given regularization parameter.

In this work a 3D vascular network reconstruction algorithm extending on the methodology outlined in (Goyal et al 2013) was used. We present an enhancement of this work by improving the vascular detection and extraction pipeline in numerous stages. The pipeline of this new extraction model is shown in figure 3. This modified pipeline incorporated the selected denoising algorithm as a first pre-processing step and several other changes to improve code efficiency and reconstruction accuracy.

Firstly, the multiscale Frangi vesselness descriptor (Frangi et al 1998) was expanded by complementing it with the addition of the Sato vesselness measure (Sato et al 1998). The result was an increased response and continuity at vessel junctions in comparison to the vesselness filters originally used. In addition, the sphere-fitting vessel radius estimation (Hildebrand and Rüegsegger 1997) was replaced in our new implementation by Rayburst sampling (Rodriguez et al 2006). The Rayburst sampling algorithm, which, when combined with a 2D/3D core that follows the vessel centerline, could be used for diameter estimation by locally adjusting the ray exit thresholds at each step thus improving accuracy.

To take the partial voluming of the vessel lumen into account, we calculate the vessel diameter for each casted ray by evaluating the full width at half maximum (FWHM) of the ray profile instead of simply taking a hard threshold. The diameter at a particular centerline position was then computed from the 3D sampling core by using the median lower band diameter method (Wearne et al 2005). Furthermore, intensity spikes in the ray profile that would otherwise corrupt the FWHM algorithm were filtered out by morphologically opening the signal with a structuring element of size 25% of the estimated ray length.

These steps resulted in a relatively well-defined skeleton centerline with a radius at each point. However, any inaccuracies inherent to medial axis thinning skeletonization still remained. To remove these, we performed a number of post-processing steps based on the information available. Firstly, we removed spurious branches by targeting junctions whose mother and daughter vessels were both sufficiently large (>5 voxels in diameter) and their ratio exceeded a pre-specified threshold. We then continued by removing loops within the skeleton by iteratively computing Dijkstra's shortest path algorithm (Dijkstra 1959) starting from a pre-defined inlet vessel until no more segments were pruned. For large segments (>5 voxels in diameter) we also smoothed out the outlier nodes in the centerline by recomputing the shortest path from the start to the end node. The costing between nodes was defined as a function of the distance between them. This resulted in nodes furthest away from the true medial axis of the vessel being pruned out.

Finally, a number of performance enhancements were made through the use of memory maps, GPGPUs, custom multi-PC parallelization and vectorized functions. We have successfully processed volumes of up to 4000³ voxels, distributed to standard desktop machines, limited only by the number of machines available. This pipeline is fully automatic, but manual adjustments can be made after the thresholding and during the network pruning steps when required.

3.1.1. Shepp Logan phantom

3.1.2. Vascular network

The synthetic network described in section 2.3.2 was used in order to assess the efficacy of the denoising algorithm. The 3D binary network geometry of this synthetic phantom is shown in figure 2. Initially the network was reconstructed using the FDK algorithm with Poisson noise added in the sinogram domain yielding PSNR = 8.2 dB. In figure 6(a) a cross-section corrupted by Poisson noise is displayed together with the corresponding cross-sections denoised via Gaussian smoothing (window size 10 and the SD of the Gaussian window 5), median filter (10-by-10 neighborhood matrix parameter), wavelet with soft thresholding (using wavelet Daubechies 6, level of decomposition 9 and threshold 120), BM3D (σ = 100) and ITV approaches (α = 100 and TOL = 10⁻³). The 3D vascular extraction algorithm described in section 2.6 was performed for seven cases: noise-free data, FDK Poisson noise, Gaussian smoothing, median filter, wavelet with soft thresholding, BM3D and ITV reconstructions. The 3D reconstructions were obtained by using the same parameters in the vascular extraction algorithm for all the cases.

A comparison between the network topologies was performed by quantifying the difference between the 3D noise-free data and the 3D networks extracted with the vascular reconstruction pipeline. First, in figure 8(a) the ideal synthetic vascular reconstruction (red) is shown overlapped with the Poisson noise FDK reconstruction (white). In the same fashion, the overlapped ideal network and the reconstructed denoised networks via Gaussian smoothing, median filter, wavelet with soft thresholding, BM3D and ITV methods are shown in figures 8(b)–(f) respectively. The red structures displayed in figure 8 showing overlapped reconstructions represent the missing structures in the reconstructions. By applying the denoising methods less information is lost and the topology of the network is preserved. The initial FDK NMSE of 47.6% was decreased by all the five denoising approaches, and the best result was obtained in the ITV reconstruction (NMSE = 31.3%) (table 5). The reduction in the noise level can be also observed in figure 6 where the central cross-section of the synthetic vascular phantom is displayed for the different techniques. Interestingly, while the median filter (figure 6(c)) was almost as effective at removing background noise as the ITV technique, evidences of over-filtering is apparent in the form of severely diminished contrast in the smaller vessels.

Zoom In Zoom Out Reset image size

Table 5.  Normalized mean square error (NMSE), Jaccard index (JI) and Sørensen–Dice index (SDI) for the 3D binary vascular volumes in the synthetic phantom.

Method	NMSE (%)	JI	SDI
FDK reconstruction	47.6	0.77	0.87
Gaussian filter	43.8	0.82	0.89
Median filter	41.8	0.83	0.9
Wavelet	40.7	0.84	0.91
BM3D	38.8	0.85	0.92
ITV	31.3	0.91	0.95

The ITV technique produced a Jaccard index of 0.91 and Sørensen–Dice index of 0.95, compared to the FDK baseline values of 0.77 and 0.87 respectively. According to these metrics the ITV denoising method gives the smallest segmentation errors compared to the other four denoising techniques (see table 5). In table 6 the correlation coefficients obtained between the residuals images corresponding to each denoising method and the first derivative of the noise-free data are displayed. The best value was achieved using the ITV method (0.3) followed by BM3D, wavelet, median and Gaussian filtering.

The 3D NPS estimated was calculated for the following VOIs: 50 × 50 × 50, 150 × 150 × 150, 250 × 250 × 250 and 375 × 375 × 375 voxels. The mean values of the axial NPS are displayed in table 4 (line 2). Very good results were obtained with the ITV method. The representative planes from the 3D NPS for w = 0 are shown in figure 7(b). The global noise is significantly reduced by BM3D but better results are reconstructed using ITV.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.2.1. Physical phantom

The baseline FDK reconstruction of the central transverse slice of the phantom with the acquisition and reconstruction parameters described in section 2.4 is displayed in figure 9(a). The corresponding slice denoised using Gaussian-smoothing (window size 15 and the SD of the Gaussian window 9), median filter (15-by-15 neighborhood matrix parameter), wavelet with soft thresholding (using Wavelet Daubechies 6, level of decomposition 9 and threshold 70), BM3D (σ = 30) and the ITV (α = 100 and TOL = 2 × 10⁻³) algorithms are shown in figures 9(b)–(f) respectively. Two homogeneous regions of the phantom have been chosen in order to evaluate the performance of the algorithm. These two zones are displayed in figure 9(f) as a blue square (region A) and a red circle (region B). The SD, SNR and CNR are reported in table 7 for all the denoising methods. Compared to the baseline FDK reconstruction, the application of the ITV algorithm increased the SNR in region A by a factor of 6.8 and by a factor of 12.3 in region B, respectively. Thus, this resulted in a 8.2-fold higher CNR.

Zoom In Zoom Out Reset image size

Table 6.  Correlation coefficients obtained between the residuals images corresponding to each denoising method and the first derivative of the free-noise data for the 3D binary vascular volumes in the synthetic phantom.

FDK	Gaussian smoothing	Median filtering	Wavelet	BM3D	ITV
0.013	0.042	0.044	0.035	0.12	0.3

Table 7.  Numerical results (standard deviation (SD), signal-to-noise ratio (SNR) and contrast-to-noise ratio (CNR)) obtained for the physical phantom.

Method	Region	SD	SNR (dB)	CNR (dB)
FDK	A	34.7	2	0.4
	B	32.5	1.6
Gaussian	A	10.9	6.1	1
	B	9.8	5.2
Median	A	8	8.5	1.6
	B	7.5	6.7
Wavelet	A	7.9	8.6	1.6
	B	6.7	8
BM3D	A	8.1	8.3	1.7
	B	7.3	7
ITV	A	5.2	13.4	3
	B	2.6	19.8

In the MR physical phantom case, the mean value in the axial direction for the NPS along w = 0 was established using the following VOIs: 500 × 500 × 150, 750 × 750 × 175, 1000 × 1000 × 125 and 2000 × 2000 × 200 voxels. The results are presented in table 4 (line 3). In figure 7(c) the 3D NPS representative planes are displayed for all denoising methods. ITV achieved the lowest mean NPS value outperforming the other methods.

3.2.2. Ex vivo rat heart

In figure 10(a) the FDK reconstruction of the central transverse and sagittal slices are displayed together with the Gaussian smoothed slice in figure 10(b). The central cross-section denoised by the median filter is shown in figure 10(c) and the corresponding wavelet with soft thresholding, BM3D, ITV results in figures 10(d)–(f) respectively. The Gaussian smoothing parameters were set to 10 for window size and 4 for the SD of the Gaussian window. In the median filter case a 10-by-10 neighborhood matrix parameter were found to be optimal. The wavelet parameters were level of decomposition 9 and threshold 40, whereas BM3D parameter was σ = 150 and ITV parameters α = 1000 and TOL = 5 ×  10⁻⁴ have been found to be optimal.

Zoom In Zoom Out Reset image size

On this set of central cross-sections the SNR and the CNR were measured in a homogeneous area of the heart. As highlighted with red rectangles in figure 10(a). The results obtained for the initial FDK reconstruction and the denoising methods are reported in table 8. The ITV method performed better than the other denoising algorithms by increasing the SNR by a factor of 7.5 and the CNR by a factor 10.1 compared to the initial FDK reconstruction. The improvement of the denoising methods can also be seen in figures 10(g) and (h) where the ITV method yielded a smoother profile than the other methods. Profiles corresponding to a transverse and sagittal section of the heart were plotted in figures 10(g) and (h) respectively.

Table 8.  Signal-to-noise ratio (SNR) and contrast-to-noise ratio (CNR), obtained in a homogeneous area of the ex vivo rat heart using different denoising methods.

Method	SNR (dB)	CNR (dB)
FDK	2.9	1.6
Gaussian	6.1	4.2
Median	9.1	6.4
Wavelet	16.5	10.6
BM3D	14.7	8.9
ITV	21.5	16.2

The 4 VOIs used to compute the mean NPS value in the axial direction for the ex vivo rat heart were: 50 × 50 × 50, 75 × 75 × 75, 125 × 125 × 125 and 150 × 150 × 150 voxels. The results are displayed in table 4 (line 4). The lowest mean NPS value was obtained for the images denoised with ITV, but comparable results were achieved with the BM3D and the wavelet techniques. The representative planes from the 3D NPS are displayed in figure 7(d). The original FDK noise was reduced by all denoising methods, but the best NPS estimate was achieved with ITV.

Using the algorithm described in section 2.6 a 3D vascular network was reconstructed from a 300-slice thick slab in the μCT  data corresponding to the apex of the heart.

The first processing step of our 3D reconstruction algorithm is the vesselness filter, followed by binary mask generation. For the small vessels (i.e. <25 μm diameter) in this ex vivo μCT  dataset it can be seen that the vesselness filter does not distinguish vessel-like structures from the noise structures in the image stacks for the FDK, Gaussian smoothing or median filter algorithms. This had two major consequences for the final 3D reconstruction of the coronary vascular network: (i) spurious microstructures associated with noise were detected (see below for further comments) and (ii) the computational time was increased by a factor of 20. On the other hand, the masks obtained after the vesselness filtering step using wavelet or BM3D methods were significantly improved, on par the ITV images in terms of noise level. Moreover, the wavelet CNR and SNR were higher than the BM3D values but the small vessels were not preserved. In the BM3D thresholded image the big vessels were frequently fused or discontinuous. The improvement in 3D model reconstruction after incorporating the denoising models varied significantly amongst the different filters (see top row of figure 11). Whereas the wavelet method only allowed extraction of the largest vessels, BM3D filter led to numerous spurious interconnections between segments. Many of these arose from the noisy mask generation which contained several cusp points, that are detected as segments by the skeletonization algorithm and can easily be seen on comparison with volume rendering of the image data. It is also noteworthy that 3D reconstruction was not possible for baseline FDK, Gaussian-smoothed or median filtered images due to the high level of remaining noise. In contrast, the ITV reconstruction captured the bifurcating topologies of both the small and the large vessels of the network without introducing spurious connections.

Zoom In Zoom Out Reset image size

In previous studies, the quality of the extracted network could be evaluated by acquiring a second scan of a sub-region at a higher imaging resolution (e.g. 25 μm versus 4 μm in Lee et al 2007). However, given this dataset was already close to the resolution limit, the gain in feature delineation expected from a second scan would have been minimal. As the lack of ground truth makes this evaluation difficult, we further provide a qualitative comparison between the denoised data and the binarized segmentations resulting from applying the vesselness filters. Maximum intensity projection (MIP) and raw slice data are both used, as the former makes the vascular structures more apparent while the latter gives a truer representation of the denoising filter performances.

Shown in the second row of figure 11 are representative slices from the various denoising filter outputs. The corresponding MIPs of the bottom half of the ROI are shown in the third row. The Otsu-thresholded outputs of the vesselness filter are shown in the final row, noting that the outcome of the thresholding is strongly dependent on the degree of overlap between the intensity distributions of the vessels and noise clusters. The combined results strongly suggested that the wavelet filter has a tendency to remove the finest vessel features, but maintains the intensity level of the noise relatively high such that it carried through to the final vessel mask generation. The segments resulting from the noise were however disconnected from the main network as the smaller vessels segments were missing and were subsequently removed by the model extraction code. In contrast, BM3D method performed best in terms of preserving the vessel boundaries, but at the cost of relatively high noise levels which persisted through the vascular processing pipeline, culminating in the spurious connections observed. The ITV filter was less effective at preserving the sharp boundaries but was able to suppress the noise to such a level that almost none was detected by the vesselness and thresholding steps, allowing a greater extent of the network to be recovered. Therefore in our tests with the real data, ITV filter was able to achieve the best compromise between separating the intensity distributions of noise from vessels while preserving vessel contrast.

A key limiting factor in the current study is a lack of a ground truth dataset in the real data evaluation. We have partially addressed this issue by employing a number of different phantoms. However, in a real data many more complications may be present such as inhomogeneous background intensity distribution and contrast extravasation as well as complex noise characteristics which differ from the simulated data. Ways of overcoming these difficulties may include using a synchrotron μCT  to obtain a higher quality reference data, or to employ 3D printing of a known network structure, which can then be imaged. Neither option was possible at the time of this work, as the resolution of the 3D printing techniques has not reached a sufficient level for microvascular network reproduction and the synchrotron facilities were not readily accessible. The future introduction of these resources will allow a quantification of the qualitative enhancements demonstrated in the current work.

Further, we address the potential criticisms regarding the selection of the denoising filters here, by clarifying the rationale behind the selection criteria. The broader motivation of this work has been to identify an optimal denoising methodology that can be widely taken up by research communities. In many of the laboratories that utilize μCT, the primary focus is on applying these scanners to answer a scientific question rather than to conduct in-depth image processing analyzes. While more sophisticated schemes may generate superior outcomes, their complexity will be a barrier for applications in many general research settings. Thus we have prioritized filters that have practical ease of use, with limited number of parameters (all filters in this work have a single parameter). In addition, the implementations of these filters can be readily found in standard programming languages and web resources. To assist with this goal, our Matlab implementation of the TV filters can be found at http://havmgroup.com/tools/. Data and algorithms will be made available upon request.