Translator Disclaimer
12 August 2013 Wavelet-based data and solution compression for efficient image reconstruction in fluorescence diffuse optical tomography
Abstract
Current fluorescence diffuse optical tomography (fDOT) systems can provide large data sets and, in addition, the unknown parameters to be estimated are so numerous that the sensitivity matrix is too large to store. Alternatively, iterative methods can be used, but they can be extremely slow at converging when dealing with large matrices. A few approaches suitable for the reconstruction of images from very large data sets have been developed. However, they either require explicit construction of the sensitivity matrix, suffer from slow computation times, or can only be applied to restricted geometries. We introduce a method for fast reconstruction in fDOT with large data and solution spaces, which preserves the resolution of the forward operator whilst compressing its representation. The method does not require construction of the full matrix, and thus allows storage and direct inversion of the explicitly constructed compressed system matrix. The method is tested using simulated and experimental data. Results show that the fDOT image reconstruction problem can be effectively compressed without significant loss of information and with the added advantage of reducing image noise.

1.

Introduction

Fluorescence diffuse optical tomography (fDOT) is an imaging modality that aims at reconstructing three-dimensional (3-D) distributions of fluorescent markers embedded within biological tissues.1 This technology has facilitated monitoring of molecular activity, tumor growth, response to drug therapy, etc., mostly in small animals.

In fDOT, a near-infrared (NIR) light source at excitation wavelength is illuminated onto the subject under study at different positions. Fluorophores absorb the excitation light and re-emit part of the energy at a longer wavelength. In most modern noncontact fDOT systems, light intensity measurements (CW) at the emission, and possibly, excitation wavelengths are collected by a charged-coupled device (CCD) camera placed opposite to the light source, which is rotated around the object of study.2

The image reconstruction consists of the inversion of a linear operator (sensitivity matrix), mapping the fluorescence yield of the fluorophore, which is linearly related to its concentration, to the measured data. However, due to the diffusive nature of light propagation in biological tissue the image reconstruction is an ill-posed inverse problem, meaning that noise in the data may give rise to significant errors in the reconstructed image.

The ill-posedness of the problem can be reduced by using large data sets.3,4 The use of CCD cameras produces large amounts of data, typically of the order of nmeas107. However, this makes the problem large scale and more challenging to solve. Furthermore, the parameters to reconstruct are so numerous, nrecon106, that the sensitivity matrix is extremely large. The large matrix size poses a problem for storage, as well as computation time required for inversion. Therefore, specialized algorithms that can handle problems with large dimensions are required.

Markel and Schotland3,5 proposed a series of algorithms for the reconstruction of images from extremely large data sets, which exploit symmetry properties of simple geometries. Konecky et al.6 reconstructed diffuse optical tomography (DOT) images using a fast inversion method, which combines plane wave illumination7 with analytic image reconstruction methods. The reconstruction method exploits the block structure of the linear operator that couples the measurements to the optical properties of the object of study, and hence, instead of inverting a large matrix, the problem is reduced to inverting a set of relatively small matrix blocks. Lukic et al.8 used a similar method, but replaced point sources by structured illumination (Fourier encoding), and hence, data are directly measured in the Fourier space. Since DOT suffers from low spatial resolution, only a few low Fourier components retain relevant information. In this approach, the inverse problem is fully Fourier encoded, i.e., sources, detectors, and solution space are represented in the Fourier space. Ducros et al.9 used a similar approach, but based on wavelet encoding methods. Ripoll10 proposed a hybrid approach, where only the measurements are Fourier encoded, while the sources and solution space are kept in the real space. However, most of these approaches are based on analytical expressions (Green’s functions) for the slab geometry.

For complex geometries, an alternative is to use a matrix-free Krylov-subspace approach, in which the storage of the sensitivity matrix is replaced by vector-matrix products, i.e., solving the forward and adjoint problems where the solutions are computed numerically using the finite element method (FEM).11 The matrix-free method relies on iterative Krylov-subspace methods, such as generalized minimal residual method, which may suffer from slow convergence. Therefore, the method overcomes the storage problems of large data sets at the expense of computational efficiency.

Multigrid methods offer the capability of reconstructing images at lower spatial scales by solving smaller instances of the problem, which allows approximate solutions to be computed quickly at low resolution and progressively refined. Zhu et al.12 developed a wavelet-based multigrid approach. The rows and columns of the sensitivity matrix are wavelet transformed and the system at a coarser resolution is solved. Solutions at coarser resolutions are used as the initial guess for solving systems at finer resolutions. It is also possible to select regions of interest at smaller scales based on coarse level reconstructions. However, in order to maintain reconstruction resolution, this approach requires reconstructions at many scales, and thus, does not result in significant reduction in computation time. Furthermore, it requires the full sensitivity matrix to be calculated and stored.

Another approach is to use the wavelet-Galerkin method to generate the forward model, instead of using the more time consuming FEM (Ref. 13). However, unlike FEM, this method is only applicable to simple geometries, but the fictitious domain method is used to overcome this limitation.14 In the wavelet-Galerkin method, the wavelet scaling functions are used as a basis on a regular grid. This allows a more compact representation of the problem and the use of multigrid techniques. As before, solutions can be obtained at progressively higher wavelet resolution scales, terminating when appropriate. Although this is an efficient way to compute the full matrix for large solution spaces, it is not suitable for large data sets.

Wavelets were also used in Ref. 15 to reduce the fDOT forward model computation time by projecting the original FEM matrices into a series of wavelet bases.

This paper proposes a method that reduces significantly the size of the sensitivity matrix by compressing its rows and columns using wavelets, allowing it to be stored explicitly while maintaining the spatial resolution and accuracy of the reconstructed images and low-computational time. This method is an extension of the work in Ref. 16, where wavelet compression is applied to a very large data sets in order to reduce data space dimensions. Here, a similar approach is used to reduce the solution space dimensions. Therefore, given the reduced dimensions of the image reconstruction problem, the solution can be found using direct inversion methods. The performance of the method is assessed using simulations and experimental phantom data. The effect of incorporating a priori structural information into the compressed image reconstruction is also analyzed.

2.

Methods

2.1.

Forward Problem

In CW-fDOT, the forward model is described by a set of coupled diffusion equations in a domain Ω:17

Eq. (1)

·κeΦe+μaeΦe=0,

Eq. (2)

·κfΦf+μafΦf=Φef,
with boundary conditions on Ω

Eq. (3)

Φe+2Rκen^·Φe=J,

Eq. (4)

Φf+2Rκfn^·Φf=0,

Eq. (5)

Γe,f=κe,fn^·Φe,f,
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 κ=(1/3)(μs+μa)1, where μa and μs are the absorption and reduced scattering coefficients, respectively. The superscript e and f indicates the excitation and emission wavelengths λe and λf, respectively. The fluorescence yield coefficient f is related to the quantum yield of the fluorophore η and its absorption coefficient at λf. For CCD camera based measurements we have

Eq. (6)

ye,f=M[Φe,f]=ΘP(ΩΣ)Γe,f,
where M is a measurement operator that gives the data, Θ is the unknown source and detector coupling coefficients and the operator P represents the projection from the domain boundary Ω to the camera Σ.

The forward problem is computed numerically using the FEM on a tetrahedral mesh based on the geometry being considered, and then mapped into a Cartesian voxel grid. In the FEM, the domain Ω is divided into P elements, joined at N nodes. The solution of the diffusion equations is approximated by the piecewise function Φe,f=j=1NΦje,fuj, where uj is the basis functions. In the FEM framework, Eqs. (1)–(5) can be expressed as

Eq. (7)

KeΦe=Q,

Eq. (8)

KfΦf=H,
where

Eq. (9)

Ke,f=Ω(κe,fui·uj+μae,fuiuj)dΩ+Ω12Ruiujd(Ω),

Eq. (10)

Q=Ω12RJuid(Ω),

Eq. (11)

H=ΩΦefuidΩ.
Therefore, the solution Φe,f can be found by inverting the matrix Ke,f.

Normalizing the measured fluorescence photon density M[Φf] by the measured excitation photon density M[Φe] reduces the effects of Θ.18 The normalized forward problem in fDOT is given by

Eq. (12)

y^=yfye=M[ΦeΦf*]M[ΦeΦe*]=Jf,
where Φe,f* is the solution to the adjoint diffusion equation for a source located on the boundary Ω,19 f is the fluorescence yield coefficient, and J is the Jacobian or sensitivity matrix. For a set of Ms source-detector positions, where the detector is a camera with Mx×My pixels, it follows that the measurements ye,f are vectors of size Ms×Mx×My and J is a matrix of dimensions (Ms×Mx×My)×N, where N is the dimension of the solution space that, once mapped into a regular grid, has dimensions Nx×Ny×Nz.

2.2.

Image Reconstruction Problem

The image reconstruction in fDOT consists in solving the problem

Eq. (13)

f=minimise[12y^Jf2+αΨ(f)],
where α is the regularization parameter and Ψ is a regularizing functional that represents a priori information. The previous equation can be solved using zero-order Tikhonov regularization, i.e., Ψ=(1/2)f2, or iteratively using the split operator method with anisotropic diffusion regularization and structural a priori information (ADSP) introduced in Ref. 20:

Eq. (14)

fk+1/2=fk+JT(JTJ+λI)1(y^Jfk),

Eq. (15)

fk+1=fk+1/2+ΔtL(fk+1/2)fk+1/2,
where the first step is the Levenberg–Marquardt method, where λ is the damping factor that changes at each iteration; and the second step is the nonlinear anisotropic diffusion method,21 where Δt is the time step and L is the anisotropic diffusion function given by L(f)=·[S(|fst|)g(|f|)]. Here, g(|f|)=1P(|f|X) is the exceedance function, which is the probability of an edge of interest being present and can be calculated using the normalized cumulative histogram (NCH) of the image gradient.20 The NCH indicates the probability P of a gradient taking on a value less than or equal to the value X that the bin represents, i.e., P(|f|X). In a multimodality framework, S(|fst|) is a weighting factor related to the structural information.20 It is an edge detection function that stops the diffusion process across edges. The additive operator splitting scheme is used to discretize the nonlinear anisotropic diffusion equation.20

2.3.

Fast Wavelet Transform

The fast wavelet transform (FWT) decomposes a signal or function f0 into different frequency subbands.22 It uses a low-pass filter h and a high-pass filter g to obtain the approximation Wc and detail Wd coefficients. The approximation coefficients represent the approximation of the signal at a resolution 2j, where j is an integer that specifies the resolution level. The detail coefficients contain the details of the original signal, i.e., the high-frequency information. For a signal f(n)=W(j+1,n), with n samples and at a starting scale j+1

Eq. (16)

Wc(j,n)=h(n)*Wc(j+1,n)2,

Eq. (17)

Wd(j,n)=g(n)*Wd(j+1,n)2.

Thus, the FWT consists of a convolution, followed by downsampling by a factor of 2 (2), i.e., keeping the even index samples. The wavelet transform can be implemented as a decomposition filter bank, where the initial signal goes through a series of filters. The synthesis filter bank can be used to perform the inverse transform (IFWT) and reconstruct the signal f0. First, the signal is upsampled by a factor of 2 (2), i.e., adding a zero between samples, followed by convolution with the inverse filter. Note that the filters are related to each other by g(n)=(1)nh(1n). Figure 1 shows a one-level decomposition and synthesis filter bank.

Fig. 1

One-level 1-D FWT and IFWT.

JBO_18_8_086008_f001.png

The separability property of the two-dimensional (2-D) wavelet transform means that performing a 2-D wavelet transform is equivalent to performing two one-dimensional (1-D) transforms, i.e., one 1-D transform along with the columns of the image f(m,n) (n-axis) and another 1-D transform along with the rows of f(m,n) (m-axis) (Fig. 2). Analogously, performing a 3-D wavelet transform is equivalent to performing three 1-D transforms. The same applies to the IFWD in higher dimensions.

Fig. 2

One-level 2-D FWT.

JBO_18_8_086008_f002.png

2.3.1.

Data compression

Data compression is used to reduce the dimensions of the data for computational efficiency,16 while simultaneously reducing the redundancy of the data. We apply a four-level 2-D FWT (W2-D) to the projection image yie,f captured by a CCD camera, where i=1,2,,Ms. The obtained wavelet coefficients are W={Wc(j3,x,y),Wd1(jk,x,y),Wd2(jk,x,y),Wd3(jk,x,y):j=7and0k3}, where Wd1, Wd2, and Wd3 are the detail coefficients that passed through a high-pass filter at least once. We keep the coefficients that are larger than a threshold TM, set to the Mψ largest coefficients (i.e., Mψ=W>TM), which contain most of the relevant information, and hence, the dimensions of the vectorized compressed data y˜ are reduced to Ms×Mψ. For example, if WT=128×128 is the total number of wavelet coefficients, then TM=Wsort(P×WT), where Wsort are the wavelet coefficients sorted in descending order and P is the percentage of coefficients to be removed. If P=0.98, then Mψ=328 and the compressed data y˜ has dimensions 18×328. Similarly, the size of the row compressed Jacobian J˜ is reduced to (Ms×Mψ)×N. The resulting compressed forward problem is y˜=J˜f˜.

2.3.2.

Solution compression

In this paper, we extend the work in Ref. 16 by applying a similar wavelet compression to each row of the Jacobian, thus reducing the dimensions of the solution space. Each row ri of the Jacobian represents the sensitivity of one measurement y^ (matrix J) or compressed measurement y˜ (matrix J˜) to changes in f. Both f and ri represent 3-D images of dimensions Nx×Ny×Nz, and therefore, a 3-D FWT (W3-D) is performed on the rows of the sensitivity matrix. Two methods are investigated, one where the sensitivity matrix is sparse and another where the matrix is fully compressed.

Sparse matrix: A highly sparse matrix J̆ can be obtained by keeping the wavelet coefficients V larger than a threshold TN, set to the NΨth largest coefficient in each row, and setting the remaining coefficients to zero

Eq. (18)

J̆ij={Vij,ifVij>TN0,elsewhere,
where the subscript i is the row index and j is the column index. For example, if WT2=64×64×64 is the total number of wavelet coefficients, then TN=Vsort(P×WT2), where Vsort are the wavelet coefficients sorted in descending order. If we consider P=0.98, then NΨ=5243.

The resulting forward problem, using data and solution compression, is y˜=sparse(J̆)f̆. The solution is a sparse image of size Nx×Ny×Nz and NΨ nonzero elements. The matrix J̆ has dimensions (Ms×Mψ)×(Nx×Ny×Nz) and (Ms×Mψ)×NΨ nonzero elements. The solution is obtained by solving the inverse problem and performing a 3-D IFWT (W3-D1) to f̆.

Fully compressed matrix: Alternatively, if only NΨ wavelet coefficients are kept, then the size of the compressed Jacobian J¯ is reduced to (Ms×Mψ)×Nψ. The forward problem is y˜=(J¯)f¯, where f¯ is a vector of size NΨ. Therefore, the solution is a compressed image of size NΨ, which is converted into a sparse Nx×Ny×Nz image and inversely wavelet transformed. Consider Ai=W3-D{ri} and AM=mean(Ai) (average along the rows), the wavelet coefficients kept correspond to the jN location of the NΨth largest AM elements: J¯i,j=Ai,jN. Unlike the sparse matrix method, where the jN location corresponds to the location of the NΨth largest coefficients in each row and, therefore, the jN locations are different for each row, in the fully compressed method the jN locations are the same in all the rows of the wavelet transformed matrix A.

2.4.

Computation Procedure for Compression of the Solution Space

The FWT is used to efficiently compute the compressed matrix J̆, without explicitly computing the matrix or a wavelet transform matrix. The matrix J̆ is calculated row by row; and hence, an FWT is performed on each row of the sensitivity matrix and only the NΨth largest components are kept. To compute the fully compressed matrix J¯, one needs to first explicitly compute the wavelet transformed matrix A of size (Ms×Mψ)×(Nx×Ny×Nz). The computational procedure is described in Algorithms 1 and 2 for the sparse matrix and fully compressed matrix, respectively.

Algorithm 1

Sparse matrix.

Initialize an empty sparse matrix J̆ of dimensions (Ms×Mψ)×(Nx×Ny×Nz)
fori=1,...,(Ms×Mψ)do
ri=M[ΦeΦf*]/M[ΦeΦe*]
R3-D=W3-D{ri}
R1-DR3-D
RssortR1-D in descending order (and store their location j)
Rt ← keep the NΨth largest components (and locations jN)
J̆i ← set elements jN of J̆i to their corresponding Rt values.
end for
Solve the inverse problem for f̆1-D:
Set tolerance value ε and maximum iteration number Nk
whiley˜J̆f̆1-D2<εk<Nkdo
f̆1-Dk+1/2=f̆1-Dk+J̆T(J̆TJ̆+λI)1(y˜J̆f̆1-Dk)
f̆1-Dk+1/2=f̆1-Dk+1/2+ΔtL(f̆1-Dk+1/2)f̆1-Dk+1/2
end while
f̆3-Df̆1-D
f=W3D1{f̆3-D}

Algorithm 2

Fully compressed matrix.

Initialize an empty matrix A of dimensions (Ms×Mψ)×(Nx×Ny×Nz)
Initialize an empty matrix J¯ of dimensions (Ms×Mψ)×NΨ
Initialize an empty vector f¯s of size (Nx×Ny×Nz)
fori=1,,(Ms×Mψ)do
ri=M[ΦeΦf*]/M[ΦeΦe*]
R3-D=W3-D{ri}
R1-DR3-D
AiR1-D
end for
Rssortmean(Ai) in descending order and store locations jN corresponding to the NΨth largest elements
J¯ ← elements jN of A.
Solve the inverse problem for f¯1-D:
Set tolerance value ε and maximum iteration number Nk
whiley˜J¯f¯1-D2<εk<Nkdo
f¯1-Dk+1/2=f¯1-Dk+J¯T(J¯TJ¯+λI)1(y˜J¯f¯1-Dk)
f¯1-Dk+1/2=f¯1-Dk+1/2+ΔtL(f¯1-Dk+1/2)f¯1-Dk+1/2
end while
f¯s ← set elements jN of f¯s to their corresponding f¯1-D values.
f¯3-Df¯s
f=W3-D1{f¯3-D}

2.5.

Evaluation

The performance of our data and solution compression method is evaluated through simulations and a phantom study. The matrix J is calculated using TOAST, which is an FEM based software developed at University College London (UCL). Images are reconstructed (a) using zero-order Tikhonov regularization and solving the linear problem directly and (b) iteratively using our ADSP method.20 Solutions are computed on a tetrahedral mesh and mapped into a regular 64×64×64 grid covering the domain of the mesh.

In the following sections, we refer to the data compression method as DC, sparse matrix method as Sx and the fully compressed matrix method as FCx, where x represents the number of coefficients NΨ.

Wavelet filters: After experimenting different wavelet filters, we choose to use the Battle–Lemarié wavelet filters, which appear to be suitable for this type of application. The Battle–Lemarié wavelets are symmetric, orthonormal, and smooth. These wavelets are attractive since they have good time and frequency localization properties and can approximate smooth solutions.22 Other popular wavelets, such Daubechies, are not smooth enough and do not possess good frequency localization.23

Battle–Lemarié wavelet filters are built from polynomial splines of order 2p+1. Let ϕ(ω) be the Fourier transform of the wavelet scaling function ϕ(t). The low-pass quadrature mirror filter H(ω) can be obtained from the following relation in the Fourier domain22

Eq. (19)

ϕ(2ω)=H(ω)ϕ(ω),whereH(ω)n=h(n)einω.

Lemarié has shown that the scaling function can be written as

Eq. (20)

ϕ(ω)=1ωnΣ2n(ω),wheren=2+2p,
and where

Eq. (21)

Σn(ω)=k=1(ω+2kπ)n.
Alternatively, one can calculate Σn(ω) by computing the n2 derivative of the formula:

Eq. (22)

Σ2(ω)=14sin2(ω/2).

For an approximation built from cubic splines, i.e., p=1, and thus, n=4, it follows from the previous equations that

Eq. (23)

ϕ(ω)=1ω4Σ8(ω)andH(ω)=Σ8(ω)28Σ8(2ω).

Figures of merit: The relative error (RE) and contrast-to-noise ratio (CNR) are used to evaluate the quality of the reconstructed images. The RE between the true image ftrue and the reconstructed image frecon is defined as RE=ftruefrecon2/ftrue2. The CNR is calculated as follows:24

Eq. (24)

CNR=μROIμB[wROIσROI2wBσB2]1/2,
where μROI is the mean value of the region of interest (ROI); μB is the mean value of the background (B); σROI and σB are the standard deviations of the ROI and background, respectively; and wROI=AROI/(AROI+AB) and wB=AB/(AROI+AreaB) are noise weights, where AROI and AB are the ROI and B areas, respectively. The ROI is defined as the region where voxels have attained at least half maximum intensity and the remaining voxels represent the background B.

2.5.1.

Simulation

The geometry used in the simulation is similar to the phantom geometry. The matlab package Iso2mesh25 was used to generate a tetrahedral mesh from the X-ray computed tomography (XCT) image of a cylindrical phantom (see Sec. 2.5.2 and Fig. 3). The FEM mesh has 8449 nodes, 49,797 elements and dimensions 28×28×42mm3. Similar to the phantom, the optical properties are μa=0.01mm1 and μs=0.8mm1, and the fluorescent targets are two vertical cylinders placed parallel to each other (Fig. 3) with contrast f=2. Projection images of size 128×128 pixels are calculated for 18 evenly spaced source-camera positions over a full 360 degrees range at z=0mm (Fig. 3). Data consist of fluorescence and excitation projections with 2% Gaussian random noise. The data space is compressed by setting Mψ=32 and solutions are computed with NΨ={8192,4096,1024,256}. For the reconstructions using zero-order Tikhonov regularization the damping factor is initially set to λTik=atrace(JcJcT) and for reconstructions using ADSP regularization it is set to λADSP=10atrace(JcJcT), where a=3×104 and Jc are the different types of compressed Jacobian matrices.

Fig. 3

Geometry used in the simulations. The FEM mesh is generated from X-ray CT images of the experimental phantom. The asterisk markers indicate the source positions.

JBO_18_8_086008_f003.png

2.5.2.

Phantom

Data are acquired using the rotating fDOT–XCT system described in Ref. 2. The optical and XCT images are acquired sequentially. The solid phantom consists of a mixture of Agar, Intralipid, and ink. The stock solution is 4.56 ml Intralipid 20% (Sigma–Aldrich Co. LLC, St. Louis, USA), 0.256 ml India ink, 100 ml deionized water and 2.8 g Agar (Sigma–Aldrich Co. LLC, St. Louis, USA). The phantom is a cylinder with diameter =28mm, height L=65mm and optical properties μa=0.01mm1 and μs=0.8mm1, which were measured using a spectrometer (USB4000-VIS-NIR, Ocean Optics, London, UK). Two translucent tubes with an inside diameter of 3 mm (and outside diameter 4 mm) are inserted parallel to its long axis and filled with an optically matched fluid containing a fluorescent dye (Alexa 750) at 500 nM (molar) concentration. Excitation and fluorescence data images of size 512×512 (resized to 128×128) are recorded at 18 angular positions (same than the simulations). The FEM mesh is the same than the one used in the simulation study. We use the compression Mψ=32 and NΨ={1024,256} (NΨ=256 when ADSP regularization is used). The damping factor is initially set to λTik=atrace(JcJcT) and λADSP=10atrace(JcJcT), where a=8×103.

3.

Results

Profile plots across the reconstructed fluorescence targets (z=0mm and y=0mm) are shown in Figs. 4 and 5 for images obtained using Tikhonov regularization for the sparse and fully compressed matrix cases, respectively. Similarly, Figs. 6 and 7 show the profile plots for images obtained using ADSP regularization. Figures 4 and 6 also show the profile plots for the data compression case.

Fig. 4

Profile plots across the fluorescence images (z=0mm and y=0mm) reconstructed from simulated data using Tikhonov regularization and data compression, a sparse matrix with different compression levels, and the true solution.

JBO_18_8_086008_f004.png

Fig. 5

Profile plots across the fluorescence images (z=0mm and y=0mm) reconstructed from simulated data using Tikhonov regularization and a fully compressed matrix with different compression levels, and the true solution.

JBO_18_8_086008_f005.png

Fig. 6

Profile plots across the fluorescence images (z=0mm and y=0mm) reconstructed from simulated data using ADSP regularization and data compression, a sparse matrix with different compression levels, and the true solution.

JBO_18_8_086008_f006.png

Fig. 7

Profile plots across the fluorescence images (z=0mm and y=0mm) reconstructed from simulated data using ADSP regularization and a fully compressed matrix with different compression levels, and the true solution.

JBO_18_8_086008_f007.png

Since the profiles across the fluorescent inclusions do not differ significantly for the different types of compression, particularly when ADSP regularization is used, only results obtained without solution compression and high levels of solution compressions are shown. Therefore, reconstructions obtained using ADSP are shown for data compression and solution compression with NΨ=256. Reconstructions obtained using Tikhonov regularization and solution compression with NΨ={1024,256} are displayed for comparison with the case where only data compression is used, to show the effect of the solution compression on the noise in the reconstructions.

The fluorescence distribution images reconstructed from simulated data using zero-order Tikhonov regularization are shown in Figs. 8Fig. 9Fig. 10Fig. 1112. The reconstructions are displayed as 3-D isosurfaces (70% of the maximum value) together with the structural volumetric image. Additionally, a 2-D cross section of the reconstructed volume is taken at z=0mm and thresholded at 70% of the maximum value. The thresholded images are overlaid onto the structural image. Figures 13Fig. 1415 show the images reconstructed using ADSP regularization and thresholded at 50% of the maximum value. The threshold value is lower since these images are not as contaminated with noise as the previous ones.

Fig. 8

Fluorescence yield f reconstructed from simulated data using Tikhonov regularization and data compression in (left) 3-D and (right) 2-D cross section at z=0mm.

JBO_18_8_086008_f008.png

Fig. 9

Fluorescence yield f reconstructed from simulated data using Tikhonov regularization and S1024 in (left) 3-D and (right) 2-D cross section at z=0mm.

JBO_18_8_086008_f009.png

Fig. 10

Fluorescence yield f reconstructed from simulated data using Tikhonov regularization and S256 in (left) 3-D and (right) 2-D cross section at z=0mm.

JBO_18_8_086008_f010.png

Fig. 11

Fluorescence yield f reconstructed from simulated data using Tikhonov regularization and FC1024 in (left) 3-D and (right) 2-D cross section at z=0mm.

JBO_18_8_086008_f011.png

Fig. 12

Fluorescence yield f reconstructed from simulated data using Tikhonov regularization and FC256 in (left) 3-D and (right) 2-D cross section at z=0mm.

JBO_18_8_086008_f012.png

Fig. 13

Fluorescence yield f reconstructed from simulated data using ADSP regularization and data compression in (left) 3-D and (right) 2-D cross section at z=0mm.

JBO_18_8_086008_f013.png

Fig. 14

Fluorescence yield f reconstructed from simulated data using ADSP regularization and S256 in (left) 3-D and (right) 2-D cross section at z=0mm.

JBO_18_8_086008_f014.png

Fig. 15

Fluorescence yield f reconstructed from simulated data using ADSP regularization and FC256 (left) 3-D and (right) 2-D cross section at z=0mm.

JBO_18_8_086008_f015.png

Similarly, the images reconstructed from phantom data using zero-order Tikhonov regularization are shown in Figs. 16Fig. 17Fig. 18Fig. 1920 and using ADSP regularization are shown in Figs. 21Fig. 2223. Their respective profile plots are shown in Figs. 24 and 25.

Fig. 16

Fluorophore concentration (nM) reconstructed from phantom data using Tikhonov regularization and data compression in (left) 3-D and (right) 2-D cross section at z=0mm.

JBO_18_8_086008_f016.png

Fig. 17

Fluorophore concentration (nM) reconstructed from phantom data using Tikhonov regularization and S1024 in (left) 3-D and (right) 2-D cross section at z=0mm.

JBO_18_8_086008_f017.png

Fig. 18

Fluorophore concentration (nM) reconstructed from phantom data using Tikhonov regularization and S256 in (left) 3-D and (right) 2-D cross section at z=0mm.

JBO_18_8_086008_f018.png

Fig. 19

Fluorophore concentration (nM) reconstructed from phantom data using Tikhonov regularization and FC1024 in (left) 3-D and (right) 2-D cross section at z=0mm.

JBO_18_8_086008_f019.png

Fig. 20

Fluorophore concentration (nM) reconstructed from phantom data using Tikhonov regularization and FC256 in (left) 3-D and (right) 2-D cross section at z=0mm.

JBO_18_8_086008_f020.png

Fig. 21

Fluorophore concentration (nM) reconstructed from phantom data using ADSP regularization and data compression in (left) 3-D and (right) 2-D cross section at z=0mm.

JBO_18_8_086008_f021.png

Fig. 22

Fluorophore concentration (nM) reconstructed from phantom data using ADSP regularization and S256 in (left) 3-D and (right) 2-D cross section at z=0mm.

JBO_18_8_086008_f022.png

Fig. 23

Fluorophore concentration (nM) reconstructed from phantom data using ADSP regularization and FC256 in (left) 3-D and (right) 2-D cross section at z=0mm.

JBO_18_8_086008_f023.png

Fig. 24

Profile plots across the fluorescence images (z=0mm and y=0mm) reconstructed from phantom data using Tikhonov regularization and DC, S256, FC256, and the true solution.

JBO_18_8_086008_f024.png

Fig. 25

Profile plots across the fluorescence images (z=0mm and y=0mm) reconstructed from phantom data using ADSP regularization and DC, S256, FC256, and the true solution.

JBO_18_8_086008_f025.png

Figures 26 and 27 show the RE and CNR, respectively, of the images reconstructed from simulated and phantom data using Tikhonov and ADSP regularization for the different compression methods. Figure 28 shows the condition number of the matrices in the data compression, compressed sparse, and fully compressed format.

Fig. 26

RE of images reconstructed from (left column) simulated data and (right column) phantom data using (top row) Tikhonov regularization and (bottom row) ADSP regularization, for the different compression cases.

JBO_18_8_086008_f026.png

Fig. 27

CNR of images reconstructed from (left column) simulated data and (right column) phantom data using (top row) Tikhonov regularization and (bottom row) ADSP regularization, for the different compression cases.

JBO_18_8_086008_f027.png

Fig. 28

Condition number of the data compressed, sparse, and fully compressed matrices.

JBO_18_8_086008_f028.png

The reconstructions are faster when a fully compressed matrix is used, and the sparse matrix method is faster than using data compression only, with a ratio of 9.3:3.9:1.

4.

Discussion and Conclusion

In fDOT, the data and solution spaces can be quite large, and hence, the sensitivity matrix may be too large to store. Even if this matrix can be stored, it may not be possible to directly invert it, and thus, iterative methods must be employed. However, these are computationally expensive and may take a long time to converge. The method we proposed reduces the dimensions of the sensitivity matrix allowing the storage and direct inversion of the matrix and, for the sparse matrix case, without requiring the construction of the full matrix. Both sparse and fully compressed matrix representations allow storage of large amounts of data generated from CCD-based fDOT systems, rather than keeping redundant data and solution components.

The performance of our method was tested on simulated and experimental phantom data. Our simulations show that our solution compression approach can help to suppress noise in the reconstructed images (Figs. 8Fig. 9Fig. 10Fig. 1112). This result is not surprising, since wavelet transforms are a common and efficient technique for denoising. The use of structural a priori information in the image reconstruction greatly improved the quality of the images (Figs. 6, 7, and 13Fig. 1415). The images reconstructed using the different types of compression are qualitatively very similar. However, the fluorescence yield was underestimated when both data and solution compression were used with ADSP regularization. The images reconstructed from experimental data are similar (Figs. 16Fig. 17Fig. 18Fig. 1920). The target locations and dimensions were accurately estimated by all the methods when structural information was incorporated into the image reconstruction (Figs. 21Fig. 2223). The profiles plots in Figs. 24 and 25 also show that the use of structural information can greatly improve the quality of the reconstructions.

Our results show that the sensitivity matrix can be effectively compressed into a wavelet space by selecting the most significant wavelet coefficients. The sensitivity matrix, which in our studies has dimensions 576×262,144, can effectively be represented in a fully compressed form with dimension 576×256 or sparse form with size 576×262,144 and 576×256 nonzero elements.

Figure 26 shows that the RE decreases with increasing level of compression, i.e., the quality of the reconstruction improves, which is presumably a consequence of the decrease of the condition number (Fig. 28). The fully compressed matrix with NΨ=1024 has a higher condition number, and hence, one should expect a higher RE. Nevertheless, the denoising characteristic of the wavelets seems to compensate for this effect. Figure 27 shows that the CNR of the images reconstructed from simulated data improves when solution compression is used, whereas the CNR of the images obtained from phantom data decreases. Overall, the simulation results using Tikhonov regularization with data and solution compression are quite similar, but the lowest RE was achieved with the FC256 method and the highest CNR using FC1024. For the phantom study, the best results were obtained with the FC256. The images reconstructed using ADSP and S256 had the lowest RE and highest CNR.

The sparse matrix compression method has the advantage that the problem becomes less ill-conditioned and it does not require the explicit construction of the full matrix. The fully compressed matrix method provides images that are visually quite similar, but it underestimates the fluorescence yield more than the latter when ADSP is included. This method requires the construction of the full matrix and the wavelet coefficients kept in the compressed form are all in the same jN locations in all the rows of the matrix. These are not necessarily the most relevant coefficients, which may explain the lower qualitative accuracy. Nevertheless, this method allows a faster inversion process and, when Tikhonov regularization is used, provides images with the lowest RE and highest CNR.

The disadvantage of using ADSP regularization is that the solution depends on the parameter λ value and the number of iterations until convergence is reached, i.e., it depends on the convergence criteria, and thus, there is still potential for improvement. Nevertheless, the results obtained using data and sparse matrix compression are quite satisfactory.

In this work, we map from the FEM mesh to a regular grid, which is then transformed into the wavelet space, with no regard to the irregular boundary of the object. In future work, we will attempt to apply a similar compression directly on the FEM mesh, so that the reconstruction is performed on the mesh and to avoid error propagation resulting from these mappings. Furthermore, the boundary of the object introduces sharp edges in the images, which may lead to artifacts in the compressed representation.

There is a large diversity of available wavelet filters, and hence, we performed a simple qualitative comparison of the suitability of some of the most popular ones (e.g., Haar, Daubechies, Coiflet, Symmlets, Battle–Lemarié, etc.) for use in fDOT. However, improved results can be expected if wavelet filters more adequate for the fDOT problem are designed and the ideal number of wavelet coefficients to keep is identified.

In summary, image reconstruction using data and solution compression can greatly reduce the dimensions and ill-conditioning of the problem, allowing fast and less noisy reconstructions.

Acknowledgments

This work has also been supported by funding from ECs seventh framework programme FMT-XCT under Grant Agreement No. 201792.

References

1. 

V. Ntziachristos, “Fluorescence molecular imaging,” Annu. Rev. Biomed. Eng., 8 1 –33 (2006). https://doi.org/http://dx.doi.org/10.1146/annurev.bioeng.8.061505.095831 ARBEF7 1523-9829 Google Scholar

2. 

R. B. Schulzet al., “Hybrid system for simultaneous fluorescence and X-ray computed tomography,” IEEE Trans. Med. Imag., 29 (2), 465 –473 (2010). https://doi.org/http://dx.doi.org/10.1109/TMI.2009.2035310 ITMID4 0278-0062 Google Scholar

3. 

V. MarkelJ. Schotland, “Symmetries, inversion formulas and image reconstruction for optical tomography,” Phys. Rev. E, 70 056616 –056619 (2004). https://doi.org/http://dx.doi.org/10.1103/PhysRevE.70.056616 PLEEE8 1063-651X Google Scholar

4. 

Z. Wanget al., “Experimental demonstration of an analytic method for image reconstruction in optical tomography with large data sets,” Optic. Lett., 30 (24), 3338 –3340 (2005). Google Scholar

5. 

V. MarkelJ. Schotland, “The inverse problem in optical diffusion tomography: III. Inversion formulas and singular value decomposition,” J. Opt. Soc. A, 20 (5), 890 –902 (2003). https://doi.org/http://dx.doi.org/10.1364/JOSAA.20.000890 JOSAAH 0030-3941 Google Scholar

6. 

S. D. Koneckyet al., “Imaging complex structures with diffuse light,” Opt. Express, 16 (7), 5048 –5060 (2008). https://doi.org/http://dx.doi.org/10.1364/OE.16.005048 OPEXFF 1094-4087 Google Scholar

7. 

M. XuM. LaxR. Alfano, “Time-resolved Fourier optical diffuse tomography,” J. Opt. Soc. Am. A, 18 (7), 1535 –1542 (2001). https://doi.org/http://dx.doi.org/10.1364/JOSAA.18.001535 JOAOD6 0740-3232 Google Scholar

8. 

V. LukicV. MarkelJ. Schotland, “Optical tomography with structured illumination,” Opt. Lett., 34 (7), 983 –985 (2009). https://doi.org/http://dx.doi.org/10.1364/OL.34.000983 OPLEDP 0146-9592 Google Scholar

9. 

N. Ducroset al., “Full-wavelet approach for fluorescence diffuse optical tomography with structured illumination,” Opt. Lett., 35 (21), 3676 –3678 (2010). https://doi.org/http://dx.doi.org/10.1364/OL.35.003676 OPLEDP 0146-9592 Google Scholar

10. 

J. Ripoll, “Hybrid Fourier-real space method for diffuse optical tomography,” Opt. Lett., 35 (5), 688 –690 (2010). https://doi.org/http://dx.doi.org/10.1364/OL.35.000688 OPLEDP 0146-9592 Google Scholar

11. 

A. Zacharopouloset al., “A matrix-free algorithm for multiple wavelength fluorescence tomography,” Opt. Express, 17 (5), 3025 –3035 (2009). https://doi.org/http://dx.doi.org/10.1364/OE.17.003042 OPEXFF 1094-4087 Google Scholar

12. 

W. Zhuet al., “A wavelet-based multiresolution regularized least squares reconstruction approach for optical tomography,” IEEE Trans. Med. Imag., 16 (2), 210 –217 (1997). https://doi.org/http://dx.doi.org/10.1109/42.563666 ITMID4 0278-0062 Google Scholar

13. 

A. Landragin-Frassatiet al., “Application of a wavelet-Galerkin method to the forward problem resolution in fluorescence diffuse optical tomography,” Opt. Express, 17 (21), 18433 –18448 (2009). https://doi.org/http://dx.doi.org/10.1364/OE.17.018433 OPEXFF 1094-4087 Google Scholar

14. 

J. BaccouJ. Liandrat, “On coupling wavelets with fictitious domain approaches,” Appl. Math. Lett., 18 (12), 1325 –1331 (2005). https://doi.org/http://dx.doi.org/10.1016/j.aml.2005.02.028 AMLEEL 0893-9659 Google Scholar

15. 

A. Landragin-Frassatiet al., “Model reduction using wavelet multiresolution technique applied to fluorescence diffuse optical tomography,” Appl. Opt., 48 (36), 6878 –6892 (2009). https://doi.org/http://dx.doi.org/10.1364/AO.48.006878 APOPAI 0003-6935 Google Scholar

16. 

T. J. RudgeV. Y. SolovievS. R. Arridge, “Fast image reconstruction in fluorescence optical tomography using data compression,” Opt. Lett., 35 (5), 763 –765 (2010). https://doi.org/http://dx.doi.org/10.1364/OL.35.000763 OPLEDP 0146-9592 Google Scholar

17. 

S. R. ArridgeJ. C. Shotland, “Optical tomography: forward and inverse problems,” Inverse Probl., 25 (12), 123010 (2009). https://doi.org/http://dx.doi.org/10.1088/0266-5611/25/12/123010 INPEEY 0266-5611 Google Scholar

18. 

A. SoubretJ. RipollV. Ntziachristos, “Accuracy of fluorescent tomography in the presence of heterogeneities: study of the normalized Born ratio,” IEEE Trans. Med. Imag., 24 (10), 1377 –1386 (2005). https://doi.org/http://dx.doi.org/10.1109/TMI.2005.857213 ITMID4 0278-0062 Google Scholar

19. 

S. R. Arridge, “Optical tomography in medical imaging,” Inverse Probl., 15 (2), R41 –R93 (1999). https://doi.org/http://dx.doi.org/10.1088/0266-5611/15/2/022 INPEEY 0266-5611 Google Scholar

20. 

T. Correiaet al., “Split operator method for fluorescence diffuse optical tomography using anisotropic diffusion regularisation with prior anatomical information,” Biomed. Opt. Exp., 2 (9), 2632 –2648 (2011). https://doi.org/http://dx.doi.org/10.1364/BOE.2.002632 BOEICL 2156-7085 Google Scholar

21. 

P. PeronaJ. Malik, “Scale-space and edge detection using anisotropic diffusion,” IEEE Trans. Pattern Anal. Mach. Intell., 12 (7), 629 –639 (1990). https://doi.org/http://dx.doi.org/10.1109/34.56205 ITPIDJ 0162-8828 Google Scholar

22. 

S. Mallat, “A theory for multiresolution signal decomposition: the wavelet representation,” IEEE Trans. Pattern Anal. Mach. Intell., 11 (7), 674 –693 (1989). https://doi.org/http://dx.doi.org/10.1109/34.192463 ITPIDJ 0162-8828 Google Scholar

23. 

I. Daubechies, Ten Lectures on Wavelets, SIAM, Philadelphia (1992). Google Scholar

24. 

X. Songet al., “Automated region detection based on the contrast to noise ratio in near-infrared tomography,” Appl. Opt., 43 (5), 1053 –1062 (2004). https://doi.org/http://dx.doi.org/10.1364/AO.43.001053 APOPAI 0003-6935 Google Scholar

25. 

Q. FangD. Boas, “Tetrahedral mesh generation from volumetric binary and gray-scale images,” in Proc. IEEE Int. Symposium on Biomedical Imaging, 1142 –1145 (2009). Google Scholar
© 2013 Society of Photo-Optical Instrumentation Engineers (SPIE) 0091-3286/2013/$25.00 © 2013 SPIE
"Wavelet-based data and solution compression for efficient image reconstruction in fluorescence diffuse optical tomography," Journal of Biomedical Optics 18(8), 086008 (12 August 2013). https://doi.org/10.1117/1.JBO.18.8.086008
Published: 12 August 2013
JOURNAL ARTICLE
14 PAGES


SHARE
Advertisement
Advertisement
Back to Top