1.

Introduction

Optical tomography (OT) is a fast growing field in which near-infrared light is used to image the distribution of optical properties inside the human body. Optical properties of interest are, for example, the absorption coefficient μ_a, the reduced scattering coefficient μ_s ^′, or the diffusion coefficient D=c/(3μ_a+3μ_s ^′), where c is the speed of light in the medium. The instrumentation for making light transmission measurements that are necessary for OT is nowadays widely available.^{1 2 3 4 5 6} Furthermore, several algorithms that transform these measurements into useful cross-sectional images have matured to such a degree that first clinical trials are underway, especially in breast imaging.^{7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30} However, a major difficulty in OT remains that the image reconstruction problem is ill-posed or underdetermined. In other words, there are many distributions of optical properties inside the medium under investigation that lead to the same set of detector readings on the surface of the medium.

Depending on the underlying structure and concepts of any particular reconstruction scheme, the problem of ill-posedness can be approached in different ways. Most of the currently employed reconstruction schemes fall in one of two classes, which we refer to as the linear-perturbation approach and the nonlinear-gradient method (see also Ref. 30). In both cases the goal is to minimize the difference between intensities measured on the boundary of the medium and some prediction for those measured values. Therefore the image reconstruction problem may be interpreted as an optimization problem in which an objective function is minimized. If we use the χ² error norm to calculate the difference between measurements and predictions we can define the objective function as

In the linear perturbation approach, it is assumed that we have an estimate ζ ₀ that is close to the true distribution ζ. In this case we can perform a Taylor expansion of P_s,d around ζ ₀. Neglecting nonlinear terms, we obtain

More recently several groups (Saquib et al.;,²⁴ Hielscher et al.;,^{25 27 28} Arridge and Schweiger,²⁶ Roy and Sevick,²⁹ and Ye et al.;³³) have developed so called gradient based iterative reconstruction (GIIR) schemes that do not solve Eqs. (4) or (5) to obtain an update of ζ. Instead, starting from an initial guess ζ₀ subsequent distributions ζ _k+1 are obtained by calculating

The question arises: what does ill-posedness mean using the GIIR approach? In the perturbation approach ill-posedness is identical to an ill-conditioned matrix JJ ^T. In the GIIR scheme, which interprets image reconstruction as a multidimensional minimization problem, ill-posedness means that a global minimum is not well defined. A global minimum is not well defined when either many minima with similar values exist or one minimum ζ _min , which is surrounded by many other ζ, result in almost identical objective functions. Examples of such objective functions in optical tomography have been described, for example, by Arridge,³⁰ Schweiger and Arridge³⁴ and Hielscher et al.;³⁵ In general, the reconstruction result will strongly depend on the initial guess ζ ₀. These phenomena are well known in the field of general optimization theory. A typical way to overcome this problem is to add penalty terms that provide additional constraints on the solution space. The penalty term may push the solution into the right area of the solution space and the minimization of the resulting objective function may provide a better-defined minimum. While the use of penalty terms has been studied for a variety of problems,^{36 37 38 39 40 41 42 43 44} their use in GIIR schemes for optical tomography has not been explored.

The goal of this paper is to introduce and study the effects of various penalty functions on the reconstruction process in OT. In particular, we seek to derive penalty function from a priori knowledge about the systems under investigation. A priori knowledge is, in general, information that does not depend on the difference between predicted and measured data. For example, we know that optical properties are larger than zero. Furthermore in the near-infrared wavelength region μ_a<10 cm ⁻¹ and μ_s ^′<100 cm ⁻¹. In addition, it is often known how many different tissue types, and therefore how many different optical properties, are present in the interrogated medium. For example, in the brain we find cerebrospinal fluid, gray, and white matter. If magnetic resonance imaging data are available, one may even know the location of these tissues. The questions arise: how can this additional a priori knowledge be appropriately cast into a penalty function and can using this information improve image quality?

To illustrate the effects of penalty functions derived from a priori knowledge about the system, we will focus in this work on information regarding the composition of a given tissue volume. We consider the following cases: First, it is assumed that one knows that n tissue types with n different optical properties are present in the medium. However one neither knows the location nor their respective volume fraction of any given component in the tissue sample. Second, we will assume to have prior knowledge of all tissue types present and the volume percentage they occupy. However, just as in the first case, one does not know where the different tissue types are located inside the medium. This latter point corresponds to knowing the histogram of the medium to be reconstructed. We will derive the appropriate penalty terms for both cases. Numerical examples that demonstrate the effect of these penalty terms on the quality of the reconstructed images are shown and discussed. Furthermore, the relationship between penalty functions in gradient-based schemes and algebraic regularization mechanisms used in linear-perturbation schemes is discussed and put into context with our results.

2.

Mathematical Background

2.2.

Penalty Functions

The incorporation of additional knowledge about the image to be reconstructed can be achieved by including a penalty term into the definition of the objective function. Instead of only trying to minimize the difference between the measurements and predictions, as is done if the objective function equals the χ² error norm defined in Eq. (1), one can define the objective function as

Eq. (7)

$$\Phi (\zeta )={\chi }^2(\mathbf{\text{ζ}})+\omega \cdot \Pi (\mathbf{\text{ζ}})\text{,}$$

where Π(ζ) denotes the penalty term and ω is the coupling parameter, also often referred to as the hyperparameter. The additional term ωΠ(ζ) is designed to provide a better-defined minimum in the solution space and push the gradient scheme toward more probable solutions. The latter is achieved by penalizing certain distributions ζ that are unlikely given our prior information of the system under investigation. Ideally, we expect Π(ζ) to constrain the N-dimensional minimum search to a small subspace of desirable solutions.

The penalty term should be continuously differentiable to be useful in GIIR schemes. Therefore the derivative of the penalty term with respect to the optical properties should not have any discontinuities. In general, smooth functions are much better suited for any gradient based minimum search. We will address these restrictions further when we derive the respective penalty functions.

2.2.1.

Tissue-Type Penalty Functions

In this work we seek to provide penalty functions that are rooted in a particular kind of prior knowledge, namely the different types of tissues that are most likely to be found in the sample under investigation. The optical properties of these tissues are known within certain error margins. Areas of the reconstructed image that show any other than the expected optical parameters should therefore be penalized the more they differ from the expected values. Under these considerations we define the following penalty function:

Eq. (8)

$${\Pi }^{\text{tt}}(\mathbf{\text{ζ}})=\sum _{x\in S}(1-\sum _{k=1}^K\text{exp}(-\frac{{(a_k-{\zeta }_x)}^2}{2{\sigma }_k^2}))\text{.}$$

Here K is the number of different tissue types k in the system S. The parameter a_k is the most likely optical property of tissue k, and ζ_x is the reconstructed optical property at pixel x. The parameter 1/σ_k can be interpreted as the confidence that a_k is the exact value. For 1/σ_k→∞ (⇔σ_k→0) the Gaussian becomes a delta function that is 1 only if ζ_x=a_k and zero otherwise. Therefore, only if ζ_x=a_k does the penalty term disappear. In this case Eq. (8) is not continuously differentiable. For any nonzero value of σ_k there will be a range of ζ_x for which the penalty term changes smoothly from 1 to 0 and Eq. (8) becomes continuously differentiable. An example is shown in Figure 1, where Π ^tt is plotted for the case that three different tissue types with optical properties D₁, D₂, and D₃ are present and σ_k=0.1 cm ² ns ⁻¹ for all three D values. The derivative of this penalty function with respect to a pixel value ζ_x is easily calculated as

Figure 1

Example of a tissue-type penalty function Π ^tt [Eq. (8)] that assumes three different types of tissues (k=3), with D₁=0.52 cm ² ns ⁻¹, D₂=0.92 cm ² ns ⁻¹, and D₃=1.44 cm ² ns ⁻¹.

Eq. (9)

$$\frac{\partial {\Pi }^{\text{tt}}\left(\mathbf{\text{ζ}}\right)}{\partial {\zeta }_x}=-\sum _{k=1}^K\frac{a_k-{\zeta }_x}{{\sigma }_k^2}\text{exp}(-\frac{{(a_k-{\zeta }_x)}^2}{2{\sigma }_k^2})\text{.}$$

2.2.2.

Histogram Penalty Function

If not only the various types of tissue in the sample are known but additionally their respective volume fraction, we are able to provide a most likely histogram H⁰ to the reconstruction scheme. In general, a histogram maps the N pixel values of an image represented by ζ=(ζ₁,ζ₂,…,ζ_x,…,ζ_N) onto L discrete intervals or bins of width Δζ_l=ζ_l−ζ_l−1(l=1,2,…,L). Therefore, the histogram H(ζ_l,ζ) associated with an image S or given set ζ of optical properties is defined by the following sum over all pixels x of the image

Eq. (10)

$$H({\zeta }_{\mathrm{l}},\mathbf{\text{ζ}})\equiv \sum _{x\in S}\delta ({\zeta }_{\mathrm{l}},{\zeta }_x),\mathrm{l}=1,\dots ,L\text{,}$$

where

Eq. (11)

$$\delta ({\zeta }_{\mathrm{l}},{\zeta }_x)\equiv \{\begin{array}{cc}1\hfill & \text{for}{\zeta }_{l-1}<{\zeta }_x\le {\zeta }_{\mathrm{l}}\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}\text{.}$$

The sum in Eq. (10) sorts all pixel values ζ_x into L bins of the histogram, since the δ-function only contributes to the sum H(ζ_l,ζ) if a pixel lies within the corresponding interval.

We can now define a penalty term Π ^Hist , which evaluates the histogram H(ζ_l,ζ) of a reconstructed image relative to the expected histogram H⁰(ζ_l,ζ). Here we choose the χ² error norm between both functions and define the penalty term

Eq. (12)

$${\Pi }^{\text{hist}}\equiv {\chi }_{H^0,H}^2=\sum _{l=1}^L\frac{{(H^0({\zeta }_{\mathrm{l}})-H({\zeta }_{\mathrm{l}}))}^2}{H^0({\zeta }_{\mathrm{l}})}\text{,}$$

where we use the short notation H(ζ_l) for H(ζ_l,ζ).

However, if the definition of the delta function in Eq. (11) is considered, we observe that Eq. (12) is not continuously differentiable, which makes it unsuitable for gradient based schemes. To overcome this problem we consider a representation of the delta function that uses a Gaussian

Eq. (13)

$$\delta ({\zeta }_{\mathrm{l}},{\zeta }_x)=\underset{\sigma \to 0}{\text{lim}}\sqrt{1/(2\pi {\sigma }^2)}\text{exp}(-\frac{{({\zeta }_{\mathrm{l}}-{\zeta }_x)}^2}{2{\sigma }^2})\text{.}$$

With this we can calculate the derivative and obtain

Eq. (14)

$$\frac{\partial {\Pi }^{\text{hist}}}{\partial {\zeta }_x}=-\sum _{l=1}^{L}2\sqrt{1/\left(2\pi {\sigma }^2\right)}\frac{\left(H^0\right({\zeta }_{\mathrm{l}})-H({\zeta }_{\mathrm{l}}\left)\right)}{H^0\left({\zeta }_{\mathrm{l}}\right)}\frac{({\zeta }_{\mathrm{l}}-{\zeta }_x)}{{\sigma }^2}\text{exp}(-\frac{{({\zeta }_{\mathrm{l}}-{\zeta }_x)}^2}{2{\sigma }^2})\text{,}$$

where we omitted the lim. By choosing an appropriate nonzero value for σ and calculating the image histogram using Eq. (10) and the penalty term Eq. (12), we obtain a penalty term that is continuously differentiable and can be applied within the GIIR framework. Typically, we chose σ, so that the singular peaks overlapped (see Figure 2) to avoid forbidden intervals in the target histogram H⁰.

Figure 2

Ideal histogram H_δ for tissue that contains 15 tissue type 1 (D₁=0.520 cm ² ns ⁻¹), 75 tissue type 2 (D₂=0.92 cm ² ns ⁻¹), and 10 tissue type 3 (D₃=1.44 cm ² ns ⁻¹) and smoothed histogram H. The smoothed histogram results from convolution of H_δ with a Gaussian function of width σ=0.15 [Eq. (10) and Eq. (13)].

Using Eq. (13) with a nonzero σ to generate H(ζ_l) results in a smoothing operation on the histogram and is equivalent to convolving the exact histogram with a Gaussian of width σ. It may seem that the histogram thus obtained is blurred to an extent that spoils the goal of reconstructing only certain types of tissue. However, since we perform the same convolution on both, the histogram H of the reconstructed image and the target histogram H⁰ before they are compared [Eq. (12)] becomes minimal, if and only if the unconvolved functions are identical. Hence, let (H_δ ⁰,H_δ) be accurate histograms [Eq. (10)] and (H⁰,H) be their convolution with Eq. (13) (with nonzero σ), then

Eq. (15)

$$H_{\delta }^0=H_{\delta }\iff H^0=H\text{.}$$

Consequently, by minimizing the χ² norm between two blurred histograms, we are at the same time minimizing the difference between the underlying sharp histograms. However, using the blurred histogram provides us with a continuously differentiable penalty term that is better suited for use in gradient based minimization schemes.

2.2.3.

The Hyperparameter

Each of the penalty functions has to be coupled to the error norm χ² [Eq. (1)] with a hyperparameter. This parameter fixes the relative strength of the penalty term in the minimization scheme and describes the confidence one has that the additional information is correct. A good starting point is to choose ω in a way which ensures that the gradients of both the χ² term and the penalty function are of similar magnitude. If we interpret the χ² term and the ωΠ term as potentials, the derivatives ∂χ²/∂ζ and ω∂Π/∂ζ can be interpreted as two forces that pull the pixel values in certain directions. By defining a new hyperparameter

Eq. (16)

$$\gamma =\left|\frac{\omega \partial \Pi }{\partial \mathbf{\text{ζ}}}\right|/\left|\frac{\partial {\chi }^2}{\partial \mathbf{\text{ζ}}}\right|\text{,}$$

we can adjust the relative strength of these forces, originating from the χ² potential and the Π potential, respectively. If γ=1 then ω equals the ratio of the gradients and both forces are equally strong. If γ>1 the force due to penalty function ωΠ is stronger than due to the χ² term, and if γ<1 the χ² term has a stronger influence. We will later see how the choice of γ influences the reconstruction results.

3.

Results

3.1.

Problem Setup

To test how assumptions about the medium can improve the image reconstruction we consider the following example. Given is an 8×8 cm medium that contains two objects (see Figure 3). The speed of light in this medium is considered to be constant c=22 cm ns ⁻¹. The optical properties of the background medium are given by μ_a=0.1 cm ⁻¹, μ_S ^′=8 cm ⁻¹ (⇒D=c/(3μ_a+3μ_s ^′)=0.905 cm ² ns ⁻¹). The optical properties of the inclusions are given by μ_a=0.1 cm ⁻¹, μ_S ^′=14.0 cm ⁻¹ (⇒D=0.520 cm ² ns ⁻¹), and μ_a=0.1 cm ⁻¹, μ_S ^′=5.0 cm ⁻¹ (⇒D=1.438 cm ² ns ⁻¹), respectively. The medium is discretized into a 40×40 mesh with Δx=0.2 cm. Four sources, one centered at each side, surround the medium. For each source 20 detector readings are available, four on each side and one on each corner. These detector readings were simulated by using the time-dependent finite-difference forward model and adding Gaussian noise with a signal to noise ratio of 30 db to the result. The simulated detector readings each consists of 100 time-dependent fluence rates, with Δt=0.01 ns. Therefore we have 4×20×100=8000 data points. From these points the χ² term is calculated as

Eq. (17)

$${\chi }^2=\sum _s\sum _d\sum _t{\left(\frac{M_{s,d}^t-P_{s,d}^t(D,{\mu }_a)}{P_{s,d}^t(D,{\mu }_a)}\right)}^2\text{,}$$

where D and μ_a are vectors of length 40×40=1600, which contain the diffusion and absorption coefficients throughout the medium.

Figure 3

(a) Composition of two-dimensional example problem. An 8×8 cm domain is divided into a 40×40 grid. Source positions are indicated with white circles and detectors with black circles. (b) Reconstruction without penalty functions after 20 iterations. The initial guess for the first iteration is a homogeneous medium with D=1 cm ⁻² ns ⁻¹.

To quantify the quality of the image reconstruction, we furthermore define the relative image error as

Eq. (18)

$$I_E\equiv \frac{1}{N}\sum _n^N\sqrt{\frac{{(D_n^{\text{target}}-D_n^{\text{wp}})}^2}{{(D_n^{\text{target}}-D_n^{\text{wop}})}^{2\text{’}}}}$$

where N is the number of pixels in the image, D_n ^target is the correct diffusion coefficient at pixel position n, D_n ^wop is the diffusion coefficient at pixel position n reconstructed without using a penalty term, and D_n ^wp is the diffusion coefficient at pixel position n reconstructed using a penalty term. Defining the image error this way allows us to quantify the improvement in image quality when using a penalty function as compared to not using a penalty function. A value of I_E=1 indicates that the penalty term yields neither an improvement nor degradation in image quality.

3.3.

Reconstruction with Tissue-Type Penalty Term

Next we consider how these results can be improved when the information is used that only three different tissues types with given optical properties are present in the medium [Eq. (8)]. However, no knowledge is assumed about where these three different regions are located, or how much of the image volume is occupied by a particular tissue type. Figures 4(a)–4(c) show reconstructions obtained with the penalty term Π ^tt for three different parameters γ. In this example we chose D₁=0.520 cm ² ns ⁻¹, D₂=0.92 cm ² ns ⁻¹, and D₃=1.44 cm ² ns ⁻¹. Furthermore, we set the width of the Gaussians to σ_k=0.1 cm ² ns ⁻¹ [Eq. (8)]. The initial guess is a homogeneous medium with D=1 cm ² ns ⁻¹ and μ_a=0.1 cm ⁻¹. Reconstructions were terminated after 20 iterations. It can be seen that the penalty term has an effect on the image for γ=0.02 [Figure 4(a)]. Increasing γ to 0.08 does produce better results [Figure 4(b)], however, a further increase to γ=0.1 yields reconstructions with almost constant D values across the entire image [Figure 4(c)]. All pixel values are close to the background value D=0.92 cm ² ns ⁻¹. Figure 5 shows the dependence of the image error I_E on the parameter γ. It can be seen that only for a small region of 0.01<γ<0.1 does the image quality improve (I_E<1). For all other values γ>0.1 the image error is larger than without the penalty term (I_E>1).

Figure 4

Image reconstruction with tissue-type penalty terms [Eq. (8)] after 20 iterations. The initial guess is a homogeneous medium with D=1 cm ⁻¹: (a) γ=0.02, (b) γ=0.08, (c) γ=0.1, and (d) γ=0.005 n², where n equals the number of iterations.

Figure 5

Image error I_E as function of the hyperparameter γ for the tissue type prior Π ^tt (open circles) and histogram prior Π ^hist (filled circles). Both curves are normalized to 1 for the result of the reconstruction without penalty term (γ=0).

Rather than keeping γ fixed for all iterations, one can also change this parameter dynamically. Figure 4(d) shows the result of increasing γ with the number of iterations n, so that γ=0.005 n². After 25 iterations both the location and D values of the inclusions are correctly recovered. Deviations between the reconstruction and the original image are on the order of single pixels, and the image error is relatively small (I_E=0.18).

3.4.

Reconstruction with Histogram Prior

Next, we consider the case where we have knowledge about the different tissue types in the medium and in addition know their respective volume fractions. Therefore, we assume to have knowledge about the histogram of the cross-sectional image of the tissue sample. Figures 6(a)–6(d) show reconstructed images using the histogram prior with increasing hyperparameter γ. As in previous cases, the initial guess is a homogeneous medium with D=1 cm ² ns ⁻¹ and μ_a=0.1 cm ⁻¹. Choosing γ<0.1 does not result in a significant improvement of the reconstruction and yields images similar to the unbiased case γ=0 [Figure 3(b)]. For γ⩾0.2 the two objects become more localized and show plateaus of constant D values. For γ=0.5 we observe increasingly sharp boundaries separating the heterogeneities from the background medium. Although their shape does not exactly match the original objects they reflect the correct size or rather volume fraction of these objects, as enforced by the histogram prior. Furthermore, the reconstructed D values remain constant across most of the areas covered by the cubes, rather than showing a pronounced peak, as observed without the penalty term. This allows for a much better extraction of the absolute optical parameters from the reconstruction. As the strength of the prior is increased further, we observe only a weak dependence of the image quality on γ [Figures 5, 6(c) and 6(d)]. For large values of γ, the reconstruction tends to exaggerate sharp edges and produces artifacts, while sacrificing the correct shape of the objects. As expected, the diffusion coefficient eventually assumes one of the three values imposed by the histogram. In Figure 6(d), the prior information is weighted 10⁶ stronger than the χ² term, so that fitting the image’s histogram to the prediction becomes the main objective. Nevertheless we still obtain reasonably good agreement with the original image.

Figure 6

Image reconstruction with histogram penalty terms [Eq. (12)] after 20 iterations. The initial guess is a homogeneous medium with D=1 cm ⁻¹: (a) γ=0.2, (b) γ=0.5, (c) γ=1.0, and (d) γ=10⁶.

4.

Discussion

In general we found that the effectiveness of a penalty term depends as much on the right choice of the penalty term as on the right choice of the hyperparameter, which adjusts the influence of the penalty term during the reconstruction process. The results for the tissue-type penalty terms Π ^tt show that an improvement in the reconstruction result can only be obtained for a relatively small range of hyperparameters (0.01<γ<0.1). For values of γ<0.01, the additional penalty term has no effect, while for γ>0.1 the quality of the reconstruction is actually worse than without penalty function. Invoking the earlier discussed interpretation of the χ² term and penalty term as potentials and there derivatives as forces, we see that if γ is chosen too strong the tissue-type penalty term traps all pixel values in the minimum closest to the initial guess (D=1 cm ² ns ⁻¹, Figure 1). Forces due to the χ² term are too weak to lift a pixel out of the valley centered around D=0.92 cm ² ns ⁻¹.

Choosing the hyperparameter dynamically can improve the convergence of the GIIR scheme with Π ^tt -type penalty terms, as shown in Figure 4(d). Giving only little weight to the penalty term in the early stages of the reconstruction process allows the algorithm to find areas with increased or decreased values of D based on information from the χ² term. Increasing the weight of the penalty term later in the reconstruction process leads to a sorting of each pixel in one of the three “potential” wells of the penalty term. However, even the use of a dynamically adjusted penalty term depends strongly on the appropriate choice of the initial value of γ and the rate at which it is increased. One is forced to find these parameters empirically for each particular reconstruction problem. The same drawback was previously encountered for total variation minimization penalty terms (see Ref. 17).

The histogram penalty terms Π ^hist provides good reconstruction results for all γ>0.01. Surprisingly even the very large γ=10⁶ provides a good reconstruction. The information contained in the histogram prevents the scheme from converging toward singular solutions, such as trapping all pixels at the same value of D. If more and more pixels assume only one of the expected values, the force acting on them diminishes and eventually becomes repulsive. Pixel values are pushed by penalty-term forces to alternative values to fit the overall histogram distribution. When the histogram of the reconstructed image approaches the “true” histogram the force acting on individual pixels diminishes and eventually becomes zero. At the global minimum of Π ^hist we may change single pixel values without increasing the penalty, by raising and lowering several pixels simultaneously. This is different in the scheme that uses the tissue-type penalty term Π ^tt . To interchange two pixel values that are already located at the bottom of two different potential valleys (Figure 1), one always needs a force that pushes each pixel value over the potential hill in-between.

The weak dependence of the reconstruction on γ and the fact that sensible results are still obtained as the penalty term becomes dominant, are important features of the histogram regularization. This makes it a very stable method to incorporate additional information into the reconstruction process. However, the full histogram information may not always be available, and only the tissue-type penalty term, which contain less information, may have to be used.

An additional effect of applying the histogram penalty function is a faster conversion of the reconstruction algorithm. Typically only 20 iterations are necessary to produce qualitatively correct images, compared to approximately 50 iterations without this prior information. As we increase γ to values ≫1, as little as 10 iterations suffice to produce the sharply contrasted results shown in Figure 6(d).

Given that the target histogram is known, one could be tempted to start with an image that has the right histogram. Therefore, rather than using a homogenous medium as the initial guess, one could chose a heterogeneous image as the initial guess, which has the correct volume fraction of certain tissue types. In Figures 7(a) and 7(b), a scrambled version of the exact image was used as input, so that during the first iteration Π ^hist =0 is already minimal. The hyperparameter γ was set to 0 [Figure 7(a)] and 0.5 [Figure 7(b)], respectively. While for both cases the two objects are somewhat recovered, many of the initial single pixels at wrong positions remain, even when the histogram prior is switched on [Figure 7(b)]. Overall the image quality is much poorer than starting from a homogenous initial guess. In general, we found that starting with a homogeneous guess provides better images than starting with a wrong heterogeneous initial guess.

Figure 7

Image reconstruction starting from nonhomogeneous initial guess. Histogram of initial guess equals that of the target medium, however pixels with different D are randomly placed: (a) reconstruction without penalty term; (b) result with histogram penalty term (γ=0.5).

While the presented examples only show one type of medium, similar results have been obtained for a variety of differently structured media. The main results, namely that additional a priori information can be expressed in appropriate penalty terms, and that more information leads to better reconstructions and less sensitivity to a “correct” hyperparameter γ, is of general validity. Another example is shown in Figure 8 where instead of two compact elliptical inhomogeneities two curved lines with different diffusion coefficient cross a 4×4 cm medium. The medium is discretized into a 40×40 mesh with Δx=0.1 cm. The optical properties of the background are given by μ_a=0.1 cm ⁻¹, μ_S ^′=8 cm ⁻¹, (⇒D=0.905 cm ² ns ⁻¹), while the optical properties of the inclusions are given by μ_a=0.2 cm ⁻¹, μ_S ^′=13.0 cm ⁻¹ (⇒D=0.556 cm ² ns ⁻¹), and μ_a=0.05 cm ⁻¹, μ_S ^′=4.5 cm ⁻¹ (⇒D=1.61 cm ² ns ⁻¹), respectively. The initial guess for all reconstructions is D=1.0 cm ² ns ⁻¹. Figure 8(a) shows the original image and Figures 8(b), 8(c), and 8(d) provide the reconstruction results after 40 iterations without a penalty term, with tissue-type penalty term, and histogram penalty term, respectively. Again we can see how penalty terms improve the image quality and that the use of more information leads to better results. The tissue-type penalty term only leads to improved images for γ<1.0, while the histogram penalty term can be applied with a much wider range of γ values. What type of penalty functions and choices of hyperparameters will be most suitable for various other applications such as optical tomographic imaging of the brain, breast, limb, joint, and other body parts remains to be determined for each case. Here we have limited ourselves to the derivation of the appropriate framework for such studies.

Figure 8

Example of reconstruction of two curved lines after 40 iterations. The initial guess is a homogeneous medium with D=1 cm ² ns ⁻¹: (a) original medium, (b) reconstruction without penalty term (γ=0), (c) reconstruction with tissue-type penalty term [Eq. (8), γ=0.8)], and (d) reconstruction with histogram penalty term [Eq. (12), γ=0.8].

Finally, it is interesting to compare the penalty term approach of gradient based schemes, with regularization schemes employed in the linear-perturbation approach (see Sec. I). It can be shown that the minimization of the functional in Eq. (7), which contains the χ² term as well as a penalty term ω⋅Π, is equivalent to solving a matrix equation for Δζ of the form ^{17 30 32}

For the tissue-type penalty term we obtain

For the histogram penalty function we obtain:

Therefore, (R ^hist )_ij has nondiagonal entries, which means that coupling between pixels occurs. However, calculating R ^hist for the examples used in this work, we find that the values of diagonal elements are in general 3 orders of magnitude higher than the values of nondiagonal elements. Therefore, adding either the tissue-type or histogram penalty function to the χ² term in a GIIR scheme, can also be interpreted as adding diagonally dominant matrices to the ill-conditioned matrix JJ ^T obtained in the linear perturbation approach to the image reconstruction problem.

The advantage of the penalty-term approach in GIIR schemes over regularization schemes employed in linear perturbation methods seems to be that penalty terms have an immediate physical interpretation. Even though useful regularizers can obviously be derived without the concept of penalty terms as shown for example by Pogue et al.;, ³² the concept of penalty terms may provide an intuitive conduit through which we can incorporate additional information into the reconstruction process. Since penalty terms can be considered as potentials that bias the χ² potentials according to our prior knowledge about the system they may supply more insight into the underlying physical assumptions about the system to be reconstructed. Furthermore, matrix diagonality as demanded in the linear-perturbation approach to regularization, is clearly not the only factor to get reasonable reconstruction results. As shown in this work, the choice of a good hyperparameter is very critical, and much more work needs to be done to identify robust schemes for finding appropriate hyperparameters.

5.

Summary

In this study we addressed the problem of ill-posedness of the image reconstruction problem in optical tomography. Ill-posedness is caused by the fact that different spatial distributions of optical properties inside the medium can lead to similar detector readings on the surface of the medium. The question arises as to how can one nevertheless find the correct distribution of optical properties inside the medium from measurements on the surface of the medium. In this work we approach the problem within a gradient-based iterative image reconstruction (GIIR) scheme. Using this scheme the image reconstruction problem is interpreted as an optimization problem in which an appropriately defined scalar objective function is minimized using the gradient of the objective function with respect to the distribution of optical properties. In this context ill-posedness is reflected by the absence of a unique global minimum. To overcome this problem we added penalty terms that are derived from a priori information about the system to the commonly used least-square-error term, which compares numerically predicted and actual detector readings.

Specifically we derived and tested so-called tissue-type-penalty terms and histogram-penalty terms. The tissue-type penalty term is derived from the knowledge that a fixed number of different tissues are present in the medium. The optical properties of these tissues are known within certain error margins. If a pixel in the reconstructed image takes on a value that is different from the expected value, a penalty is added to the objective function. The larger the difference the larger is the penalty. The histogram penalty term is derived from the knowledge of the exact histogram of the correct image. Therefore, in addition to the knowledge about the number of different tissue types, one also knows their respective volume fractions in the image.

Performing numerical studies we found that both types of penalty terms can improve the image quality as well as the convergence rate of the iterative image reconstruction process. However, a crucial variable is the coupling or hyperparameter that adjusts the relative strength of the penalty term with respect to the least-square-error term. The tissue-type penalty term only provides image improvement for a narrow range of hyperparameters. If the hyperparameter is not chosen correctly the addition of this penalty term can have a detrimental effect on the image quality. On the other hand, the histogram-penalty term is much less sensitive to the right choice of the hyperparameter and in almost all cases leads to improved reconstruction results.

Finally, we compared the penalty terms derived for the gradient-based iterative image reconstruction scheme, with regularization schemes employed in the linear-perturbation approach. We showed that both tissue-type and histogram-penalty terms lead to diagonally dominate matrixes in the linear-perturbation approach.