Block matching 3D random noise filtering for absorption optical projection tomography

Absorption and emission optical projection tomography (OPT), alternatively referred to as optical computed tomography (optical-CT) and optical-emission computed tomography (optical-ECT), are recently developed three-dimensional imaging techniques with value for developmental biology and ex vivo gene expression studies. The techniques' principles are similar to the ones used for x-ray computed tomography and are based on the approximation of negligible light scattering in optically cleared samples. The optical clearing is achieved by a chemical procedure which aims at substituting the cellular fluids within the sample with a cell membranes' index matching solution. Once cleared the sample presents very low scattering and is then illuminated with a light collimated beam whose intensity is captured in transillumination mode by a CCD camera. Different projection images of the sample are subsequently obtained over a 360° full rotation, and a standard backprojection algorithm can be used in a similar fashion as for x-ray tomography in order to obtain absorption maps. Because not all biological samples present significant absorption contrast, it is not always possible to obtain projections with a good signal-to-noise ratio, a condition necessary to achieve high-quality tomographic reconstructions. Such is the case for example, for early stage's embryos. In this work we demonstrate how, through the use of a random noise removal algorithm, the image quality of the reconstructions can be considerably improved even when the noise is strongly present in the acquired projections. Specifically, we implemented a block matching 3D (BM3D) filter applying it separately on each acquired transillumination projection before performing a complete three-dimensional tomographical reconstruction. To test the efficiency of the adopted filtering scheme, a phantom and a real biological sample were processed. In both cases, the BM3D filter led to a signal-to-noise ratio increment of over 30 dB on severe noise-affected reconstructions revealing original—noise-hidden —image details. These results show the utility of the BM3D approach for OPT under typical conditions of very low light absorption, suggesting its implementation as an efficient alternative to other filtering schemes such as for example the median filter.


Introduction
Optical imaging represents a set of techniques to study biological processes over time, often at the cellular and molecular levels.Until recently, optical scattering was a major impediment to imaging samples thicker than a few hundreds microns limiting their possible applications at the cellular level or for very small organisms' investigations.While some new techniques such as multiphoton microscopy (Denk et al 1990) or optical coherence tomography (Huang et al 1991) actively or passively reject scattered photons offering therefore higher resolution and better penetration depths, other tomographic methods such as fluorescence molecular tomography (FMT) (Ntziachristos 2006) or mesoscopic fluorescence tomography (MFT) (Razansky et al 2009) implement advanced mathematical models oftissue photon propagation to reconstruct fluorescence and absorption imaging contrast at the centimeter and millimeter scale, respectively.Both these tomographic techniques imply the formulation of a mathematical inverse problem in analogy to the other more conventional tomographic imaging modalities.While x-ray computed tomography (X-CT) implements straight-line integral projections (Radon 1986), FMT and MFT rely instead on diffusive models and depending on the light path shape, the reconstruction algorithms can involve analytical (Davis 1996), iterative (Herman and Lent 1976), deterministic (Roerdink 1998) or stochastic (Sheng and Ying 2005) procedures.Unfortunately, the price to pay in exchange of the deeper penetration depth is a substantial decrease in resolution.Optical projection tomography (OPT) (Sharpe et al 2002) is a recently introduced imaging technique which allows us to image small biological samples, such as insects (McGurk et al 2007), embryos (Sharpe 2003) or organs from small animals (Alanentalo et al 2008, Vinegoni et al 2010) ex vivo and at high resolution.Its applications spread from anatomical and histological analysis (Oldham et al 2007) to tissue proteins expression distribution (Sharpe 2003), from developmental biology to gene functions and recently to inflammation disease studies (Vinegoni et al 2010).One of the most attractive features of OPT is its capability to image entire biological samples up to a few centimeters in size at very high spatial resolution.To obtain such a degree of accuracy, the sample under investigation is made optically transparent through a chemical clearing process.In this way, its scattering and absorption properties are highly reduced, making the light diffusive contribution negligible.The sample is then illuminated with a light collimated beam intensity of which is captured in transillumination mode by a CCD camera.Fluorescence or absorption projection images of the cleared sample are taken every degree over a 360° full rotation in an X-CT analog fashion.Through the use of an optical imaging system with high telecentricity, only photons travelling parallel to the optical axis are projected on the pixels of the imaging camera.The technique therefore represents an optical analog of X-CT and the mathematical reconstruction problem, at least for the case of absorption OPT, can be solved in a similar fashion using the parallel-beam filtered backprojection (FBP) algorithm.Fluorescence tomographic absorption reconstructions in contrast utilize a Born-normalized method that relies on a normalized transillumination approach (Vinegoni et al 2009b).
While OPT has been demonstrated for both fluorescence (Sharpe 2003, Oldham et al 2007) and molecular imaging (Vinegoni et al 2010) we here concentrate on absorption reconstructions.Such reconstructions are valuable because they provide morphological features in addition to the functional information acquired through the respective fluorescence measurements; unfortunately they inevitably suffer of a reduction in the obtained image quality due to different physical limitations.In fact the chemical clearing process while providing for both light scattering removal and index matching, at the same time often leads to a negligible absorption coefficient in the treated sample.As a result, the transmitted optical signal will not fill the dynamic range of the CCD camera and because the signal-to-noise ratio (SNR) of an image depends on both the power of the true signal and the power of the noise, the acquired image will therefore present a very low SNR.The main disadvantages of working with such low SNR images reflect first in a loss of all the anatomical features embedded within the noise and second in the noise propagation during FBP, which make the reconstructions not suitable for further processing.
One approach to reduce the noise component consists in increasing the number of acquisitions per projection while taking the average.The drawbacks are twofold: first, this strategy increases the total acquisition time limiting all instances of high-throughput imaging, and second a prolonged illumination can lead to the bleaching of the fluorescence signal that will reflect on the eventual concomitant fluorescence reconstructions.
In this paper, we propose an alternative approach focusing on the implementation of a high performance filter based on a block matching for 3D collaborative filtering (BM3D) (Dabov et al 2007) for absorption OPT image quality improvement.Specifically, wedeal with the removal of the noise contributions that are due to the random fluctuations superimposed on the noise-free images.Note that in this study the CCD noise source contribution was assumed to be normally distributed, a condition easily obtainable with common lamp illuminators due to the sample high transparency.To test the performances of the BM3D filter for OPT, we used a phantom with specifically designed optical properties; this feature allows for a better testing of the filter effects under controlled SNR conditions.After the filter characterization we applied the algorithm for the case of a biological sample (mouse embryo), and the obtained reconstructions were compared with the corresponding histological sections.In addition a comparison with a standard filter such as the median filter was also used in support of the BM3D filter image improvement.

Material and methods
The experimental setup is shown in details in figure 1.A white light source beam (HL250-AY; AmScope), after passing through a beam expander and a combination of one or two interference filters, was used to illuminate the sample under investigation.In one configuration, only one filter (D525/50 m; Chroma, VT) centered at 525 nm and with a narrow full width at half maximum (FWHM) of 50 nm was present.In the other arrangement, a combination of a shortpass (FES0650; Thorlabas, NJ) and a long-pass (FEL0450; Thorlabas, NJ) filters was used in order to select a broadband (about 200 nm) of the optical spectrum.In addition filters with different optical densities (NDxxx; Thorlabs NJ) were used in order to keep the CCD values 10% below saturation.A beam expander with a combined two-lenses Galilean telescope (Thorlabs, NJ) and a diffuser (10DIFF-VIS; Newport, CA) were employed to make the source homogeneous over the sample's field of view.In order to collect the projections, the sample was rotated over 360° by way of a rotational stage (PR50; Newport, CA).The unitary step of rotation was equal to 1° with an absolute accuracy of 0.05°.Images were taken with a high resolution (1392 × 1024 pixel) 12 bit camera (Pixelfly QE; PCO Germany).
To test the filter we imaged both a phantom and a whole body mouse embryo as well.The phantom allowed us to control the effects of the filtering process while the mouse embryo provided for a real-life biological sample.
The phantom consisted of a pure agarose cylinder designed in order to provide images with different degrees of dynamic range and SNRs.This is a necessary feature in order to compare the reconstructions obtained from the BM3D-filtered projections acquired under the low SNR condition (figure 2(f)), with the reconstructions calculated starting from the high SNR images (figure 2(c)).A cylinder with a diameter of 8 mm was fabricated with a 1% agarose gel in which red india ink was diluted.Once solid the phantom was dehydrated in ethanol and cleared by immersion in a BABB solution (1:2 benzyl alcohol, benzyl benzoate ratio).Due to the porosity of the agare, the dye molecules diffused slightly within the entire cylinder creating a radially diffusion gradient.Because the ink presents an absorption peak at around 520 nm (see figure 2(b)), the presence of a narrow bandpass optical filter centered at 525 nm and placed in front of the illumination lamp (figure 2(a)) spectrally selects the component of the white light spectrum that is measured in transillumination mode and that is the most absorbed by the sample, creating a projection image with a high dynamic range and a high SNR.The replacement of the filter with a combination of a long-pass and a short-pass filters (figure 2 (d)) increases instead the bandwidth of the white light spectrum detected by the CCD giving rise to an image with a reduced dynamic range and low SNR that mimics a biological sample imaging condition (i.e.very small absorption coefficient).
In this way it is possible to image the same phantom and artificially modify the dynamic range and the associated SNR of the acquired images while keeping constant the images background and without changing any acquisition parameters (such as CCD's gain, integration time, etc) but only acting on the illumination lamp intensity by way of neutral optical density filters.Reconstructions obtained starting from high SNR projections were then compared to reconstructions calculated using the data post-processed with the use of the BM3D filter.Note that while the total amount of absorption changes when switching from one configuration to another, the overall profile remains the same allowing for a direct comparison between low and high SNR reconstructions.Finally, because in real-world case situations we are particularly interested in rendering the morphological features that are present within the biological samples, we therefore decided to disperse absorptive microparticles in the agarose phantom to mimic the typical image high frequency components present in common biological tissue.With this approach, sharp edges are introduced within the phantom and both low and high frequencies in the image power spectrum can be simultaneously simulated.
For what concern the biological sample, a mouse embryo at 10.5 days of development after fertilization was fixed in a 4% PFA (paraformaldehyde) solution at 4 °C temperature for 1 h followed by embedding in a 1% agarose gel and successively dehydrated through a series of immersions in 20-100% ethanol solutions.The clearance process was performed by immersion in a BABB solution for 2 days.At the end of the chemical treatment, the sample presented very low absorption with associated low SNR projections.

Denoising filter
In optical projection tomography, transillumination images are typically acquired at high light source intensities to minimize the camera noise contribution.For these levels of illumination, the camera noise can be modeled as normal distributed, additive, with zero mean and with unknown variance σ 2 .The captured image z(x) is therefore equal to where x are the image's pixels spatial coordinates, while y(x) and n(x) represent the noise free image signal and the camera noise, respectively.Although σ 2 is unknown, a number of accurate methods are available in order to estimate such a parameter from noise affected images.Herein we employed a wavelet domain estimation method (Donoho and Johnstone 1994).
To reduce the amount of noise present in the images, we used a high performance spatially adaptive Block matching 3D filter (BM3D) which is valid for normally distributed noise removal.BM3D is based on the assumption that a noise-free image spectrum of similar image fragments group can be better approximated as a combination of a few spectrum elements than a single image fragment (Katkovnik et al 2009).This guarantees that the signal energy is distributed to a small amount of coefficients, leading to a high signal sparse representation.Real images cannot always be well represented by a fixed 2D transformation and due to their huge varieties such transformations can only achieve sparse representations of certain patterns.BM3D filtering permits us to enhance the sparsity representation of an image by grouping, increasing the noise removal while at the same time preserving image features like localized details (e.g.sharp edges) and it has recently been demonstrated to produce very high-quality images (Dabov et al 2007), with few image artifacts.This is a critical factor for an OPT implementation because due to the FBP operation, any image distortion could result amplified in the final processed data (Walls et al 2005).
The basic operating principle behind the non-local filtering technique is based on the fact that 3D transformations of high correlated data can reach a highly sparse representation by grouping similar 2D image fragments into a stack in a 3D manner to form a so-called group.The signal is therefore transformed into a new domain where it can be well approximated as a linear combination of a few basis elements.In this way, the separation between the true signal and the noise components can be easily achieved by amplitude thresholding.This so-called transformed signal shrinkage preserves the few high-magnitude transform coefficients which should be related to the noise-free image power spectrum components, while at the same time discarding the low-magnitude ones which are essentially due to the noise.This approach is known as collaborative filtering (Dabov et al 2007) and is graphically illustrated in figure 3.Here in the whole noise clearance process involves the application of the BM3D filter to each projection image and the application of a ramp filter prior to the computation of the filter backprojection reconstruction as illustrated in figure 4.

Algorithm
The BM3D algorithm was introduced by Dabov et al and is schematically illustrated in figure 4(a).More in-depth details can be found in Dabov et al (2007).The first step is composed of several cascade blocks: grouping by a pre-filtered version of match blocking is performed on the noisy image in order to build a 3D stack of similar image fragments; discrete cosine transformation (DCT) is then used to map the blocks into a domain where the signal is sparse; shrinkage by hard thresholding is employed to attenuate the noise; the inverse 3D transform is used to map the processed data back to the image space returning noise-free image fragments; finally fragments are then repositioned to their original positions and since they may overlap, they are fused together by a weighted average (aggregation).The resulting image is named 'coarse' or 'basic estimate'.Step 2 is similar to step 1 and the major modification is in the way the groups are actually filtered.A Wiener filter is used during this second phase and the 'coarse estimate' image obtained in step 1 is exploited as a pilot spectrum.After Wiener filtering, 'inverse 3D', 'repositioning' and 'aggregation' processing blocks are sequentially performed as in step 1 but with different weights.More details about the theory of the BM3D filter can be found in the original paper of Dabov et al (2007) and a brief summary is reported in the appendix.
After BM3D, in order to obtain the absorption maps, a 1D ramp filter is applied to the image rows which are orthogonal to the rotation axis before applying the backprojection reconstructing algorithm; this processing leads to the final noise-free reconstruction estimation shown in figure 4(c).
The filtering code is written in the Matlab code and is in part based on the code made available at 'http://www.cs.tut.fi/~foi/GCF-BM3D/~ref_software'.

Results
The noise removal filter was tested first on a phantom as described above and applied directly on the acquired projections before obtaining any reconstruction via FBP (figure 4(d)).In order to do this, we tomographically reconstructed the absorption maps of the cylindrical phantom starting from the images collected at high values of SNRs (i.e.narrow bandpass filter inserted, figure 2(a)).
An axial reconstruction with an absorption profile taken along the cylinder's diameter AA′ (dashed line) is shown in figure 5(a) as a reference.As can be clearly seen, both low and high frequency components are in place.The low frequencies correspond to the phantom's radial decreasing optical absorption (shown in detail at the bottom of the reconstruction), while the high frequencies components are due to the presence of the sharp dots corresponding to the dispersed microparticles.The absorption map is then used as a reference in order to understand if the filter processing scheme removes any essential feature present in the noisy images.
Figure 5(b) shows instead an axial tomographic reconstruction taken at the same height of the phantom as in (a) with a wide bandpass filter inserted in the optical path (figure 2(e)).Due to the reduced SNR, the noise degrades the image severely making it very difficult to obtain quantitative information and to distinguish any underlying structure.Figure 5(c) shows the same axial tomographic reconstruction as in (a) and (b) but this time obtained after the application of the BM3D filter to each single projection.In the resulting reconstruction, the circular edge can be seen to be well preserved and it is worth to note the very low oversmoothing.Both the tomographic reconstruction and the absorption profile are very close to the reference one (figure 5(a)).For comparison, we show in figure 5(d) an axial reconstruction obtained after the application of a 2D median filter.By using this filter, the results are similar to the BM3D case but with a pronounced oversmoothing.In fact, while the absorption profiles are comparable, the axial reconstruction appears faded.This is more evident when considering the high frequency components of the image.Figure 6 shows a magnification of the axial reconstructions of figures 5 (a), (c) and (d) over an area where several microparticles are present for the high SNR (a), the BM3D filter (b) and the 2D median filter (c) case, respectively.It is evident that the phantom's perimeter is sharper in figure 6(b) than that in figure 6(c).For the median filter in addition, all spots corresponding to the absorptive microparticles are strongly attenuated and wider rendered.Worth to note are the artifacts present over the reconstructed spots on the median filtered data.A ring-like shape seems to be visible around the filtered spots, a phenomenon occurring while processing noisy data.Moreover, ring dimensions are correlated to the size of the filter window (i.e. the wider the filter window, the larger the ring diameter).
The filter performance was then tested on a biological specimen (i.e.low dynamic range within the projection image), specifically a mouse embryo.The sample, when chemically cleared, presented very low absorption optical properties, a feature quite common in early stages embryos.This contributes in making it very difficult to obtain absorption OPT reconstructions making it necessary to post-process the acquired data with different filtering options.
Figure 7 shows a comparison between filtered and unfiltered data for a mouse embryo.Figure 7(a) indicates the level at which the reconstructed sections were computed.Figures 7(b)-(d) correspond to the absorption reconstructions of the noisy, the BM3D filtered and the median filtered projections, respectively.As clearly evinced from the reconstructions, the BM3D filter was able to reduce the noise without distorting the original features of the noise-free image.
Sharp edges and indistinguishable embryo's structures are preserved and further enhanced.Moreover, no apparent artifacts were introduced by the filtering procedure.Looking at the median filtered results, the noise has been removed as well but at the expenses of the high frequency image components, which have been hardly attenuated.For this reason the resulting image, like in the phantom's case, appears severely blurred.A histological section at approximately the same level is provided for a further comparison indicating a high morphological correlation with the BM3D filtered reconstruction.
In order to quantitatively characterize the action of the BM3D filter on the tomographic reconstructions we compared the high SNR absorption maps with the reconstructions obtained from the low SNR projections with and without BM3D filtering (figures 5(c) and (b), respectively).In addition, to test its goodness with respect to other common filters, we compared the BM3D with a classic 2D median filter (figure 5(d)).The median filter window size was selected in such a way to obtain the same standard deviation of the error component for both filtered reconstructions.This is defined as the deviation of the absorption map, obtained from the BM3D filtered projections, from the absorption maps obtained from the projections with high values of SNR, after proper rescaling.The window was found to be 15 × 15 pixels wide.
For what concern the biological specimen, because there is no a priori knowledge of the absorption profile we computed the median filter window size using the MAD wavelet domain method (Donoho and Johnstone 1994), applied to both filtered reconstructions, as an error estimator.In a similar fashion as for the phantom's case, specimens' images were separately filtered with a standard 2D median filter for comparison.The filter window size (11 × 13 pixels wide) was chosen in such a way to obtain an error of the median filtering approximately equal to the one of the BM3D case.Filtering times were approximately 5 s for the median and 17 s for the BM3D filters.These values refer to only one projection (1392 × 1024 pixels) as processed by a Core 2 Duo processor at 2.53 GHz.A complete (360 projections) BM3D filtering requires approximately 1 h 40 min, while a FBP tomographic reconstruction takes approximately 8 h.Using graphic accelerated cards (Vinegoni et al 2009a) such a computational time could be reduced to a few seconds.
Table 1 summarizes the estimated SNR dB and the noise standard deviation for both the phantom and the mouse embryo for the cases of the noisy reconstructions and for those having applied the median and the BM3D filter, respectively.The BM3D filter was able to reduce the noise standard deviation by about 38 times in the phantom and over 30 times in a real imaged biological sample.The SNR increased about 30 dB after filtering for both cases.Although both the BM3D and the median filters are able to drastically reduce the noise, the major asset of the BM3D filter relies in its ability to preserve the images' features as illustrated above.

Conclusions
In this study we demonstrate the use of a random noise removal algorithm to increase image quality reconstructions in absorption optical projection tomography.The block matching 3D filter when applied separately to all optical projections was able to increase the signal-to-noise ratio over 30 dB while preserving original features.
The significant improvements obtained in the images' quality make this filter highly suitable for absorption OPT specifically in cases where little absorption and low dynamic range within the projection image are present.This is particularly true for early mammalian embryos.Due to the long computational time necessary to filter all the projections, future work will be aimed at implementing the filter on a low-cost parallel processing card making it suitable for high throughput imaging (Vinegoni et al 2009a).Because the fluorescence signal can be very faint, particularly for those cases where the CCD integration time has to be kept low in order to speed up the acquisition process or to reduce the photobleaching, we also expect the filter to significantly improve the fluorescence OPT reconstructions.

Figure 1 .
Figure 1.(a) Experimental setup.WLS, white light source; Sh, shutter; BE, beam expander; F, front sample band pass filters; S, sample; ODF, neutral optical density filter; TL, telecentric lens; CCD, imaging camera.(b) The imaged animal embryo is embedded in an agarose cylinder for rotational imaging.(c) The imaging phantom consisted of an agarose cylinder with radially diffused red ink.

Figure 2 .
Figure 2.Comparison between the two imaging modalities (a), (b) exploited to obtain high (c) and low SNR (f) projections.The use of a narrow band pass filter BPF1 centered at 520 nm and with 50 nm bandwidth (b) allowed to image only in the high absorption condition leading to a high dynamic range within the acquired projections (c).Instead, the use of the combination of a short-pass SP and a long-pass LP filters (d) led to a higher bandwidth (e) with a reduced dynamic range (f).

Figure 3 .
Figure 3. Graphical representation of the hard thresholding collaborative filter acting on a noise image.

Figure 4 .
Figure 4. (a) Block scheme of the BM3D filter.The first step allows us to compute a coarse noise-free image.The second step employs the coarse estimate in the grouping and in the Wiener collaborative filtering.Example of a specimen projection before (b) and after (c) BM3D noise removal.(d) Block scheme for the absorption OPT backprojection algorithm.

Figure 5 .
Figure 5. Axial tomographic reconstructions of the cylindrical phantom with red ink inserted.All reconstructed sections (a)-(d) are orthogonal to the phantom's axis and are all taken at the same height.The absorption profiles along the cylinder's diameter (dashed red lines) are here presented at the bottom of each reconstruction.(a) High SNR reconstruction obtained by way of the use of a narrow spectral optical band pass filter (figure 2(a)) for dynamic range increasing.(b) Low SNR reconstruction obtained with the configuration as indicated in figure 2(d); the circular shape is visible because of image range rescaling; the attenuation profile, scaled as in (c) and (d), is completely hidden by the noise.Reconstructions obtained after BM3D (c) and 2D median filter denoise (d) applied on the low SNR dataset.Horizontal and vertical scales' units are cm and cm −1 respectively.

Figure 6 .
Figure 6.Comparison between a high SNR (a), BM3D filtered (b) and median filtered (c) reconstruction.Details (spots and sharp edges) are well preserved with BM3D filtering but not with the median filter approach.

Figure 7 .
Figure 7.Comparison between different filters for the reconstruction of a mouse embryo.(a) Tomographic reconstruction of a mouse embryo obtained with absorption OPT using the BM3D filter.The plane indicates the position of the axial slice taken into account for the filters' comparisons.(b) Native reconstruction by parallel-beam FBP (360 projections and 1° step angle) without noise suppression.(c) Reconstruction obtained after BM3D filter denoising.The filter leads to very good reconstruction, preserving the high frequencies of the image spectrum, contrary to the 2D median filter shown in (d).(e) Histological section at approximately the same height.All reconstructions are presented with a negative colormap for visual details enhancing.White bar, valid for (b)-(e), corresponds to 0.4 mm.

Table 1
Estimated noise standard deviation and SNR dB for both the phantom and mouse embryo.