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

Micro-computed tomography systems are widely used for high-resolution, non-destructive analysis of internal microvascular networks. When the scale of the targeted vessel approaches the imaging resolution limit, the level of noise becomes a limiting factor for accurate reconstruction. Denoising algorithms provided by vendors are often suboptimal for enhancing SNR of fine (vessel) features. Furthermore, the performance of existing methods has not been systematically analyzed in the context of final network reconstruction and graph model extraction. This work evaluates several standard and state-of-the-art noise reduction techniques using both in silico and physical phantoms, and ex vivo rat coronary data for their ability to improve vascular network analysis. We compared five noise reduction approaches, including vendor-supplied (Gaussian smoothing), conventional (median filter) and advanced (i.e. wavelet filter with soft thresholding, block-matching collaborative filtering (BM3D), and isotropic and anisotropic total variation denoising) techniques. The latter two methods were chosen for their reported ability to preserve fine details, a prerequisite for a successful microvascular extraction. The full evaluation pipeline included the reconstruction from projection images, denoising, vascular segmentation and graph model extraction to be performed on all simulated and real image data sets. SNR, CNR and 3D NPS were quantified from denoised images, and where the ground truth was known, Sørensen–Dice coefficients, Jaccard index metrics were calculated as measures of segmentation error. The performance of the image denoising algorithms where the ground-truth was available has been assessed by computing the correlation coefficients between the residual images (obtained between the noise-free data and the denoised data) and the first derivative of the noise-free data were computed. Overall, simpler denoising techniques including the median and wavelet filters and the vendor-supplied implementations have been found to perform inadequately for segmentation of fine vessel features, particularly on real images. BM3D technique performed well in most of our tests, however isotropic total variation (ITV) was the optimal choice for noise reduction and feature preservation in real data as shown by the extracted network models. Globally, ITV increased the SNR from 10.2 to 31.7 dB in a Shepp Logan phantom, doubled SNR and CNR values in a scanned physical phantom compared with BM3D, enabled the smallest vessels to be fully recovered in an in silicon phantom and achieved a near-ideal outcome in the rat coronary data.


Introduction
X-ray micro-computed tomography (μCT) has emerged as a powerful tool to study biological samples due to its non-destructive nature and the ability to investigate three-dimensional internal structures of specimen with high resolution (Jorgensen et al 1998, Holdsworth and Thornton 2002, Schambach et al 2010, Ritman 2011. Contrast-enhanced imaging of the vascular networks is a major application area of μCT, as the detailed extent of the vascular topology and morphology achievable with μCT are of key interest for studying function, adaptation, energetics, disease and remodeling of the organ. To date, images of small animal heart (Lee et al 2007), kidney (Nordsletten et al 2006), liver (Wan et al 2000), tumours (Badea et al 2006 ) and brain (Dorr et al 2007) vascular structure-function relationships have been investigated using μCT. In order to allow efficient and reproducible processing of such data, we have previously developed automated processing pipelines for reconstruction of full-organ vascular networks from 3D images (Lee et al 2007, Goyal et al 2013. They were successfully applied to synchrotron μCT, clinical CT, MR and cryomicrotome images, and accurately reconstructed the vascular networks in the presence of moderate to low noise levels (figure 1).
For this type of application, where a resolution in the order of micrometers (μm) is required, commercial μCT systems are a ubiquitous alternative to the synchrotron systems (Brunke et al 2008), although they are limited by less favorable characteristics in terms of gross noise (Tang et al 2011). Most commercial μCT scanners are cone-beam micro-computed tomography (CBμCT) systems. By rotating the studied object over 180°or 360°in discrete steps, a set of transmission x-ray images are acquired in each angular position. When the acquisition process is complete, orthogonal slices through the object are obtained by applying a reconstruction algorithm such as the Feldkamp (FDK) method (Feldkamp et al 1984). The performance and the quality of a CB μCT system are strongly influenced by multiple factors, such as photon scatter, scanner geometry, scanner misalignment, x-ray spectrum, detector pixel size, voxel size, number of projections, reconstruction algorithm, absorption and detection of x-ray photons, electronic noise and quantum noise from random generation (Tu et al 2006, Zhang and Ning 2008, Baek and Pelc 2011. Another important aspect of an x-ray μCT imaging system is the blurred effect of the edges created when the scale of a vessel is reached. This effect is known as partial volume effect (PVE) and can affect the surrounding voxels of the edges of the small vessels and create false connections between two adjacent structures. Even if the PVE is decreased by using a higher spatial resolution, some inverse filtering techniques have been developed in PET and SPECT in order to correct these artifacts (Erlandsson et al 2012).
Recently, a μCT technique to reduce this blurring effect was proposed in (Kline et al 2014) by placing the specimen (a rat liver) very close to the detector array. Is well known that PVEs depend on several factors: type of scanner (spatial resolution, x-ray source focal spot and detector), the size of the specimen, the range of the attenuation coefficients, motion, and other temporal effects. In our study, by applying post-reconstruction techniques like the wavelet regularization, BM3D and ITV are known to preserve fine details, the blurring effect associated with PVE and with the reconstruction algorithm should be decreased.
Moreover, the image quality of the x-ray system can be measured quantitatively via a modulation transfer function (MTF). MTF is a metric describing the ability of the imaging system how accurate the smallest details of a studied object are reproduced (Judy 1976). MTF is strongly dependent on the x-ray imaging system. Different techniques to measure MTF were developed for clinical CT system using standardize accreditation phantoms (Silverman et al 2009, Richard et al 2012, Friedman et al 2013 and for μCT system (Rueckel et al 2014). In μCT vascular images the values in the opacified vessels are decreased and the small vessels diameter is expanded because the MTF drops progressively at high frequency (Kline and Ritman 2012).
While the majority of the undesired effects caused by these factors can be controlled for, random noise remains an unresolved key issue that impacts negatively on the quality of the μCT images. Thus, the measurement of various parameters such as volume and texture which is performed directly on the Figure 1. Examples of the 3D extraction of full-organ vascular networks using the proposed automated processing pipelines described in section 2.6 from cryomicrotome images for: (a) dog heart, (b) pig heart (Goyal et al 2013), (c) pig heart and from μCT images for (a) rabbit kidney (Alastruey et al 2014). reconstructed images can be compromised if a significant noise is present. Furthermore, segmentation and more advanced image analysis such as the automated extraction of micro-structures can be severely hampered. These applications necessitate the use of denoising algorithms as a pre-processing step. Two different denoising schemes are possible in a reconstruction work flow: • Pre-reconstruction denoising is applied directly on the projection data and explores the statistical properties of the noise available in the sinogramdomain (Wang et al 2008). Different approaches were proposed for (μ)CT, including iterative numerical algorithms (Elbakri and Fessler 2003), the penalized-likelihood approach (Rivière 2005), nonlinear (Zhong et al 2004) and edge-preserving noise filters (Wang et al 2005). It was demonstrated that both directed noise and streak artifacts are reduced by these methods. However, if the noise variance is constant over the projection set, these methods do not perform well and suffer from an additional loss of spatial resolution in the reconstructed images (Matrecano et al 2010).
• In order to explore the image characteristics postreconstruction denoising is preferred (Rust et al 2002). In this domain, a statistical noise model is no longer available, although structural properties of the objects may be exploited instead. Experimental studies showed a band-pass behavior in this domain, i.e. the noise is non-white and nonstationary. These techniques can be directly applied on the reconstructed images without having access to the projection data. For this reason the postreconstruction approach may be preferred for general research application settings in which the vendor software cannot be modified.
In commercial μCT systems, often the range of implemented denoising methods tends to be limited in comparison to clinical CT systems, e.g. to simple linear smoothing filters (symmetrical, asymmetrical or Gaussian Bruker 2011). Since the power of the noise is concentrated in high and medium frequencies (Boedeker et al 2007), these types of filters cause blurring of fine structures and loss of sharp edges. In order to preserve details, more advanced techniques such as wavelet filter with soft thresholding (Donoho 1995), shrinkage schemes (Chambolle et al 1998, Beck andTeboulle 2009b) and nonlinear methods have been proposed. Various classes of algorithms belong to the nonlinear category, including non-local means denoising (Buades et al 2005), total variation (TV) denoising (Rudin et al 1992), or simpler ones such as the median and bilateral filters which perform well against impulsive noise (Tomasi and Manduchi 1998). The non-local means based BM3D filter (Dabov et al 2007) has been previously evaluated on μCT data of bone microstructure, which demonstrated its superiority over wavelet shrinkage schemes (Matrecano et al 2010). However, to date there is a lack of a systematic and quantitative comparison of various denoising schemes on vascular images, and their influence on automated vessel extraction.
Recently, an alternative formulation of the TV framework was proposed that led to an efficient numerical solution, overcoming the long standing issue of slow convergence in TV-based denoising methods (Micchelli et al 2011). The novelty of this approach was in treating the denoised image as the proximity operator of the TV calculated from a given noisy image. The convergence of the algorithm was accelerated by a specialized fixed-point iteration scheme. While these TV models were only tested on standard 2D images with white Gaussian noise profile, the geometric sensitivity of the TV-based methods to the underlying features, paired with the typical images targeted in this work (comprising contrast-enhanced vessel lumen volume against a relatively uniform background), suggests a potential benefit to be gained in the hitherto unreported application to μCT data.
Despite the advances in denoising methods as outlined above, their efficacy in enhancing vascular model extraction remains unknown to date. In this study, we conduct a comparative systematic evaluation of major classes of post-reconstruction denoising methods for the ultimate goal of automated vascular image analysis from contrast-enhanced μCT datasets. A representative selection of five techniques are assessed, including the linear Gaussian filter, median filter, wavelet denoising with soft thresholding (Beck and Teboulle 2009b), BM3D (Dabov et al 2007) and the TV denoising (Micchelli et al 2011). Simulated and experimental image data are employed to quantify SNR, CNR and mean error characteristics. The denoised images are then processed by an automated vascular processing pipeline to assess the impact of the denoising methods on the reconstructed vascular networks.

Methods and materials
We begin by briefly reviewing the recently proposed TV denoising models and BM3D, before providing details on simulations, experiments and the 3D vascular network reconstruction methodology.

Total variation μCT image denoising
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: where Î  u d is the desired true solution computed at each iteration, z is the observed or noisy data, α>0 is the regularization parameter and =  | | u u TV 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 a |·| TV , evaluated on the noisy image z (Micchelli et al 2011). While the proximity operator a |·| 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: where j is a convex function on  m and B is a m×d matrix defined as: where Ä denotes the Kronecker product, I N the N×N identity matrix and D the N×N matrix, defined as: while the anisotropic total variation (ATV) is defined as the composition of the norm · 1   : By viewing the TV-norm as a composition map, the explicit form of the proximity map of the l 1 -or l 2norm 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.
2.2. Block-matching and 3D filtering μCT image denoising 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 2 the noisy image   z X : 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 σ 2 . A fixed size N 1 × N 1 reference block denoted by Z x ref is extracted from the noisy image z within a neighborhood X and a set of blocks similar with Z x ref are found after computing the l 2 -distance: where ·   is the l 2 -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 ref are found: is a parameter defining the maximum distance d for which a selected block is similar with the reference block Z x ref . This parameter is chosen in an empirical manner. Moreover, the cardinality | | S x ht ref will always be greater than 1 because is formed by stacking the matched noisy blocks. By using hardthresholding in the 3D transform domain the collaborative filtering of Z S x ref 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 1 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: where ϒ ′ is the hardthresholding operator with threshold λ2Dσ and T 2D ht the normalized 2D linear transform. It has been suggested by the authors to compute T 2D ht using DCT instead of biorthogonal wavelet transform and also to increase N 1 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.

Simulations 2.3.1. Shepp Logan phantom
A 3D Shepp Logan phantom consisting of a series of ellipsoids was first simulated in order to evaluate the denoising approaches (Kak and Slaney 1989). 600 projections over 180°were calculated about the longitudinal axis in a cone-beam geometry at 60 kVp. The nominal resolution used for this simulation was 1 μm. After applyingthe cone-beam volumetric FDK reconstruction algorithm, transverse cross-sections of 512×512 pixels were obtained. The reconstructed central slice of the phantom was used for comparison as the 'ground-truth' reconstruction. In addition, different levels of Poisson and Gaussian noise were added in the sinogram domain and reconstructed using the FDK algorithm. For each level of noise and for each filter, the optimal parameters were determined empirically and remained fixed for denoising of successive slices. The simulations were implemented in Matlab R2013b (Mathworks, Natick, MA, US) running under Windows 7 on a laptop with a 3 GHz Intel Core i7 CPU and 16 GB RAM.

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.

Physical phantom
A plastic phantom of diameter 36 mm containing a various internal structures surrounded by copper sulphate solution was scanned over 180°with 0.15°r otation step at 80 kVp (source current 124 μA) at a nominal resolution (pixel size) of 2 μm, using an Al +Cu filter. The frame averaging parameter was set to 6 and the random movement parameter to 10. While setting these two parameters to large values resulted in long acquisition times (i.e. 14 h), this compromise was made in order to acquire the best image quality. The projection size was 2672×7744. The reconstruction parameters were: beam hardening correction 43% and ring artifact 20%.The beam hardening correction, ring artefact reduction and random movement parameter were set within the NRecon reconstruction software (Bruker 2011). The cross-section size was 7744×7744 and 411 slices were reconstructed in 47 min using proprietary NRecon software (Bruker, Kontich, Belgium).

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°r otation 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.

Assessment of reconstructed image accuracy
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 2 norm 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 For the experimental data, the overall quality of the reconstructed images was quantitatively assessed in terms of their signal-to-noise (SNR) and contrastto-noise ratio (CNR), defined on the tomographic slices of the object as where  A and  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 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 whereas the Sørensen-Dice coefficient is defined as (Dice 1945) 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): x y z where  is the 3D Fourier transform operator u, v and w are spatial frequency variables corresponding to spatial variables x, y and z.
x y 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-anderror 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.

3D vascular model extraction pipeline
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 3 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.

Results and discussion
3.1. Simulated data 3.1.1. Shepp Logan phantom In order to decrease the normalized mean square error and to improve the signal to noise ratio of the reconstructions the central cross-section of the simulated phantom was denoised using three more advance techniques listed previously (wavelet, BM3D and TV). Moreover, only one result for the TV method is given since ITV and ATV methods yielded very similar results.
An important aspect of the ITV/ATV algorithm is the stopping criteria (equation (15)). For this reason several TOL values were tested for the simulated data where the ground-truth is known and the convergence of the algorithm to the optimal solution could be validated. By stopping the ITV/ATV algorithms when the TOL value of 10 −4 is reached we ensure that the algorithm converged to the best solution. Is possible to set this value for lower values but the algorithm will stagnate and the convergence time will increase. Moreover, if the TOL value is larger then the algorithm will be stop before finding the optimal solution.
In order to assess the performance of the ITV algorithm versus the ATV algorithm in tables 1 and 2 different metrics for the Shepp Logan phantom corrupted with Poisson noise and Gaussian noise respectively are presented. Moreover, using the same stopping criteria (i.e. same TOL) it can be observed that the performance of the two mentioned algorithms is strongly dependent on the regularization parameter. By carefully choosing this parameter the convergence of the algorithm is guaranteed. On the other hand, the optimal solution is found after several iterations increasing the processing time. This aspect can be critical when large data sets are involved. For the Shepp Logan phantom the optimal ITV parameter was set to α=35.8 (table 1).
Using Matlab implementations, the wavelet Daubechies 6 package was found to perform best compared to other wavelet packages for the level of decomposition 4 and the soft threshold 32. In all our data sets an optimal BM3D parameterσ (the value of the noise) was found measured via SD and after refine by visualizing the results. In this case, the optimal BM3D parameter was found to be σ=20.
The residual images between the noise-free Shepp Logan central cross section and the denoised images are displayed in figure 4. The performance improvement was quantitatively assessed by computing the correlation between the each residual image and the first derivative of the noise-free Shepp Logan. The highest correlation coefficient was found for the ITV method (0.57), followed by the Gaussian smoothing (0.33), wavelet filter (0.3), median filter (0.24) and finally by the BM3D (0.11) (table 3). This finding suggests that the ITV method preserve the edges of the phantom better than the other denoising algorithms.
In all cases, superior NMSE and PSNR were obtained compared with the wavelet and BM3D techniques. The lowest value for the NMSE was 3.3% given by the ITV method, with 9.6% for wavelet with soft thresholding and 5.9% for BM3D. In the same fashion the best PSNR value was 29.6 dB obtained by ITV, and Table 1. Normalized mean square error (NMSE), peak-signal-to-noise ratio (PSNR) and execution time, obtained after applying ITV and ATV algorithms for the Shepp Logan phantom with Poisson noise (PSNR=10.2 dB) using the same stopping criteria (TOL). The convergence of the ITV and ATV algorithms is strongly dependent of the regularization parameter α. In order to converge to the optimal solution the convergence time is increased. 20.3 dB and 24.6 dB for wavelet and BM3D methods respectively. After applying these six denoising methods (including ATV) the closest reconstruction to the 'ground truth' image ( figure 5(a)) in the Shepp Logan case was obtained using ITV technique, as shown in figure 4(f). This is confirmed in the profile plots ( figure 5(b)), which depict cross-sections through a homogeneous region of the phantom. It can be observed that the borders of the analyzed region are well preserved, and only small deviations from the ground truth are present.
Moreover, for the ITV/ATV methods very similar results were obtained when Gaussian noise was added (data not shown). In this case the simulation time was approximately the same as for Poisson noise, but the value of the regularization parameter α was smaller. The best results were obtained using the ITV algorithm as summarized in table 2. The PSNR of the initial Shepp Logan slice was 10.2 dB and after 15.6 s a denoised image with PSNR=31.7 dB and NMSE=2.6% was obtained.
All presented methods reduced the global error level of the reconstructed images compared to the noisy baseline case, but the best results were obtained using the ITV/ATV approach. The values of the NMSE and the PSNR (tables 1 and 2) were significantly improved by the TV denoising technique. ITV and the ATV models gave similar results, but the ITV method converged faster (by around 10%) for a given level of performance, and was able to converge to a better solution. Thus, the ITV algorithm was subsequently used in the following sections.
Another important metrics computed for all the denoised images was the 3D NPS (equation (18)). The NPS measurements in the radial (axial) direction were obtained by choosing four different VOIs in each denoised volume. Depending on the original size of the analyzed volume, the size of the VOIs (voxels) were established. For the Shepp Logan phantom, the following VOIs size were used: 50×50×50, 75×75×75, 125×125×125 and 150×150× 150 (voxels). The 4 VOIs overlapped in all the denoised volumes by 30%. Finally, the NPS axial (w=0) mean value for each VOIs was calculated, and by averaging the results of the 4 VOIs the values displayed in table 4 were obtained. For the Shepp Logan phantom the lowest mean NPS value was obtained for the images denoised with BM3D (line 1), but very similar results were achieved with the ITV method. By visualizing the representative planes from the 3D NPS displayed in figure 7(a) it can be seen that the original noise form the FDK reconstruction was reduced by all the denoising method but the best NPS estimated was achieved with BM3D.

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 Table 2. Normalized mean square error (NMSE), peak-signal-to-noise ratio (PSNR) and execution time, obtained after applying ITV and ATV algorithms for the Shepp Logan phantom with Gaussian noise (PSNR=10.2 dB) using the same stopping criteria (TOL). The convergence of the ITV and ATV algorithms is strongly dependent of the regularization parameter α. In order to converge to the optimal solution the convergence time is increased. (σ=100) and ITV approaches (α=100 and TOL=10 −3 ). 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, Table 4. Mean value of the 3D NPS in the axial direction computed by averaging four different VOIs (with an overlap of 30%).  . Representative plane in the radial (axial) direction along w=0 from the 3D NPS measured for (a) Shepp Logan phantom, (b) synthetic 3D vascular network (in the ITV case the maximum of gray level for the 2 NPS plane was set to 10 −9 compared with the rest of the methods where it was 10 −7 ), (c) physical MR phantom and (d) apex of ex vivo rat heart. For each data set six representative 2D NPS planes corresponding to: FDK reconstruction, Gaussian smoothing, median filtering, wavelet filtering, BM3D and ITV denoising techniques were illustrated. It can be see that the initial noise is reduce gradually by all the the denoising techniques. 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.
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.

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 −3 ) 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.
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.

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 −4 have been found to be optimal.
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.
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.

Vascular network extraction
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 vessellike 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 Figure 9. Central transverse slice of the phantom obtained using (a) modified FDK algorithm, (b) Gaussian smoothing filter, (c) median filter, (d) wavelet with soft thresholding algorithm, (e) BM3D and (f) ITV denoising method. Diagonal streaking artefacts became more apparent with the median, BM3D and ITV techniques, presumably due to their edge-preserving characteristics. Table 7. Numerical results (standard deviation (SD), signal-tonoise ratio (SNR) and contrast-to-noise ratio (CNR)) obtained for the physical phantom.

Method
Region 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.
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 ]. Signal profiles of the ex vivo rat heart reconstructed with FDK (light blue), Gaussian smoothing (light green), median filter (magenta), wavelet with soft thresholding (dark green), BM3D (orange) and ITV (red) denoising algorithms in the transverse section (g) and in the sagittal section (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 ( 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.

Study limitations
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.

Conclusions
In this paper, we compared five major post-reconstruction denoising approaches suitable for μCT data acquired from commercial scanners. The methods were evaluated quantitatively on simulated and experimental data from a commercial μCT system (Sky-Scan 1172). Of the different techniques evaluated, the ITV algorithm had the best performance in terms of SNR and CNR for both simulated and experimental data.
The ITV method decreased the network reconstruction errors globally compared to the Feldkamp reconstructions including the smoothing option available as standard in the commercial software of the μCT system, but also compared with wavelet and BM3D filtering. Applied directly on the Feldkamp solutions, the ITV denoising algorithm improved the uniformity of the reconstructed tomographic images, and, importantly, facilitated an automated and accurate extraction of the 3D vascular network of an ex vivo rat heart. This 3D network reconstruction was not possible using any of the μCT data sets reconstructed with FDK, Gaussian smoothed or median filtered, due to the high level of noise. The ITV method's ability to preserve small structures in the 3D reconstruction was demonstrated in a synthetic vascular network and ex vivo data.