Patch-based anisotropic diffusion scheme for fluorescence diffuse optical tomography—part 2: image reconstruction

Fluorescence diffuse optical tomography (fDOT) provides 3D images of fluorescence distributions in biological tissue, which represent molecular and cellular processes. The image reconstruction problem is highly ill-posed and requires regularisation techniques to stabilise and find meaningful solutions. Quadratic regularisation tends to either oversmooth or generate very noisy reconstructions, depending on the regularisation strength. Edge preserving methods, such as anisotropic diffusion regularisation (AD), can preserve important features in the fluorescence image and smooth out noise. However, AD has limited ability to distinguish an edge from noise. We propose a patch-based anisotropic diffusion regularisation (PAD), where regularisation strength is determined by a weighted average according to the similarity between patches around voxels within a search window, instead of a simple local neighbourhood strategy. However, this method has higher computational complexity and, hence, we wavelet compress the patches (PAD-WT) to speed it up, while simultaneously taking advantage of the denoising properties of wavelet thresholding. Furthermore, structural information can be incorporated into the image reconstruction with PAD-WT to improve image quality and resolution. In this case, the weights used to average voxels in the image are calculated using the structural image, instead of the fluorescence image. The regularisation strength depends on both structural and fluorescence images, which guarantees that the method can preserve fluorescence information even when it is not structurally visible in the anatomical images. In part 1, we tested the method using a denoising problem. Here, we use simulated and in vivo mouse fDOT data to assess the algorithm performance. Our results show that the proposed PAD-WT method provides high quality and noise free images, superior to those obtained using AD.

calculated using the structural image, instead of the fluorescence image. The regularisation strength depends on both structural and fluorescence images, which guarantees that the method can preserve fluorescence information even when it is not structurally visible in the anatomical images. In part 1, we tested the method using a denoising problem. Here, we use simulated and in vivo mouse fDOT data to assess the algorithm performance. Our results show that the proposed PAD-WT method provides high quality and noise free images, superior to those obtained using AD.
Keywords: fluorescence diffuse optical tomography, image reconstruction, anisotropic diffusion, nonlocal means, structural information, multimodal imaging, regularisation S Online supplementary data available from stacks.iop.org/PMB/61/1452/ mmedia (Some figures may appear in colour only in the online journal)

Introduction
Fluorescence diffuse optical tomography (fDOT), also known as fluorescence molecular tomography, is an optical imaging modality that uses near-infrared excitation light sources to obtain fluorescence emission measurements of biological tissue, mostly small animals (Ntziachristos 2006, Stuker et al 2011, Darne et al 2014. Detection can be performed using a charged-couple device (CCD) camera, placed opposite the source, that is rotated around the subject of study. These tomographic measurements are used to recover three-dimensional (3D) images of the fluorescence distribution. Many fluorescence dyes have been developed, and many are commercially available (PerkinElmer Inc., Life Technologies), which can be used to label proteins, small molecules, antibodies, cancer cells, in order to non-invasively monitor disease progression and response to potential therapeutics, inflammation, glucose uptake, skeletal changes, infection progression, just to name a few of many applications.
Due to the diffusive nature of light when propagating through biological tissue, the image reconstruction is an ill-posed problem and images have low spatial resolution, which is one of the main limitations and challenges to overcome in fDOT. Regularisation methods that impose a priori constraints on the solution can be used to reduce the ill-posedness of the problem and obtain a stable and meaningful solution. Furthermore, the use of a priori anatomical information is known to improve the accuracy of the reconstructed images significantly. This information can be provided by a high resolution anatomical imaging modality, such as x-ray computed tomography (XCT) .
It is common to use the so-called L2 norm regularisation, which tends to produce smooth solutions, thus blurring edges in the images. Recently, a lot of attention has been drawn to L1 norm regularisation methods, which have the ability to smooth out noise while preserving edges in images, but are computationally more complex (Freiberger et al 2010a, 2010b, Dutta et al 2012. In Correia et al (2011) we proposed to use a nonlinear anisotropic diffusion regularisation method that can incorporate anatomical information and showed that spatial localisation and size of fluorescence inclusions can be accurately estimated. This method is based on edge detection techniques, where an image gradient is calculated and any gradient above a certain threshold is considered to be an edge. However, when noise levels are significant, noise may produce large gradients and, therefore, may erroneously be considered to be an edge. On the other hand, low contrast objects may be considered to be noise, and hence smoothed in the reconstruction process.
Similar to fDOT, Positron emission tomography (PET) and single photon emission computed tomography (SPECT) suffer from low resolution and greatly benefit from the inclusion of anatomical information in the reconstruction process. Bowsher et al (2004) introduced one of the most popular methods for incorporating anatomical information into PET/SPECT image reconstruction. This method promotes smoothing of radiotracer activity values in a set η, of fixed size, which consists of the most similar neighbours in the anatomical image. It is based on the assumption that neighbouring voxels with similar intensities in the anatomical image, within a search window W, are more likely to have the same radiotracer concentration than those with different intensities. However, the number of selected neighbours affects the quality of the reconstructed images, since irrelevant voxels may be included in η and important voxels neglected, causing an over or undersmoothing effect on the activity image. In order to overcome this limitation, Kazantsev et al (2012) proposed to use spatially variant set sizes based on anatomy and showed that it gives better results than choosing a constant set size. Furthermore, they used an edge preserving potential function, which penalises less the outliers than homogeneous regions, instead of the quadratic function used by Bowsher et al (2004), which favours smooth solutions. More specifically, they used an anatomically weighted anisotropic diffusion framework with Huber as the edge-preserving function.
In the original method proposed by Bowsher et al (2004), voxels with high similarity were assigned a weighting factor equal to 1 and the remaining 0. Other, more robust functions, can be used, which proved to be more accurate than a simple binary selection (Bousse et al 2010, Kazantsev et al 2011. Recently, the nonlocal means method (NLM), initially proposed by Buades for image denoising (Buades et al 2005), was used in PET image reconstruction as a regularisation method with (Chun et al 2012, Nguyen andLee 2013) and without anatomical information (Chen et al 2008, Wang andQi 2012). In the NLM method, pixels are averaged according to similarity between patches around pixels within a search window W, instead of a simple averaging strategy within a local neighbourhood. NLM exploits the redundancy of information within an image, by assuming that patches from different regions contain similar patterns and averaging them effectively reduces noise. NLM has superior denoising performance than local-based methods, but at the expense of higher computational complexity. To overcome this drawback, several methods have been proposed to accelerate the NLM without loss of denoising performance (Buades et al 2005, Mahmoudi and Sapiro 2005, Coupé et al 2008, Brox et al 2008, Orchard et al 2008, Tasdizen 2008. As a regulariser, NLM penalises less (anatomical) image patches that are very different and penalises more similar patches. Methods based on NLM regularisation have the advantage of being edge-preserving, able to preserve features within patches and do not require an explicit selection of the number of most similar neighbours.
Here, inspired by the recent developments in PET/SPECT and the NLM method, we propose an improved version of our nonlinear anisotropic diffusion regularisation method (AD) for fDOT image reconstruction . In our previous work we considered a local neighbourhood. Here, we consider a patch-based approach where we use robust edgepreserving weighting and potential functions. Patch dimensions are reduced using wavelet compression, not only to reduce computational complexity but also to increase the robustness of the method to noise. Therefore, the new method combines the advantages of NLM, AD and wavelet shrinkage methods (Donoho and Malik 1994). Furthermore, a priori structural information (SI) can be easily incorporated into the image reconstruction method. We refer to the new regularisation method as patch-based anisotropic diffusion with wavelet patch compression (PAD-WT). Since the AD and NLM methods were initially proposed for image processing, in part 1 we performed a thorough analysis of the performance of the PAD-WT method using a 2D denoising test problem (an ill-posed inverse problem). The aim of this study was to test the hypothesis that the PAD-WT method is superior to NLM, AD and that wavelet transform (WT) can be used to efficiently compress patches. Here, in part 2, we use PAD-WT as a regularisation function in fDOT image reconstruction and assess the efficacy of the method using simulated and in vivo mouse fDOT data.

Methods
A glossary of the notation and acronyms is presented in tables 1 and 2. ( where J − is the excitation source flux, Φ is the photon density, R is a boundary term that incorporates the refractive index mismatch, Γ is the boundary measurement on ∂ Ω, and n is the outer normal vector. The diffusion coefficient is given by where a µ is the absorption coefficient and s µ′ is the reduced scattering coefficient. The superscript e,f indicates the excitation and emission wavelengths, respectively. The fluorescence yield coefficient f is related to the quantum yield of the fluorophore and its concentration. For CCD camera measurements we have: where y is the measured data, M is a measurement operator that gives the data, Θ are the unknown source and detector coupling coefficients and the operator P represents the projection from the domain boundary ∂ Ω to the camera Σ (figure 1). The forward problem is computed numerically using the finite element method (FEM) on a tetrahedral mesh based on the geometry being considered, and then mapped into a Cartesian voxel grid (Schweiger and Arridge 2014). In the FEM, the domain Ω is divided into E elements, joined at S nodes. The solution of the diffusion equations is approximated by the piece- where u j are the basis functions. In the FEM framework equations (1)-(5) can be expressed as: K H,  (Soubret et al 2005). The normalised forward problem in fDOT is given by: where f Φ * is the solution to the adjoint diffusion equation for a source located on the boundary ∂ Ω at the detector position (Arridge 1999) and A is the Jacobian or sensitivity matrix. For a set of M s source-detector positions, where the detector is a camera with M M x y × pixels, it follows

Image reconstruction
The image reconstruction in fDOT consists in solving the problem: where α is the regularisation parameter and Ψ is a regularising functional that represents a priori information. (13) can be solved using quadratic or L2 norm regularisation, i.e. f 1 2 2 ∥ ∥ Ψ = , or iteratively using the split operator method with anisotropic diffusion regularisation introduced in Correia et al (2011) and summarised in algorithm 1: The solution is updated at each iteration k by a two-step algorithm. The first step, equation (14), is the Levenberg-Marquardt method, where λ is the damping factor and δ is the step length; and the second step, equation (15), is the nonlinear anisotropic diffusion method (Perona and Malik 1990), where Nit in specifies the number of inner or AD iterations,

Anisotropic diffusion. Equation
6, 2 neighbours per direction) and g is an edge-preserving function such as Huber, Tukey, Total Variation, Perona-Malik, etc (Correia et al 2011). Here, we use the Perona-Malik function defined by The parameter T is the threshold and can be selected using the normalised cumulative histogram (NCH) of the gradient (Perona andMalik 1990, Correia et al 2011). This method is commonly used in edge detection problems. The NCH indicates the probability ℘ of a gradient taking on a value less than or equal to the value X that the bin represents, i.e.
. It increases monotonically and the smoothness/sharpness of the curve indicates how smooth/sharp the edges Algorithm 1: Two-step image reconstruction algorithm.
are. The threshold can be calculated from the NCH by setting the threshold at, for example, 90 percent.
The algorithm runs until one of three stopping criteria is met, i.e. until the relative error ε is smaller than some given tolerance value, ε increases in comparison to the previous iteration or if the maximum number of outer iterations is reached Nit out . The relative error at the kth iteration between two consecutive iterations is defined as In a multimodality framework equation (15) becomes: where x w ij ( ) is a weighting factor related to the structural image x . It is an edge detection function similar to g that reduces diffusion across edges (coming from x).

Patch-based anisotropic diffusion.
We propose the following modification to the second step of the previous image reconstruction method: where w ij is the weight that measures the similarity between two cubic neighbourhoods (patches with fixed size P N N N = × × ) centred at voxel i and j, in a search window where W is the maximum distance from voxel i (refer to part 1 for further details). The weight w ij is defined as: where i N and j N are the neighbourhoods centred at pixel i and j, respectively. The parameter h is the smoothing parameter, which is related to the noise standard deviation σ and patch dimensions (Coupé et al 2008, Tasdizen 2008, Duval et al 2011, and C(i) is a normalising constant given by: and where f ip and f jp represent the p element of the patches f i ( ) N and f j ( ) N , respectively. Therefore, the weights w ij are large when patches are similar and small when they are much different.
Note that the patch-based anisotropic diffusion (PAD) method becomes the local AD when N = 1 and 1 = W . Here, we also propose to apply a WT compression method to reduce the dimensions of the patches and speed up the PAD method. Hence, the PAD-WT is our proposed regularisation method, where the weights in equations (17) and (18) are calculated using: where ip ϖ and jp ϖ represent the p wavelet coefficient of the wavelet transformed patches f i ( ) N and f j ( ) N , respectively. In the PAD-WT method, h depends on the size of the compressed patches ℓ, i.e. the number of wavelet coefficients kept to represent patches. These parameters are related through the following relationship: h 2 , where β is a constant. The constant β can be empirically tuned and σ is estimated using the NCH.
The new patch-based regularisation step is based on the NLM denoising method (Buades et al 2005). Refer to part 1 of this two-part paper for more technical details on the NLM, AD and proposed PAD(-WT) methods.

Evaluation
Simulations and in vivo mouse data were used to compare the new image reconstruction method using PAD-WT regularisation to our previous method using local AD  and a single-step linear reconstruction using quadratic regularisation (LR), which is equivalent to solving non-iteratively the first step of our two-step algorithm (see equation (14)).
Both PAD-WT and AD regularisation steps, equations (15) and (17), respectively, were solved iteratively. Note that the solution of the first step was computed on a tetrahedral mesh generated using the Matlab package Iso2mesh (Fang and Boas 2009), and then, since AD and PAD-WT are image-based methods, mapped into a regular grid covering the domain of the mesh. The number of outer iterations was set to Nit 150 out = and we chose the number of inner iterations Nit in required per outer iteration empirically. The time step was set to 1 τ = for both AD and PAD-WT methods. The parameter λ was defined as The parameter 0 λ and remaining parameters used in the algorithms are described in the results section. All reconstructions were run on a linux PC with an Intel Xeon E5-2665 CPU @ 2.4 GHz using Matlab Mex-files.

Simulation studies
We used the digimouse atlas (Dogdas et al 2007) to simulate a mouse with: (1) a lung disease, such as pulmonary inflammation (figure 2, movie 1 (stacks.iop.org/PMB/61/1452/ mmedia)); (2) tumour-like lesions within the lungs (figure 2). The generated finite element mesh had 63 591 nodes, 256 268 elements and dimensions 28 mm × 18 mm × 26 mm. The mesh was finer within the region of interest (ROI), i.e. lungs, and coarser elsewhere (see figure 3). We used a simplified two-tissue model and assigned different optical properties to the lungs ( 0.033 a µ = mm −1 and 2.2 s µ = ′ mm −1 ) and other tissue ( 0.01 a µ = mm −1 and 1 s µ = ′ mm −1 ). We calculated the optical properties using the parameters provided by Alexandrakis et al (2005) at a wavelength of 680 nm. We also used the atlas information to segment the lungs from the digimouse XCT images, to provide a realistic basis for the simulation and structural information (figure 4). In the first simulation, we assumed that the fluorescence signal was distributed throughout the whole lung structure with contrast f = 2 a.u. (arbitrary units) (figure 2, movie 1 (stacks.iop.org/PMB/61/1452/mmedia)). In the second study, we simulated two nearly spherical fluorescent targets embedded in the lungs, with contrast f = 2 a.u. and size of approximately 1.5 mm (figure 2). The centre of the targets was positioned at z = 0 mm.
In both simulations, the source and detector (CCD camera), placed opposite the source, were rotated around the chest of the digimouse. We calculated projection images of size 128 128 × for 36 source-detector positions evenly distributed around the mouse at z = 0 mm (figure 2). Data consisted of fluorescence and excitation projections corrupted by Gaussian random noise with a standard deviation equal to 20% of the mean value, and were used to test the hypothesis that our proposed method is robust to high levels of noise. Data space was compressed as in Rudge et al (2010), Correia et al (2013) by setting M 64 = ψ . We reconstructed images using a simple LR, our previous AD method  and our new PAD-WT method, with and without structural information.
We obtained results for 10 noise realisations (N 10 r = ) to analyse the bias-variance performance of the different algorithms for varying 0 λ values. We defined bias (or mean relative error) and variance as: where the solution space S was a regular grid of size 64 64 64 × × , the bias image and variance image where [ ] * is the value of voxel i in the reconstructed image for a given noise realisation r, f i true, is the value of voxel i in the true solution and We calculated each bias-variance curve for 0 λ = 0.01, 0.05, 0.1, 0.5, 1 10 2 { } × − . We chose the optimal 0 λ values for each reconstruction method to be the one that returns the lowest bias and variance, i.e. the point closest to the axis origin.
To further evaluate the quality of the reconstructed images we calculated the contrast to noise ratio (CNR) defined as: where the numerator is the mean value within a ROI of reconstruction f r [ ] * , which we considered to be all the voxels within the lungs (1st simulation) or inclusions (2nd simulation) and the denominator is the standard deviation of ROI c , which is the complement of ROI (background).

In vivo mouse study
The fDOT-XCT system (figure 5) used in this study is described in Schulz et al (2010). The fDOT instrumentation was integrated into a commercial micro-computed tomography system (eXplore Locus, GE HealthCare). The XCT system comprises an x-ray source and detector mounted on a rotating gantry. The fDOT system consisted of two laser diodes (680 nm and 750 nm, B&W Tek) and a back-illuminated cooled CCD camera (Pixis, 512B, Princeton Instruments). A filter wheel was positioned in front of the CCD to select the wavelength to be detected. The fDOT components were mounted on the XCT gantry orthogonal to the x-ray instrumentation axis, in a transmission geometry. The hybrid system enabled acquisition of both fDOT and XCT data over a full 360° rotation.
In the in vivo experimental study we used data obtained from an adult mouse with fluorescently labelled tumours embedded in the lungs. We used a double transgenic Kras +/− and BL6 Tyr −/− mouse model that developed lung cancer. The mouse was injected 24 h before imaging with a αvβ3 integrin-targeting fluorescence imaging probe (IntegriSense 680 from PerkinElmer) (Ale et al 2013). The integrin αvβ3 significantly upregulates in growing and metastatic tumours.
The fDOT-XCT system acquired excitation and fluorescence data images of size 512 512 × (resized to 128 128 × ) at 48 angular positions around the thorax area ( figure 6). For validation of the in vivo imaging results, the mouse was cryosectioned post-mortem in a cryotome and, in order to acquire RGB and fluorescence planar images of the cryosections, imaged using a multispectral-cryoslice imager (Sarantopoulos et al 2011). The cryoslice images used in this study correspond to the slices indicated by the arrows labelled S1 and S2 in figure 6.
The We reconstructed images using LR, our previous AD method  and our new PAD-WT method. The figure of merit used to analyse the quality of the reconstructions was the Dice similarity index (DSI) (Dice 1945), which is an overlap measure: where I seg and I true represent segmented images obtained by thresholding at 25% of the maximum value the images f * and f true , respectively. A DSI value close to 1 indicates a perfect overlap between I seg and I true , whereas a value of 0 indicates no overlap.

Simulation studies
In the simulation studies, we observed that the PAD-WT method converged faster to the final solution using a high patch compression, i.e. keeping a small number of wavelet coefficients of the compressed patches. After testing different levels of compression we found that setting P /8 ℓ ⌈ ⌉ = | | gave good results, where ⌈ ⌉ ⋅ denotes the ceil function. For all reconstructions we set 0.1 δ = and calculated the parameters T and h using the NCH and setting a threshold at 97% and 90%, respectively. The constant β was set to 1. For the PAD-WT method we set Nit 5 in = , whereas for the AD method we used Nit 150 in = .
4.1.1. Simulation 1. In the first simulation study, we used window sizes 5, 7 { } = W and patches of size N = {3, 7} in the PAD-WT reconstructions. Figure 7 shows the bias-variance curves obtained for the different image reconstruction methods, without and with SI (+SI). For all methods the lowest bias and variance is obtained with 1 10 0 3 λ = × − , except for the LR method where 5 10 0 3 λ = × − . Therefore, we used this optimal 0 λ for the reconstructions. Figure 8 shows the mean CNR obtained for the different image reconstruction methods, without and with SI, using the optimal 0 λ . In table 3 we can see the computational time for one AD and PAD iteration for two different window and patch sizes. For comparison, the PAD iteration time with N = 3, 7 = W and without patch compression was 81.32 s. For the larger patch size, N = 7, the iteration time was 968.80 s. Table 3 also shows the total reconstruction time of each method, with and without SI. Iteration time of step one (equation (14)) of the reconstruction algorithm was approximately 1.2 s. The computational time of the wavelet transform of patches with N = 3 and N = 7 was approximately 0.7 s and 5 s, respectively. Figures 9 and 10 show the fluorescence distribution images reconstructed from simulated data using LR, the two-step algorithm using AD and PAD-WT, with and without the , i.e. the two best reconstructions obtained using PAD-WT. Figure 9 (movies 2-7 corresponding to figures 9(b)-(g) (stacks.iop.org/PMB/61/1452/ mmedia)) shows the estimated fluorescence distribution as 3D isosurfaces at 25% of the maximum value. Additionally, 2D cross sections of the reconstructed volumes at z = −0.5 mm (to see the two lungs as separate structures), overlaid on the corresponding atlas cross section, are shown in figure 10. Figure 11 shows the bias-variance curves of the different reconstruction methods, without and with SI, that were obtained for the second simulation study. The PAD-WT reconstructions were performed using the patch and window sizes that returned the two best results in the previous study, i.e. using N = 3, 5 = W and 7 = W

Simulation 2.
. From the bias-variance curves we found 1 10 0 3 λ = × − to be the optimal value for the AD and PAD-WT methods, except for the PAD-WT with N = 3, 7 = W and SI for which the optimal value was 5 10 4 × − . For the LR method the estimated optimal 0 λ value was 5 10 3 × − .     Figure 12 shows the mean CNR obtained for the different image reconstruction methods, without and with SI, using the optimal 0 λ . Figure 13 shows the corresponding fluorescence distribution images reconstructed from simulated data using LR, the two-step algorithm using AD and PAD-WT (N = 3 and 5, 7 { } = W ), with and without the incorporation of SI. The images are 2D cross sections of the reconstructed volumes at z = 0 mm, overlaid on the corresponding atlas cross section. Figures 14 and 15 show the normalised fluorescence distribution images obtained from experimental mouse data and validation images corresponding to slices S1 and S2, respectively. Figures 14 and 15(a)-(d) show the images reconstructed using LR, the two-step algorithm with AD and PAD-WT for N = 3, 5 = W and 7 = W , without including SI into the reconstruction. Figures 14 and 15(f)-(h) show the images obtained when SI is used in the two-step algorithm with AD and PAD-WT for the two different window sizes. The validation cryoslices can be seen in figures 14 and 15(e), showing the planar fluorescence signals overlaid on the RGB images. The fluorescence peaks are represented in yellow. The cryoslices are from the thorax area (see figure 6) and display the bones of the rib cage, heart, lungs and (fluorescently labelled) tumours of different sizes within the lungs. Figure 16 shows the profile plots across the fluorescence cryoslice images and fluorescence images reconstructed using the AD and PAD-WT methods, corresponding to the dashed lines labelled P1 and P2 in figures 14 and 15(e).

In vivo mouse study
The obtained DSI are shown in figure 17. These results indicate the amount of overlap between the validation cryoslices and images reconstructed with the different methods. The parameters used in the reconstruction of mouse data were 1 10 0 3 λ = × − , 0.01 δ = , 4 ℓ = , 1 β = , T and h were calculated from the NCH by setting a threshold at 99% and 90%, respectively. For the AD method Nit 150 in = and for the PAD-WT method Nit 5 in = . We used the homogeneous sensitivity matrix A, since the results obtained showed the best agreement with the validation data.

Discussion
We proposed to modify our split operator method for fDOT image reconstruction by using PAD-WT instead of AD regularisation. We tested, using simulated and in vivo mouse data, the hypothesis that the proposed image reconstruction method is more robust to noise. The first simulation study consisted of a mouse with a lung disease, simulating an extensive inflammatory response as in, for example, asthma or due to inhalation of toxic particles, most commonly  by smoking. In the second simulation study, two tumour-like structures were placed within the lungs. Finally, in the in vivo study we used a mouse with several tumours lesions within the lungs.
The proposed method was motivated by the effectiveness and simplicity of the NLM method for image denoising. Therefore, in part 1 we used a denoising problem to compare our PAD-WT method with the AD method and iterative NLM. Results showed that the PAD-WT method provides better denoised images than the AD or NLM methods.
In the first fDOT simulation study, the fluorescent marker was distributed throughout the whole lung structure. We observed that our previous AD method  provides noise free fluorescence images (figures 9 and 10(b)) as opposed to the simple linear reconstruction method (figures 9 and 10(a)). However, our PAD-WT provides better estimates of the size and contrast of the simulated fluorescence distribution (figures 9 and 10(c)-(d)). Furthermore, like the AD method, structural information can be incorporated to further improve the reconstructions (figures 9 and 10 (e)-(g)). From figures 7 and 8 we can see that the images generated with the PAD-WT method (for 10 noise realisations) have the lowest bias, lowest variance and highest CNR, whereas the worst results, with images highly dominated by noise, were obtained using LR. As expected, incorporating SI into the reconstruction with PAD-WT resulted in quantitatively and visually superior images (figures 7-10).
The size of the search window and patches affects the performance of the PAD-WT method. However, the larger these are the slower is the method. We tested our method for the window sizes that returned the best denoising results (refer to part 1), i.e. 5, 7 { } = W , and set the patch length to N = {3, 7}, in order to find the most suitable window and patch size for fDOT. We observed that the bias and variance is lower and CNR is higher for N = 3 and 7 = W (figures 7 and 8), i.e. the method is more robust and recovers the fluorescence source location, size and concentration more accurately. However, computational time increases with window size, but by using patch compression the iteration time for the PAD-WT method with N = 3 and 7 = W was reduced to more than half (section 4.1 and table 3). Furthermore, figures 7 and 8 show that image quality improves by using 7 = W compared to 5 = W , but yet images are visually quite similar, particularly when SI is used, at the cost of more than twice the computational time (table 3). Therefore, in terms of reconstruction quality, we consider 7 = W to be the optimal window size. Alternatively, if speed is preferred over quality, then using 5 = W is a good option. For a patch of length N = 7, as used by Buades et al (2005), even though WT significantly reduces the computational time, the method is much slower and leads to worse qualitative results.
We observed that our method converged faster for higher patch compression levels, which is probably due to the denoising properties of wavelet thresholding. Therefore, PAD (or NLM) may consider two similar patches with different noise patterns to be different, while PAD-WT compares denoised patches and is more likely to succeed in this task.
In the second simulation study, the PAD-WT method with N = 3 and 5, 7 { } = W was used to reconstruct two lesions in the lungs. Results were compared to those obtained with the AD and LR methods. We observed that the best quantitative results were once again obtained with the PAD-WT method with 7 = W (figures 11 and 12). Nevertheless, as before, the images obtained using 5 = W and 7 = W are visually very similar (figure 13). The AD method provides cleaner images with less noise compared to the LR method ( figure 13). Still, the AD method fails to filter high contrast noise, even when SI is included in the image reconstruction. This is also reflected in the larger bias values, when compared to the PAD-WT method (figure 11). In theory, one could increase the threshold T to filter high gradient isolated points due to noise. However, if we do so, the method smooths out the edges and anatomical structures and fails to recover the fluorescent targets.
We further validated our method using experimental data obtained from a mouse with fluorescently labelled tumours within the lungs. The results show that LR returns a noisy image (figures 14 and 15(a)), whereas AD reduces the reconstruction artefacts, except at a few regions near the boundary (figures 14 and 15(b)). As mentioned previously, in the presence of high contrast artefacts localised within a single or few voxels, which in fDOT often occurs near boundaries, the AD method fails to generate noise free reconstructions, even when SI is used ( figure 15(f )). Our PAD-WT method overcomes this challenge and reconstructs a noise free image with the peak of the fluorescence signal in the lungs (figures 14 and 15(c)-(d)). The fluorescent lesions in the mouse lung are clearly visible in the images reconstructed using PAD-WT with SI (figures 14 and 15(g)-(h)) and are in agreement with the validation cryoslice images (figures 14 and 15(e)). The profiles P1 and P2 across the cryoslice and reconstructed fluorescence images corresponding to slices S1 and S2, respectively, were used to validate the results (figure 16). Even though the profile plots for the LR method (figures 16(a)-(b)) are quite reasonable, it can be difficult to distinguish between artefacts and relevant information in the image. The AD method provides quite good results, but overestimates the fluorescence yield. Figures 16(c)-(d) shows that a better correspondence was achieved between the cryoslice fluorescence profiles and the fDOT reconstructions obtained when resorting to SI. Furthermore, the profiles show that the results are structurally and quantitatively more accurate for the PAD-WT. The proposed method can recover the relative contrast and size of the lesion, ranging in size from 1 to 3 mm. These results emphasise the importance of SI for the accurate recovery of fluorescence inclusions.
The highest DSI values were achieved with the PAD-WT method ( figure 17). These values show that the fluorescence signals obtained with the PAD-WT, for both window sizes, are similar to those observed in the corresponding cryoslices, even though they are not exactly, but only approximately from the same region as the fDOT images. It should be noted that the location of the fluorescence signals and anatomy is only known approximately due to the inherent difficulty in keeping these in the exact same position between in vivo and ex vivo experiments. Hence, one cannot expect a perfect overlap between I seg and I true .
Background fluorescence, due to nonspecific binding of the fluorescent agents used for imaging, can be seen in the tissue surrounding the lung area. The fluorescent agent accumulates preferentially in tumour tissue, but the agent may be present throughout the body even if at lower concentrations, providing misleading information. These signals and inherent autofluorescence of tissues contribute to the reduction of image contrast, which may preclude the detection of small lesions, with low fluorescence contrast, a situation that occurs when imaging tumours at an early development stage. The results, shown in figures 14 and 15(g)-(h), demonstrate the robustness and effectiveness of our proposed method in removing artefacts due to the presence of background fluorescence.
Some of the other proposed methods that make use of structural information require segmentation techniques to classify voxels into different tissue types (Davis et al 2007. However, the drawback of these methods is that they require additional computational time and accurate segmentation is a difficult task. Furthermore, these ignore the internal structure of organs. Both AD and PAD methods have the advantage of using noise suppression and edge preservation constraints directly in the regularisation term. Regularisation strength is controlled by weighted differences between voxel intensities or patches, in the PAD method, within a search window. If structural information is available, then weights are calculated based on similarities within the structural image. The weighted (patch-based) anisotropic diffusion regularisation combines the advantage of preserving/enhancing the edges of the reconstructed fluorescence image while reducing the noise and stopping diffusion through boundaries defined by either the optical or structural images. In our studies, the complete non-homogeneous lung structure (figure 4) was used as a priori information in the image reconstruction algorithms and not just the organ boundary as in segmentation based methods.
One of the advantages of using the AD method is the computational speed, even though it requires more inner iterations than the PAD-WT method. Nevertheless, our patch compression method with 5 = W allows reconstructions within a reasonable time. As a matter of fact, the computational time per outer iteration is longer for the AD method. Furthermore, NLM-based methods can easily be parallelised and implemented to exploit GPU (graphics processing units) acceleration (Cuomo et al 2014). Alternatively, to reduce the computational times, a semi-implicit scheme can be used, instead of the explicit discretisation used here, so that τ can be relatively large without causing numerical instability, which is particularly advantageous to reduce the number of inner iterations .
Better performance of the image reconstruction method with PAD-WT is also expected if we use a more efficient and automated selection of the parameters δ, Nit in , ℓ, T, h and 0 λ . Alternatively, a Bayesian approach (Zhang et al 2013) can be explored where these parameters (of a prior distribution, known as hyperparameters) can be updated iteratively.
In this work, we map from the FEM mesh to a regular grid with no regard to the irregular boundary of the object. In future work, we will attempt to perform the reconstruction on the mesh to avoid error propagation resulting from these mappings.

Conclusion
In this work, we propose to improve our previous image reconstruction approach for fDOT based on a split operator method with AD regularisation.
As before, a reconstruction step alternates with the regularisation step, but a patch-based AD regularisation is used instead. Our previous image reconstruction method with AD relies on image gradients within a local neighbourhood to determine the regularisation strength at each image voxel and has limited performance under high noise conditions. Using our proposed patch-based approach we obtain a more robust estimation of the AD weights, since it measures voxel similarity through patches as opposed to single voxels. Moreover, structural information can easily be incorporated into the PAD-WT method.
As the AD method, the PAD-WT method has the ability to preserve important features in the reconstructions, which may be absent in the structural images, but is considerable more efficient in noise suppression. Moreover, the use of wavelet patch compression results in faster computational times.
Our method is very flexible, it works well with or without structural information, and it can easily be extended to other imaging modalities.