## 1.

## Introduction

In optical tomographic (OT) imaging, one attempts to determine the spacial distribution of optical properties inside a biological tissue, based on measurements taken on the boundary of the tissue. This type of inverse problem is in general highly ill-posed. This means that many different optical property distributions inside the medium can lead to the same light distribution on the surface of the medium. One of the reasons for this problem is that the number of unknown parameters to be reconstructed often by far exceeds the number of boundary measurements available. Therefore, the number of known data points is substantially smaller than the number of unknown parameters in the resulting system of algebraic equations. To overcome this problem, one can try to increase the number of data points. To this end, researchers have suggested, for example, multispectral methods^{1} as well as multifrequency measurement systems.^{2} In addition, or alternatively, one can seek to reduce the number of unknowns. In this direction, various groups have suggested combining optical tomography (OT) with other imaging modalities, such as magnetic resonance imaging (MRI)^{3} or ultrasound (US).^{4} In these cases, the other imaging modalities provide the prior structural information for optical imaging. Furthermore, Schweiger and Arridge^{5} have tested different local basis functions to convert spatial optical properties to certain local basis coefficients. In this way, an image is described by a “small” number of coefficients, rather than by each pixel value in the image.

In this paper, we propose a substantially different approach. We hypothesize that the number of unknowns can be drastically reduced by parameterizing the spatial distribution of optical properties via a globe basis transform and then reconstruct the coefficients of this global basis transform. As the global basis function, we employ in this work the DCT, which has long been used in the field of image compression.^{6, 7} It has been shown that using the DCT, one can achieve high compression ratios for smooth images without losing much information. Optical tomographic images are typically very smooth due to the strong scattering of light in most biological tissues. It can be expected that only a few DCT coefficients are needed to recover the spatial variations seen in these images. By reconstructing the DCT coefficients, rather than each pixel value of an image, the number of unknown parameters can be greatly reduced as compared to conventional optical tomographic image reconstruction schemes.

To illustrate the performance of the parametric-DCT reconstruction technique, we employ in this paper the frequency-domain equation of radiative transport (FD-ERT). In Sec. 2, the pertinent aspects of the FD-ERT and related conventional reconstruction schemes are reviewed. Section 3 is devoted to the DCT and its application to OT. In Sec. 4, two numerical reconstruction examples, one containing a capsular perturbation and the other with the uniform distribution in one direction, are presented. The influence of varying initial guesses and noise levels on parametric reconstruction will be discussed. In Sec. 5, we present reconstruction results obtained with experimental data from a tissue-like phantom. Throughout the paper, we compare the performance of the novel parametric-DCT method to more conventional approaches to the reconstruction problem in OT.

## 2.

## Frequency Domain Optical Tomography with the ERT

We implemented the parametric-DCT reconstruction method by adapting an existing model-based iterative image reconstruction (MOBIIR) algorithm that uses the frequency domain ERT.^{8} In the following, we briefly review some pertinent details of this MOBIIR algorithm.

## 2.1.

### The Forward Problem

A frequency-domain ERT describes the photon density in phase space, i.e., as a function of position $x\u220a\mathcal{D}\subset {\mathbb{R}}^{3}$ and direction $\theta \u220a{S}^{2}$ (unit sphere of ${\mathbb{R}}^{3}$ ), and can be expressed as:

## 1.

## Eq. 2

${\Gamma}_{\pm}\u2254\{(\mathbf{x},\mathbf{\theta})\u220a\partial \mathcal{D}\times {S}^{2}\phantom{\rule{0.3em}{0ex}}\mathrm{s.t.}\phantom{\rule{0.3em}{0ex}}\pm \mathbf{\theta}\cdot \mathbf{\nu}\left(\mathbf{x}\right)>0\}.$^{9}

## Eq. 3

$k(\mathbf{\theta}\cdot {\mathbf{\theta}}^{\prime})=\frac{1-{g}^{2}}{{(1+{g}^{2}-2g\phantom{\rule{0.2em}{0ex}}\mathrm{cos}\phantom{\rule{0.2em}{0ex}}\varphi )}^{3\u22152}},$Solving the forward problem 1 yields the quantity of photon radiance. However, in OT experiments, the quantity that one measures is outgoing photon current. This current can be expressed as

## Eq. 4

$J\left({\mathbf{x}}_{d}\right)={\int}_{{\Gamma}_{+}}\mathbf{\theta}\cdot \mathbf{\nu}\left({\mathbf{x}}_{d}\right)u({\mathbf{x}}_{d},\mathbf{\theta})\mathrm{d}\mathbf{\theta},$With the combination of the discrete ordinates method^{10} for the angular variable, a finite-volume discretization method for the space variable,^{11} and an appropriate boundary condition, the continuous ERT can be converted to the following linear algebraic equation:

## 2.2.

### The Inverse Problem and Conventional MOBIIR Methods

After discretizing Eq. 1, one can predict the radiance field values of all the discretized elements and the photon density current values of the boundary elements. However, in real OT experiments, only a limited number of boundary measurements are available, which is typically far less than the number of the discretized mesh cells. Thus, only an underdetermined system is available to be used to reconstruct optical property maps. Classically, a least-squares problem is formulated and a regularization term is added to compensate for this underdetermination. Numerically, this requires the minimization of the following objective function:

## Eq. 6

$\mathrm{min}\phantom{\rule{0.2em}{0ex}}{\Phi}_{\alpha}({\mu}_{a},{\mu}_{s})=\frac{1}{2}\sum _{i=1}^{{N}_{s}}\sum _{j=d}^{{N}_{d}}{({\mathcal{P}}_{i,d}{\mathbf{u}}_{i}-{z}_{i,d})}^{2}+\frac{\alpha}{2}\sum _{\mathcal{C}=1}^{N}\sum _{p=a,s}[\sum _{\kappa =x,y,z}{\left({\mathcal{D}}_{\kappa}^{\mathcal{C}}{\mu}_{p}\right)}^{2}+{\left({\mu}_{p}^{\mathcal{C}}\right)}^{2}],$Summarizing the preceding solvers, we plot the flowchart of a typical MOBIIR algorithm in Fig. 1 . The code starts with an initial guess of the distribution of optical properties inside the medium. Using the forward solver, the predictions of the outgoing current on the boundary of the medium are calculated. The inverse solver compares these predicted detector readings with actual measurement data using the objective function. Subsequently, the code calculates an update of the optical property distribution using a quasi-Newton scheme. The update is used as an input for the forward solver in the next iteration. The updating will continue until the value of the objective function falls below a preset error threshold $\u03f5$ .

## 3.

## Parametric Reconstruction with the DCT

## 3.1.

### Image Compression with DCT

Various compression algorithms have been developed to reduce the memory needed to store images.^{13, 14, 15} For example, the DCT is at the core of JPEG image compression.^{7} The general formula for a two-dimensional (2-D) DCT of an image with
${M}_{x}\times {M}_{y}$
pixels is given by:

## 7.

## 8.

If we attempt to reconstruct an image that contains a sharp-edge target, more DCT coefficients have to be used. This is illustrated in Fig. 3a, where an image contains a square with sharp edges. The DCT coefficient map of this image is plotted in Fig. 3b. A reconstructed image using $6\times 6$ DCT coefficients (one-seventh of total number of DCT coefficients) is presented in Fig. 3c. Figure 3d shows its corresponding relative error map. Notice that the error map of the reconstructed image reaches 20%, which is substantially larger than the one shown in Fig. 2d. In order to decrease the error down to 0.2%, we have to use $39\times 39$ DCT coefficients. Figures 3e and 3f show the IDCT recovered image and the image errors, respectively.

## 3.2.

### OT with Parametric-DCT Method

To accurately calculate the photon propagation requires a finely discretized mesh. We found that in some cases, up to 10 discretized elements per photon mean free path are needed. Such a finely discretized mesh induces images that contain a very large number of pixels. The algorithm described in Sec. 2.2 treats each discretized element of an image as an independent variable. Therefore, the forward discretization results in a large number of unknowns for the inverse computation.

The image compression technique with DCT can help to reduce the number of these unknowns. We write the spatial distributed optical properties in a three-dimensional (3-D) domain as:

## Eq. 9

${\mu}_{a\left(s\right)}(x,y,z)\simeq \sum _{{k}_{x}=0}^{{K}_{x}}\sum _{{k}_{y}=0}^{{K}_{y}}\sum _{{k}_{z}=0}^{{K}_{z}}{A}_{a\left(s\right)}({k}_{x},{k}_{y},{k}_{z})\mathrm{cos}\left(\frac{2x+\Delta x}{2{L}_{x}}{k}_{x}\pi \right)\mathrm{cos}\left(\frac{2y+\Delta y}{2{L}_{y}}{k}_{y}\pi \right)\mathrm{cos}\left(\frac{2z+\Delta z}{2{L}_{z}}{k}_{z}\pi \right),$With Eq. 9, the objective function that we want to minimize takes on the following form:

## Eq. 10

$\mathrm{min}\phantom{\rule{0.2em}{0ex}}\Phi ({A}_{a},{A}_{s})=\frac{1}{2}\sum _{i=1}^{{N}_{s}}\sum _{d=1}^{{N}_{d}}{({\mathcal{P}}_{i,d}{\mathbf{u}}_{i}-{z}_{i,d})}^{2},$^{16}tend to smooth images by removing high-spatial-frequency components. Using a parametric reconstruction technique based on the DCT expansion, the higher frequency components are automatically filtered out and hence the use of an extra regularization term is not needed.

To minimize the objective function in Eq. 10, we employ a quasi-Newton approach with BFGS updates of the Hessian matrix.^{17} The gradient of the objective function (10) can be written as:

## Eq. 11

$\frac{\partial \Phi}{\partial {A}_{a\left(s\right)}({k}_{x},{k}_{y},{k}_{z})}=\sum ^{N}\frac{\partial \Phi}{\partial {\mu}_{a\left(s\right)}(x,y,z)}\phantom{\rule{0.2em}{0ex}}\mathrm{cos}\left(\frac{2x+\Delta x}{2{L}_{x}}{k}_{x}\pi \right)\mathrm{cos}\left(\frac{2y+\Delta y}{2{L}_{y}}{k}_{y}\pi \right)\mathrm{cos}\left(\frac{2z+\Delta z}{2{L}_{z}}{k}_{z}\pi \right),$The overall flowchart for an image reconstruction algorithm that makes use of the parametric-DCT method is presented in Fig. 4 . By comparing this flow chart with the flow chart for the conventional MOBIIR approach shown in Fig. 1, one can easily see the main differences between these two methods. Both methods employ the same forward solver and use a quasi-Newton method to update their system parameters. However, in the DCT approach, the gradient of the objective function is calculate with respect to the DCT coefficients, while in the conventional MOBIIR scheme the gradient is calculated with respect to the optical properties. Subsequently, the updated DCT coefficients are used to update the optical properties using Eq. 9; and the new optical property values become input to the forward solver in the next iteration step.

## 3.3.

### Number of DCT Coefficients

Choosing an appropriate number of DCT coefficients for the reconstruction is an important issue. The number of necessary DCT coefficients used to store and recover an image is determined by the smoothness of the image. The smoother an image is, the fewer DCT coefficients are required. As shown in Fig. 2 and Fig. 3, a sharp-edged object in an image requires more DCT coefficients. If we have some prior knowledge of the image features, we can choose a suitable number of DCT coefficients to represent the image and later use that number of DCT coefficients to perform the reconstruction. Such prior knowledge will help to improve the quality of the reconstructed images (see Sec. 4.3).

Some information about the number of DCT coefficients necessary for optical tomographic imaging can be gleaned from studies of resolution limits and minimal detectable object size in OT. These topics have been studied by several groups.^{3, 18, 19, 20, 21} These authors found that the achievable resolution depends on tissues’ size, geometry, optical properties, imaging modality and many other parameters. Pogue
^{22} argued that standard resolution testing is optimal when infinite contrast is used and hardware evaluation is the goal. However, deep tissue imaging of absorption or fluorescent contrast agents *in vivo* requires a more detailed analysis that takes tissue contrast into account. Overall, they found that for most biological tissues, a good approximation is to assume that the minimum detectable size is in the range of
$1\u221510$
of the outer diameter of the object imaged. For our analysis, this finding can be interpreted as to mean that smallest spatial variations in reconstructed optical tomographic images are in the range of
$1\u221510$
of the outer diameter.

Given this degree of spatial variation, we can estimate an approximate number of necessary DCT coefficients as follows. We define ${K}_{x}$ (or $y,z$ ) as the maximum number of DCT coefficients needed in one dimension, $\Delta x$ (or $y,z$ ) as the resolution of an OT image along $x$ (or $y,z$ ) direction, and ${L}_{x}$ (or $y,z$ ) as the length of an OT image in $x$ (or $y,z$ ) direction. Then we can write:

## Eq. 12

${K}_{x(\text{or}\phantom{\rule{0.3em}{0ex}}y,z)}=\frac{{L}_{x(\text{or}\phantom{\rule{0.3em}{0ex}}y,z)}}{\lambda \u22152}=\frac{{L}_{x(\text{or}\phantom{\rule{0.3em}{0ex}}y,z)}}{\Delta x(\text{or}\phantom{\rule{0.3em}{0ex}}\Delta y,\phantom{\rule{0.3em}{0ex}}\Delta z)}.$^{22}we see that one needs approximately 10 coefficients for each dimensions. We used these considerations as a guide throughout our studies.

## 4.

## Numerical Examples

## 4.1.

### Tissue Geometry

We illustrate the performance of the parametric-DCT reconstruction algorithm by using three-dimensional (3-D) numerical examples. The geometry of the problems considered in this work are shown in Fig. 5 . First, we consider a cylindrical computational domain defined as $D\u2254\left\{{(x,y,z)}^{T}\right|{[{(x-2.05)}^{2}+{(y-2.05)}^{2}]}^{1\u22152}\u2a7d2.05\phantom{\rule{0.3em}{0ex}}\mathrm{cm};0.0\u2a7dz\u2a7d3.0\phantom{\rule{0.3em}{0ex}}\mathrm{cm}\}$ . Inside this domain, a cylindrical capsule (perturbation) is located defined by $D\u2254\left\{{(x,y,z)}^{T}\right|{[{(x-3.0)}^{2}+{(y-3.0)}^{2}]}^{1\u22152}\u2a7d0.7\phantom{\rule{0.3em}{0ex}}\mathrm{cm};1.2\u2a7dz\u2a7d2.2\phantom{\rule{0.3em}{0ex}}\mathrm{cm}\}$ [see Fig. 5a]. The scattering coefficients and anisotropy factors of the background medium are identical to those of the cylindrical perturbation ( ${\mu}_{s}=40.0\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ , $g=0.73$ ). The absorption coefficients of the background and cylindrical perturbation are $0.1\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ and $0.2\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ , respectively. The source modulation frequency is set to $\omega =300\phantom{\rule{0.3em}{0ex}}\mathrm{MHz}$ . Using the mesh generating software GID, the computational domain is discretized into 13,798 elements with 2723 nodes. The angular domain is discretized into 24 uniformly distributed directions with a full-level symmetry. This discretization yields a total number of 13,798 unknown parameters, if we consider to reconstruct the absorption optical properties on the discretized elements only. In this case, the measurement data are generated by a synthetic forward computation rather than real experiments. Two layers of sources and detectors are placed around the boundary, as illustrated in Fig. 5a. The two layers are separated by $1.0\phantom{\rule{0.3em}{0ex}}\mathrm{cm}$ . On each layer 12 sources (S12) and 12 detectors (D12) are evenly distributed. In totally, we will generate $S24\times D24=576$ measurement data, which is far fewer than 13,798, the number of unknown discretized elements.

In a second set of studies, we consider another cylindrical domain, this time given by $D\u2254\left\{{(x,y,z)}^{T}\right|{[{(x-2.05)}^{2}+{(y-2.05)}^{2}]}^{1\u22152}\u2a7d2.05\phantom{\rule{0.3em}{0ex}}\mathrm{cm};0.0\u2a7dz\u2a7d5.0\phantom{\rule{0.3em}{0ex}}\mathrm{cm}\}$ . The cylindrical perturbation is given by $D\u2254\left\{{(x,y,z)}^{T}\right|{[{(x-3.0)}^{2}+{(y-3.0)}^{2}]}^{1\u22152}\u2a7d0.7\phantom{\rule{0.3em}{0ex}}\mathrm{cm};0.0\u2a7dz\u2a7d5.0\phantom{\rule{0.3em}{0ex}}\mathrm{cm}\}$ . Therefore, in this case both the height of the background cylinder and the perturbation are $5\phantom{\rule{0.3em}{0ex}}\mathrm{cm}$ [see Fig. 5b]. The perturbation has a uniform distribution along the $z$ direction. The discretization yields 23,221 unknown elements with 4512 nodes. We use 24 full-level symmetrical distribution directions as well. The examples chosen here only reconstruct images of absorption coefficients with fixed scattering coefficients. Scattering coefficients’ reconstruction can be implemented in a similar way.

## 4.2.

### Image Evaluation Metrics

In the following, we will compare reconstruction results of the parametric-DCT approach to those obtained through the conventional MOBIIR algorithm. To quantify the reconstruction quality and the differences between the two approaches, we use the following imaging evaluation metrics:

## Eq. 14

${\epsilon}_{\mathit{rms}}={\left[\frac{1}{N}\sum _{i=1}^{N}\left(\frac{{\mu}_{i}^{r}-{\mu}_{i}^{o}}{{\mu}_{i}^{o}}\right)\right]}^{1\u22152},$## 4.3.

### Image Reconstruction Results

## 4.3.1.

#### Reconstruction of a capsular perturbation

In this example, we varied the number of DCT coefficients to study their influence on the reconstruction results. To straddle the range of 10 coefficients as discussed in Sec. 3.3, we first fix the number of DCT coefficients in the $xy$ cross section ${K}_{x}\times {K}_{y}=10\times 10$ , and vary ${K}_{z}=5$ , 10, 15. Next, we fix ${K}_{z}=10$ and vary ${K}_{x}={K}_{y}=5$ , 10, 15. The initial absorption coefficients are set as ${\mu}_{a}=0.1\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ (background value).

Figure 6 shows the results obtained by varying the number of DCT coefficients in the $xy$ cross section. One can see that the reconstructed images with ${K}_{x}\times {K}_{y}\times {K}_{z}=5\times 5\times 10$ [Figs. 6d, 6e, 6f] reproduce the shape of the perturbation; however, the background shows somewhat larger fluctuations when compared to the results obtained with ${K}_{x}\times {K}_{y}\times {K}_{z}=10\times 10\times 10$ [Figs. 6g, 6h, 6i]. If we further increase ${K}_{x}\times {K}_{y}$ to $15\times 15$ [Figs. 6j, 6k, 6l], the image quality deteriorates.

Figure 7 shows the reconstruction results obtained by fixing ${K}_{x}\times {K}_{y}=10\times 10$ and varying the number of DCT coefficients along the $z$ direction. When ${K}_{z}=5$ , we observe [In Figs. 7d, 7e, 7f] that there are some fluctuation around the boundary. For ${K}_{z}=15$ , the reconstructed images show more fluctuation in the background. Table 1 provides the related maximum error and the normalized mean square error. Overall, one can see that the combination of ${K}_{x}\times {K}_{y}\times {K}_{z}=10\times 10\times 10$ gives the smallest errors. Therefore, we will use this number of coefficients ${K}_{x}\times {K}_{y}\times {K}_{z}=10\times 10\times 10$ in the following simulations, where we study the influence of the initial guess and noise levels on the reconstruction results.

## Table 1

Errors of the reconstructed images using the parametric-DCT method with different numbers of DCT coefficients.

Kx×Ky×Kz | 5×5×10 | 10×10×10 | 15×15×10 | 10×10×5 | 10×10×15 |
---|---|---|---|---|---|

${\epsilon}_{\mathit{max}}$ | 0.0499 | 0.0484 | 0.0632 | 0.0527 | 0.0498 |

${\epsilon}_{\mathit{rms}}$ | 0.0784 | 0.0744 | 0.0936 | 0.0749 | 0.0752 |

Next, we compare the results of the parametric-DCT reconstruction $({K}_{x}\times {K}_{y}\times {K}_{z}=10\times 10\times 10)$ with those of the conventional MOBIIR reconstruction (Figs. 8 and 9 ). Looking at Figs. 8d, 8e, 8f, we notice that using the conventional MOBIIR method the shape of the perturbation is distorted. Using the parametric-DCT method [see Figs. 8g, 8h, 8i], the shape of the perturbation is much smoother. The cross-line plots presented in Fig. 9 give the values of optical properties along the Cartesian coordinates crossing the center of the capsulized perturbation. This shows that the cross-line curves of the absorption values have a smooth bell shape when the parametric-DCT method is employed. Using the conventional MOBIIR method results in more fluctuating curves (Table 2 ).

## Table 2

Errors of the reconstructed images using the parametric-DCT method and the conventional MOBIIR method.

Capsule | Column | ||||
---|---|---|---|---|---|

Conventional | Parametric-DCT | Conventional | Parametric-DCT | ||

(10×10×10) | (10×10×1) | ||||

${\epsilon}_{\mathit{max}}$ | 0.0536 | 0.0484 | 0.0803 | 0.0563 | 0.0548 |

${\epsilon}_{\mathit{rms}}$ | 0.0770 | 0.0744 | 0.1581 | 0.1601 | 0.1435 |

## 4.3.2.

#### Reconstructions of a column perturbation

To compare experimental results with our numerical simulations, we consider a cylindrical medium with a cylindrical column perturbation [see Fig. 5b]. Given the results obtained in the previous section, we performed the reconstruction with with ${K}_{x}\times {K}_{y}\times {K}_{z}=10\times 10\times 10$ . In addition, assuming that we know that there are no variation in the optical properties along the $z$ direction, we also performed reconstructions with ${K}_{x}\times {K}_{y}\times {K}_{z}=10\times 10\times 1$ . Figure 10 shows the exact absorption maps [Figs. 10a, 10b, 10c], the reconstructed absorption maps using the parametric-DCT method [Figs. 10d, 10e, 10f, 10g, 10h, 10i], and the reconstructed absorption maps using the conventional method [Figs. 10j, 10k, 10l]. The cross-line plots of the optical properties of the reconstruction images are presented in Fig. 11 .

One can see that employing the prior information enhances the image quality [see Figs. 10d, 10e, 10f]. Fixing ${K}_{z}=1$ leads to a uniform distribution of absorption coefficients along the $z$ direction. If we ignore this prior information and reconstruct with parameters ${K}_{x}\times {K}_{y}\times {K}_{z}=10\times 10\times 10$ , we are still able to obtain smooth images. However, along the $z$ direction, the reconstructed perturbation has limited length and locates around the two source–detector planes. Information from above and below these two planes is limited. Last, we observe that the reconstructed images using the conventional MOBIIR method are much more noisy than the results obtained with the parametric image reconstruction method. The related values for the maximum errors and the root mean square errors are 0.0803 and 0.1581 for MOBIIR reconstruction and 0.0563 and 0.1601 for the parametric DCT method with $1010\times 10$ coefficients (see Table 2).

## 4.4.

### Effects of the Initial Guess on Reconstruction Results

At the heart of the inverse reconstruction for both the conventional MOBIIR and the parametric-DCT method is the quasi-Newton update, which is an optimization scheme for finding a local minimum. If the initial guess fed into the reconstruction algorithm is far away from the true minimum, it is difficult for the algorithm to converge to its real solution.^{17} To evaluate the sensitivity of the parametric-DCT reconstruction method to the choice of the initial guess, we performed numerical studies using the capsular perturbation. The number of DCT coefficients is fixed to
${K}_{x}\times {K}_{y}\times {K}_{z}=10\times 10\times 10$
.

Figure 12 presents the reconstructed cross-sectional images with initial guesses varying from ${\mu}_{a}=0.10\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ (background value) to ${\mu}_{a}=0.05\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ . As expected, we observe that the stronger the initial guess deviates from the background value, the worse the images’ qualities are. However, even an initial guess 50% lower than the actual values still leads to a reasonable result. This cannot be said when the conventional MOBIIR scheme is applied. Figures 12j, 12k, 12l show a much more noisy reconstruction than Figs. 12g, 12h, 12i. Table 3 shows the errors and mean values of the reconstructed images. We define the images mean value as ${\overline{\mu}}_{a}={\sum}_{i}^{N}{\mu}_{a}\left(i\right)\u2215N$ , with $N$ as the number of discretized elements. Comparing the parametric-DCT method with the conventional MOBIIR method at the initial guess ${\mu}_{a}=0.05\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$ , one can find that the parametric-DCT method yields smaller errors and a closer mean value to the exact image. We define the mean value offset as $\delta {\overline{\mu}}_{a}=|{\overline{\mu}}_{a}^{r}-{\overline{\mu}}_{a}^{o}|$ , with ${\overline{\mu}}_{a}^{r}$ and ${\overline{\mu}}_{a}^{o}$ as the mean values of the reconstructed and exact images, respectively. We can see the reconstructed mean values with the parametric-DCT approach are less offset, as shown in the last row of Table 3. Overall, these results show that the parametric-DCT approach is much less sensitive to the choice of the initial guess than the conventional MOBIIR scheme.

## Table 3

Errors and mean values of reconstructed images with different initial guesses in the capsular perturbation case.

Exact | Parametric-DCT | Conventional | |||
---|---|---|---|---|---|

${\mu}_{a}^{0}$ | N/A | 0.10 | 0.07 | 0.05 | 0.05 |

${\epsilon}_{\mathit{max}}$ | N/A | 0.0484 | 0.0555 | 0.0625 | 0.0954 |

${\epsilon}_{\mathit{rms}}$ | N/A | 0.0744 | 0.1596 | 0.2457 | 0.2930 |

${\overline{\mu}}_{a}$ | 0.1019 | 0.1023 | 0.0949 | 0.0903 | 0.0850 |

$\delta {\overline{\mu}}_{a}$ | 0 | 0.0004 | $-0.0070$ | $-0.0116$ | $-0.0169$ |

## 4.5.

### Effects of Measurement Noise on Reconstruction Results

To study the effect of measurement noise on the parametric-DCT reconstructions, we purposely add a certain amount of random noise to the synthetically generated data. If ${J}_{i,d}$ is a synthetic measurement corresponding to the source $i$ and the detector $d$ , then the noisy data ${\overline{J}}_{i,d}={J}_{i,d}(1+{\delta}_{i,d}\mathcal{N}(0,1){I}_{[-1,1]})$ . Here, $\mathcal{N}(0,1)$ is a standard normal distribution, and $I[-1,1]$ is an interval indicator function that is equal to 1 in the interval of $[-1,1]$ and zero otherwise. The parameter ${\delta}_{i,d}$ is the standard deviation of added noise, which we vary from 1% to 5%, which are typical values for OT measurement systems. We test the effect of this type of noise on the reconstruction results of a cylindrical domain with a capsulized perturbation (Fig. 5). Again, we chose the number of DCT coefficients as ${K}_{x}\times {K}_{y}\times {K}_{z}=10\times 10\times 10$ .

The reconstruction results are shown in Fig. 13 . As expected, we see that as the standard deviation of the added noised increases, the quality of the reconstructed images decreases. When we add 5% noise, some artificial perturbations appear in the background area, while the reconstructed absorption values decrease in perturbation area. However, we also observe that the parametric-DCT reconstruction approach outperforms the conventional MOBIIR method at the same noise level [see Figs. 13g, 13h, 13i, 13j, 13k, 13l]. Using the conventional approach, we see artificial inhomogeneities in background area and the shape of reconstructed perturbation appears distorted. The maximum errors and normalized standard deviation errors with varying noise levels are presented in Table 4 . Overall, we find that the parametric-DCT method is less sensitive to noise.

## Table 4

Errors of reconstructed images with different noise level in the capsule perturbation case.

Parametric-DCT | Conventional | |||
---|---|---|---|---|

Noise level | No noise | 2% | 5% | 5% |

${\epsilon}_{\mathit{max}}$ | 0.0484 | 0.0611 | 0.1103 | 0.3111 |

${\epsilon}_{\mathit{rms}}$ | 0.0744 | 0.1299 | 0.1596 | 0.2655 |

## 5.

## Image Reconstruction with Experimental Data

Last, we also tested our new parametric imaging reconstruction code with experimental data. The experimental data were acquired with a continuous wave dynamic near-infrared optical tomographic (DYNOT) system.^{23} In this system, light from a laser diode (wavelength
$\lambda =760\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
is delivered to 24 source fiber bundles via an optical demultiplexing switch that allowed light to be sequentially delivered to different positions. The demultiplexor consists of a mirror mounted on a rotating stepper motor, which reflected light from the incoming laser source into different source fibers depending on the angular orientation of the mirror. Twenty-four detector fiber bundles were used to direct light to the detectors. Light was collected by silicon photodiodes, which were connected to a series of programmable amplifiers and a lock-in filter. Reported signal-to-noise ratio (SNR) is in the range of 0.5% to 5%, depending on the optical properties and tissue geometries. More details concerning this system can be found in Ref. 23.

The tissue-like phantom used in this study has a similar geometry compared to the preceding numerical setup with a column perturbation [see Sec. 4.3.2, Fig. 5c]. The background medium of this phantom consists of a 1% Intralipid solution. The cylindrical perturbation is made from a mix of 1% Intralipid and ink placed in a transparent straw. The scattering coefficients of the background and the perturbation are the same
${\mu}_{s}^{b}={\mu}_{s}^{p}\sim 27.0\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$
(Ref. 24). The absorption coefficients of the background and the perturbation are
${\mu}_{a}^{b}\sim 0.023\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$
and
${\mu}_{a}^{p}>0.1\phantom{\rule{0.3em}{0ex}}{\mathrm{cm}}^{-1}$
, respectively.^{25} The anisotropic factor
$g\sim 0.73$
is close to water. The diameter of this perturbation
$d\sim 0.6\phantom{\rule{0.3em}{0ex}}\mathrm{cm}$
is a little smaller than the size of numerical perturbation used in sections above.

The reconstructed absorption maps and corresponding cross-line plots are shown in Figs. 14 and 15 , respectively. For the parametric-DCT reconstruction, We used ${K}_{x}\times {K}_{y}\times {K}_{z}=10\times 10\times 1$ DCT coefficients. We can see that the parametric-DCT approach outperforms the conventional MOBIIR method. Both the shape and absolute optical properties of the inhomogeneity are better recovered with the parametric-DCT method. This is in agreement with the findings of our numerical studies.

## 6.

## Discussion and Conclusions

We introduced a parametric image reconstruction approach to optical tomographic imaging, which makes use of the discrete-cosine transform (DCT) to represent the reconstructed images. This approach was implemented by adopting an existing model-based iterative image reconstruction (MOBIIR) scheme, which uses the frequency-domain equation of radiative transfer (ERT) as a forward model. We illustrated and analyzed the performance of the new code using synthetic as well as experimental data. We could demonstrate that the parametric-DCT reconstruction method allows us to substantially reduce the number of unknown image parameters and overall results in better quality of the reconstructed images, when compared to a conventional, nonparametric reconstruction techniques. In particular, the results show that DCT-based codes are less sensitive to noise in the data and to the choice of the initial guess needed in iterative image reconstruction schemes.

The number of unknowns in the imaging problem is reduced by expanding the reconstructed images into spatially varying 2-D or 3-D cosine functions of increasing order. If a typical image is discretized into pixels or voxels, the number of unknowns can be in the thousand and even millions, depending on the level of discretization. However, it has been shown by several groups that the detection limit of a heterogeneity is about $1\u221510$ of the outer dimensions of the medium under consideration. This effectively means that variations of spatial properties in optical tomographic images are limited to about $1\u221510$ of the outer dimensions. Our numerical results show that these variations can be captured with approximately 10 coefficients in each dimension, since the largest spatial frequency is proportional to one tenth of the base frequency. Therefore, using the parametric-DCT approach, only 10 unknowns (the expansion coefficients) per dimension need to be determined instead of the much larger number of pixels in an image.

If more and more expansion coefficients are used in the parametric-DCT algorithm, the images become increasingly similar to the images obtained with traditional nonparametric MOBIIR codes. In this case, the number of unknowns in parametric-DCT reconstruction approaches the number of unknowns in the MOBIIR method. Higher frequency components that seem to produce artifacts become visible. The reason for that is most likely that an overfitting takes places, and the many unknowns will “fit” the noise in the data. On the other hand, as we use fewer and fewer coefficients in the parametric code, the images become smoother, and high-frequency components disappear. Since fewer unknowns are used, images reconstructed with the parametric-DCT method appear in general smoother than images generated with MOBIIR codes that use a larger number of unknowns.

Another advantage of the parametric DCT approach is that these codes converge faster to the final result. For the examples considered in this paper, we observed up to 50% reduction in convergence times. The exact speed-up factor depends on the optical properties, problem geometry, and the ratio of the mesh size to the number of chosen DCT coefficients.

In this paper, we present results concerning the reconstruction of the spatial distribution of the absorption coefficients. However, the same approach can be applied to reconstruct the spatial distribution of scattering coefficients. As indicated in Eq. 9, the absorption coefficients and the scattering coefficients of an image can be expanded in the spatial-frequency domain individually. Thus, the number of DCT coefficients for the absorption coefficients and for the scattering coefficients become independent. The number of necessary DCT coefficients to obtain a good scattering coefficients map will depend on the expected spatial variations in a reconstructed image. It has been shown that the spatial variations in the scattering images are of the same order as the variations in the absorption images (see, e.g., Dehghani
^{21} and Pogue
^{22}). Therefore, one can expect that a similar number of DCT coefficients are needed for scattering images, and indeed we found the same general trends for scattering images as described in the paper for absorption images.

Last, although the examples presented in this paper were based on optical tomography using the frequency-domain ERT, the parametric-DCT method by nature is not limited to any particular light-propagation model. The parametric-DCT reconstruction techniques can also be applied to diffusion-theory-based optical tomography and other inverse problems, such as impedance tomography and microwave tomography.

## Acknowledgments

This work was supported in part by the National Institute of Biomedical Imaging and Bioengineering (NIBIB—Grant No. R01 EB001900) and the National Institute of Arthritis and Musculoskeletal and Skin Diseases (NIAMS—Grant No. 2R01-AR046255), both of which are part of the National Institutes of Health (NIH).

## References

**,” Opt. Lett., 28 (23), 2339 –2341 (2003). https://doi.org/10.1364/OL.28.002339 0146-9592 Google Scholar**

*Uniqueness and wavelength optimization in continuous-wave multispectral diffuse optical tomography***,” J. Biomed. Opt., 11 (1), 014020 (2006). https://doi.org/10.1117/1.2161199 1083-3668 Google Scholar**

*Design and implementation of a multifrequency near-infrared diffuse optical tomography system***,” Neuroimage, 23 (S1), S275 –288 (2004). https://doi.org/10.1016/j.neuroimage.2004.07.011 1053-8119 Google Scholar**

*Diffuse optical imaging of brain activation: approaches to optimizing image sensitivity, resolution, and accuracy***,” Appl. Opt., 43 (8), 1654 –1662 (2004). https://doi.org/10.1364/AO.43.001654 0003-6935 Google Scholar**

*A dual-mesh optical tomography reconstruction method with depth correction usign**a priori*ultrasound information**,” J. Electron. Imaging, 12 (4), 583 –593 (2003). https://doi.org/10.1117/1.1586919 1017-9909 Google Scholar**

*Image reconstruction in optical tomography using local basis functions***,” Commun. ACM, 34 (4), 30 –44 (1991). https://doi.org/10.1145/103085.103089 0001-0782 Google Scholar**

*The jpeg still picture compression standard***,” SIAM J. Sci. Comput. (USA), 28 (4), 1463 –1489 (2006). https://doi.org/10.1137/040619193 1064-8275 Google Scholar**

*Frequency-domain optical tomography with the equation of radiative transfer***,” Astrophys. J., 90 (16), 70 –83 (1941). https://doi.org/10.1086/144246 0004-637X Google Scholar**

*Diffuse radiation in the galaxy***,” Handbook of Numerical Analysis VII, 2nd ed.North-Holland, Amsterdam (2000). Google Scholar**

*Finite volume methods***,” Opt. Lett., 29 (6), 578 –580 (2004). https://doi.org/10.1364/OL.29.000578 0146-9592 Google Scholar**

*Algorithm for solving the equation of radiative transfer in the frequency domain***,” IEEE Trans. Image Process., 1 (1), 18 –30 (1992). https://doi.org/10.1109/83.128028 1057-7149 Google Scholar**

*Image coding based on a fractal theory of iterated contractive image transformations***,” J. Mach. Learn. Res., 7 2673 –2698 (2006). 1532-4435 Google Scholar**

*Spam filtering using statistical data compression models***,” Opt. Lett., 18 (19), 1591 –1593 (1993). https://doi.org/10.1364/OL.18.001591 0146-9592 Google Scholar**

*Resolution limits for imaging through turbid media with diffuse light***,” Med. Phys., 24 (3), 361 –368 (1997). https://doi.org/10.1118/1.597904 0094-2405 Google Scholar**

*Evaluation of spatial resolution as a function of thickness for time-resolved optical imaging of highly scattering media***,” Appl. Opt., 36 (1), 75 –92 (1997). https://doi.org/10.1364/AO.36.000075 0003-6935 Google Scholar**

*Detection and characterization of optical inhomogeneities with diffuse photon density waves: a signal-to-noise analysis***,” Appl. Opt., 42 (16), 3117 –3128 (2003). https://doi.org/10.1364/AO.42.003117 0003-6935 Google Scholar**

*Three-dimensional optical tomography: resolution in small-object imaging***,” J. Biomed. Opt., 11 (3), 033001 (2005). https://doi.org/10.1117/1.2209908 1083-3668 Google Scholar**

*Image analysis methods for diffuse optical tomography***,” Rev. Sci. Instrum., 73 (2), 429 –439 (2002). https://doi.org/10.1063/1.1427768 0034-6748 Google Scholar**

*Instrumentation for fast functional optical tomography***,” Phys. Med. Biol., 34 1927 –1930 (1989). https://doi.org/10.1088/0031-9155/34/12/015 0031-9155 Google Scholar**

*The optical properties of aqueous suspensions of intralipid, a fat emulsion***,” Phys. Med. Biol., 37 985 –993 (1992). https://doi.org/10.1088/0031-9155/37/4/012 0031-9155 Google Scholar**

*The use of india ink as an optical absorber in tissue-simulating phantoms*