1. Introduction

Frequency-domain optical tomography (FDOT) is an emerging biomedical imaging modality [1–4]. This method estimates the spatial distribution of optical properties in tissues by analyzing the amplitude and phase shift of transmitted amplitude-modulated light measured at boundary surfaces. State-of-the-art image reconstruction codes employ a frequency-domain forward model of light propagation that leads to predictions of measured values on the boundary, assuming a certain distribution of optical properties inside the medium [1,5–9]. In these so-called model-based iterative image reconstruction (MOBIIR) codes an objective function is defined that quantifies the differences between model predictions and actual measurements. The minimum of this objective function is sought by updating the parameters of the forward model. Light propagation in tissue is typically modeled either by the equation of radiative transfer (ERT) [10] or by its diffusion approximation (DA) [11].

Imaging of small tissue volumes is an area of great interest, mainly fueled by advances in small animal models of human diseases [12,13]. These models allow the study of disease genesis, progression and treatment and have already provides many new insights, especially in the case of cancer [14,15,16,]. Another promising application is imaging of fingers for the diagnosis of joint diseases [17–19].

Employing FDOT for imaging of small tissue volumes still poses a variety of challenges that have not yet been overcome. First of all, most frequency-domain reconstruction codes are based on the diffusion approximation to the equation of radiative transfer. This approximation, however, becomes less accurate when applied to small tissue-geometries and is further compromised if highly absorbing objects or fluid-filled regions, which contain, for example, cerebrospinal or synovial fluids, are considered [11]. Transport-theory-based codes can accurately model these types of tissues, and first algorithms [1,5,6,10] of this kind were recently developed. However, these codes have used the low-order spatial differencing scheme called the step (or upwind) scheme [20], which causes false scattering and can lead to large errors in predicted amplitudes and phase shifts, especially when small tissue geometries are considered. Therefore, in this study we make use of the second-order scheme [21].

But even with this code in place, a practical challenge is to find optimal source-modulation frequencies at which to perform optical tomographic imaging. In small tissue geometries, the phase-shift at low frequencies is typically very small and difficult to measure. Increasing the modulation frequency leads to larger phase shifts, but at the same time the amplitude signal decreases. Toronov et al. [22] and Gu et al. [23] have studied these tradeoffs. Gu et al. employed a transport-theory-based frequency-domain forward model. In numerical studies they have shown that for typical geometries encountered in small animal imaging, highest signal-to-noise-ratios (SNRs) can be achieved in the 400-600 MHz frequency range. What has not been done, however, is to show how these different SNR values in the measured data affect the reconstruction results. In this study we look to find source modulation frequencies that result in best image quality over a broad range of optical properties. We show results of numerical studies and experiments on tissue phantoms with well characterized optical properties that mimic small animals and human fingers.

2. Methods

2.1. Model of light propagation

In frequency-domain imaging systems, the light source is amplitude modulated at frequencies in the 50–1000 MHz range, and the demodulation and phase shift of the resulting photon-density waves are measured. It has been shown that, in general, the quality of frequency-domain reconstructions [1,5] is superior to the steady-state approach.

The frequency-domain forward problem for light propagation in turbid media can be accurately modeled by the frequency-domain equation of radiative transfer (FD-ERT) [1,5], given by

where ψ(r,Ω,ϖ) is complex radiation intensity in unit [Ω/cm²/sr.], µ _a and µ _s are the absorption and scattering coefficients, respectively, in units of [cm^-1], ϖ is the external source modulation frequency and c is the speed of light inside the medium, Φ(Ω⁺,Ω) is the scattering phase function that describes scattering from incoming direction Ω⁺ into scattering direction Ω, and S _m(r,Ω,ϖ) is the source term due to medium emission. In this work we employ the widely used the Henyey-Greenstein phase function [24], which is given by

Furthermore, to be able to consider the refractive index mismatch at air-tissue interface [7,20], we implemented a partially-reflective boundary condition as

where R(Ω′,Ω) is the reflectivity at Fresnel interface from direction Ω′ to direction Ω, ψ ⁰(r _b,Ω,ϖ) is the radiation intensity due to the external source function and subscript b denotes the boundary surface of the medium, while n⃗_b is the unit normal vector pointing outwards the boundary surface.

To find numerical solution for these equations, the spatial and angular variables have to be discretized. In particular, we employ a node-centered finite-volume discretization [21,25] in the spatial domain and use a discrete ordinate approach [20] for the angular domain. The node-centered finite-volume method takes advantage of the beneficial properties of both the finite-element and finite-volume methods, thus combining the conservation properties of the finite-volume formulation and the geometric flexibility of the finite-element approach [21]. When using the node-centered finite-volume discrete-ordinates method, the discretized form of radiative transfer equation is obtained by integrating equation (1) over the control volume with a divergence theorem as

where N _surf and N _Ω are the number of surfaces surrounding the node N and the number of discrete ordinates based on the level symmetric scheme [20], respectively, n⃗_j and ψ ^m _j denote the surface normal vector and the radiation intensity defined on the j-th surface. Also the surface intensity ψ ^m _j is related to the nodal intensity ψ ^m _N by the second-order spatial differencing scheme [21]. After discretizations for all nodes, the resulting algebraic equations can be written as follows

where each line of the matrix A contains the coefficients of the discretized equation at node number N and direction m[27]. The vector u denotes the radiation intensity vector and the vector b is the vector of implementing both the discretized boundary condition and the source term. The sparse matrix given by Eq. (5) contains complex-valued elements since we treat the frequency-domain equation of radiative transfer Eq. (1) directly, instead of separating it into two real-valued equations as found in other works [1,6]. As a result, the complex-valued sparse, linear system of equations given by Eq. (5) is iteratively solved for intensity into a discrete-ordinate direction by using a complex version of the GMRES Krylov-subspace solver [1,10,28]. In this study we employed the lower-upper-symmetric Gauss-Seidel (LUSGS) preconditioning matrix that outperforms the incomplete lower-upper (ILU) matrix [29,30].

The FD-ERT is used as a forward model within a model-based-iterative image reconstruction (MOBIIR). Therefore, the FD-ERT provides the prediction of the detector readings that are to be compared with actual measurements. Next we discuss the inverse model, which is used to obtain the spatial distribution of optical properties that best fit the measured data.

2.2. Inverse model

The inverse problem associated with optical tomography is to find the vector x=(µ ¹ _a,…,µ ^Nt _a,µ ¹ _s,…,µ ^Nt _s) that describes the spatial distribution of the optical properties inside the medium, given the measurement data obtained on the surface of the medium. This begins by defining an objective function that quantifies the mismatch between predicted and measured amplitude and phase shift of the photon density waves that travel through the medium. As mentioned in the introduction, as the source modulation frequency increases, the amplitude signal decreases and the phase shift increases.

This leads to different magnitudes of errors in the amplitude and phase-shift measurements. In other words, the standard deviation is not only position-dependent, but it is also a function of frequency [22,23]. Accordingly the objective function has to be designed so that it can treat position-and frequency-dependent magnitudes of errors with different weights. In this study, we chose the weighted least-square error norm [31] given by

where N _s and N _d are the numbers of sources and detectors used, Q _d is the measurement operator that maps the calculated radiance vector u _s into predictions Q _d u _s of measurements z _s,_d for source-detector pairs (s,d), σ _s,_d denotes the complex standard deviation at (s,d) source-detector pair, and the operator (·)* denotes the complex conjugate of the complex vector. Thus the optimization problem in Eq. (6) is characterized by N _t N _Ω forward variables and N _t(or 2N _t) inverse variables. With the objective function given by Eq. (6), less precise data is given a smaller influence over the solution, whereas precisely measured data is weighted stronger to have a large influence on the solution [31]. To obtain a stable solution, the optimization process should be stopped when the value of f(x;u) is similar to the magnitude of the errors.

Furthermore, in this study the calculation of the gradient vector of the objective function with respect to each of unknown parameters is performed with an adjoint formulation [1,5,6,32]. In the frequency domain this can be formulated as

where A ^T denotes the adjoint operator equivalent to the transpose of the sparse matrix A in Eq. (5). Equation (7) is solved by the iterative GMRES method. The gradients of the objective function is obtained by differentiating f(x;u) given by Eq. (6) with respect to x ⁱ=(µ ⁱ _a,µ ⁱ _s) at the i-th control volume δV ⁱ as

Table 1. Optical properties of the homogeneous numerical phantom for the SNR study.

View Table | View all tables in this article

which represents the gradient vector g _r=(∇_{µ a} f,∇_{µ s} f). Finally the unknown optical properties are iteratively updated by using the limited-memory version of the BFGS optimization method [32].

3. Results

3.1. Numerical studies

3.1.1. Analysis of SNR values using numerical phantoms

Before we show results on how amplitude and phase noise influence the tomographic reconstruction results, we generalize some of the findings reported by Toronov et al. and Gu et al.. This will help us later in interpreting the results found in the reconstruction studies.

The noise model introduced by Toronov et al. [22] and modified by Gu et al. [23], is given by:

where σ _AC and σ _Φ are the standard deviations of the amplitude and the phase shift respectively, and 〈·〉 represents the ensemble average of the corresponding quantity. σ _AC and σ _Φ are proportional to $\sqrt{\mathrm{DC}}$ and $\sqrt{\mathrm{DC}}/\mathrm{AC}$, respectively, i.e., the phase noise increases with the frequency. This model is used to assess the SNR for a geometry specific to our study. In previous publications it has already been shown that the optical properties of the background medium have a strong influence on the SNR value, whereas small perturbations in either absorption or scattering, have little affect on the SNR, unless their volume is larger than that occupied by the background medium.

We start by studying the frequency dependence of the SNR value by changing the intrinsic optical properties of a homogenous medium for a specific geometry. To this end we consider a 3-cm by 3-cm square phantom with a homogenous distribution of optical properties in the medium as shown in Fig. 1. The synthetic data are generated by solving the frequency-domain transport Eq. (1) and corrupting the result by Gaussian random noise according to the noise model given in Eq. (10). We have obtained the SNR values at specified locations (see Fig. 1) from running forward simulations in the range of 0.05 cm^-1≤µ _a ^≤0.5 cm^-1 for the absorption coefficient and in the range of 5.0 cm^-1≤µ ^′ _s≤20 cm^-1) for the reduced scattering coefficient (see Table 1). The SNR values for the amplitude and the phase shift are shown in Figs. 2 and 3. Note that all the values are normalized to the SNR values at 100 MHz since the proportionality constants in Eq. (10) are not explicitly known.

 

Fig. 1. Schematic of the numerical phantom used in this section. The source location is indicated by the read arrow, and the detector positions are given by the black arrows.

Download Full Size | PDF

 

Fig. 2. Normalized SNR values for the AC signal as a function of source-modulation frequency for six different source-detector pairs (see Fig. 1 for reference). The 4 graphs show 4 sets of optical properties. The absorption coefficient µ _a increases with a column, while the reduced scattering coefficient µ ^′ _s increases with a row.

Download Full Size | PDF

Comparing Figs. 2(a) and 2(b), we observe that keeping µ _a constant and increasing µ ^′ _s, the frequency dependence of the SNR becomes more pronounced. On the other hand, when we keep µ ^′ _s constant and increase µ _a from 0.05 to 0.5 cm^-1, the frequency dependence of the amplitude SNR becomes less pronounced [see Figs. 2(a) and 2(b)]. This behavior can be understood by introducing the formal solution [20] of the frequency-domain transport equation at specified angular direction Ω

Here the first term is the contribution from the boundary on the opposite side, exponentially attenuated through the distance from 0 to s. Thus the first term is always zero unless it has an external source on the opposite side. The overall source term Q _s(s ^′) is defined as

where the first term is the medium emission (e.g., fluorescence) and the second term is the in-scattering term. Thus the integrand of the second term in Eq. (11) is the contribution from the source term Q _a(s ^′) at position s ^′ due to the medium emission as well as in-scattering, exponentially attenuated through the distance from s ^′ to s. It is now obvious from the second term of Eq. (11) that the change in µ _a has a stronger impact on the amplitude signal than the change in µ _s. In other words, the increased µ _a augments the exponential attenuation term, which makes the medium’s absorption more dominant than the frequency dependent absorption. That is why the curves on Figs. 2(b) and 2(d) look somewhat frequency-independent.

 

Fig. 3. Normalized SNR of phase shifts at various frequencies. The absorption coefficient µ _a goes higher in the column direction, while the reduced scattering coefficient µ ^′ _s goes higher in the row direction.

Download Full Size | PDF

On the other hand the effect of scattering is more complicated since µ_s involves both the source term Q(s ^′) and the exponential decay exp[-∫^s _s ^′(µ _a+µ _s+iϖ/c)ds] in Eq. (11). Accordingly an increase in µ _s augments the exponential decay but, at the same time, increasingly dominates the source term due to in-scattering. As a consequence, the effect produced by a change in µ ^′ _s is moderately compromised according to the degree of scattering. This can be clearly seen from the first and second columns in Fig. 3: the change in frequency dependence is only minor for both cases.

The frequency dependence of the phase SNR value is much easier to understand. Looking at the noise model in Eq. (10) one can see that the phase SNR equals the amplitude SNR multiplied by a phase shift. As the source modulation frequency increases the amplitude SNR decreases, while the phase shift increases. Therefore one can assumed that there will be a certain frequency for which the phase SNR takes on a maximum value. And indeed, Figs. 3a and 3c show that there exist a modulation frequency between 400 MHz and 800 MHz for which the phase SNR is largest at all positions except for locations 1 and 2. We also found that the position of this maximum moves to higher frequencies when µ _a is increased [Figs. 3(b) vs. 3(a), Figs. 3(d) vs. 3(c)] and µ ^′ _s is decreased [Figs. 3(c) vs. 3(a)]. Another finding is that the SNR of the phase shift improves with frequency regardless of µ ^′ _s when µ _a is increased [compare Figs. 3(a) vs. 3(b), and 3(c) vs. 3(d) for 200 MHz and 1 GHz].

Overall we observe that the SNR value of amplitude and phase shift depend on the source modulation frequency, the optical properties of the medium and the position of the detectors for a given source. For the cases considered in this study, which are typical for small animal and finger joint imaging, the optimal source modulation frequency with respect to largest SNR values seem to lie between 400–800MHz for most source-detector pairs and optical properties. In the following section we will study the influence of source-modulation frequencies and measurement SNR values on the reconstruction results. We will start by considering synthetic data and later use experimental data to test our theoretical findings.

3.1.2. Influence of SNR values on tomographic reconstruction results

In this section we present the reconstruction results using the noise-added synthetic data, which is based on the noise model given by Eq. (10). To this end, we consider a 3cm-by-3cm square phantom with three inclusions shown in Fig. 4(a). The first two inclusions are highly scattering while the other object is highly absorbing. The optical properties of the background medium are µ _a=0.62 cm^-1 and µ ^′ _s=7.62 cm^-1, while the two highly scattering inclusions are given µ _a=0.1 cm^-1 and µ ^′ _s=25.7 cm^-1 and the purely absorbing object is given µ _a=1.2 cm^-1 and µ ^′ _s=7.62 cm^-1. We first examined the frequency dependence of the SNR value by using the forward data. The results are shown in Figs. 4(b)–4(c). The amplitude SNR decreases only moderately (approximately 20%) and the phase SNR increases with frequency in all positions. This behavior is similar to the amplitude and phase shift SNR curves as shown in Figs. 2(b) and 2(d), and Figs. 3(b) and 3(d). This implies that even at higher frequencies the loss of amplitude information is only small and the quality of the reconstruction will improve by increasing the frequency.

Next we performed reconstructions at 5 different source modulation frequencies (0, 200, 400, 600, 800 MHz). As before, the noise-added synthetic data was generated by applying the noise model in Eq. (10). Furthermore, we normalized both amplitude and phase independently for measurements on each side of the square phantom, since the experimental data as will be shown later, was obtained in this way. The 4 sources were placed around the phantom surface with the 29 detectors, and all sources were amplitude-modulated. The same setup and same procedure were repeated for the 5 modulation frequencies. The reconstruction procedure was started with an initial guess that was identical to the optical properties of the background medium.

To quantify the quality of reconstructed images, we used the correlation factor ρ(µ ^e,µ ^r) and the deviation factor δ(µ ^e,µ ^r)erdmm as introduced in reference [7]:

Here µ̄ and σ(µ) are the mean value and the standard deviation for the spatial function of the optical property that can be either the absorption coefficient µ _a or the µ ^′ _s. Similarly, µ ^e and µ ^r are the exact and reconstructed distributions of optical properties, respectively. In terms of quality of the reconstruction results, the correlation coefficient indicates the degree of correlation between exact and estimated quantities while the deviation factor describes the discrepancy in absolute values of exact and estimated quantities. Accordingly, the closer ρ(µ ^e,µ ^r) gets to 1, and the closer δ(µ ^e,µ ^r) gets to 0, the better is quality of the reconstruction.

 

Fig. 4. Problem set up and normalized SNR values; (a) schematic of a square phantom; (b) amplitude SNR values; (c) phase SNR values.

Download Full Size | PDF

 

Fig. 5. Reconstructed images of µ _a and µ ^′ _s for 5 different source-modulation frequencies using a numerical phantom. Noise was added to the synthetic data in accordance with Eq. (10).

Download Full Size | PDF

The reconstructed images are shown in Fig. 5. Note that data points within 2mm from the boundary are excluded in the calculation of ρ(µ ^e,µ ^r) and δ(µ ^e,µ ^r) in order to avoid simulation artifacts near the boundary surface. As can be seen in Fig. 5 there is considerable cross-talk between scattering and absorbing objects in the DC case (f=0 MHz). Therefore, the two scattering inhomogeneities appear in the absorption maps and the absorbing object appears in the scattering reconstruction. As the source-modulation frequency is increased the cross-talk is greatly reduced and the reconstruction separates the absorbing and scattering objects.

Figure 6(a) shows that the correlation factor ρ(µ ^e,µ ^r) is largest at 800 MHz and smallest at 0 MHz, which is in good agreement with the results shown in Fig. 5. Note that the scattering images are more accurate than the absorption images for all frequencies when the correlation coefficient is considered. The deviation factor δ(µ ^e,µ ^r) of the scattering images gradually decreases with frequency and yields the smallest value at ~800 MHz. This is to be expected because the phase SNR improves steadily with frequency and this improvement leads to better reconstructions of µ ^′ _s in the image.

 

Fig. 6. Frequency dependence of correlation (ρ) and deviation (δ) factors for the reconstructed images shown in Fig. 5.

Download Full Size | PDF

It is well-known that the amplitude signal has more impact on the reconstruction accuracy of µ _a. Hence at high frequencies one will lose more information about µ _a. However for this particular case our observation shows that the image accuracy of µ _a improves steadily with increasing frequency. One explanation for this apparent paradox could be that the errors in amplitude and phase shift differ substantially. In other words, the amplitude SNR is decreased by only 20% at 800 MHz, whereas the phase SNR is improved by 300% at 800 MHz (see Fig. 4). Thus the increase in phase SNR is large enough to compensate for the loss in the amplitude SNR; this leads to an overall improvement of the reconstruction results with increasing frequency. As shown in Fig. 6, the deviation factor δ(µ _e,µ ^r) decreases with frequency and reaches the smallest value (i.e., most accurate) at 600–800 MHz. It is also notable that the curves of δ(µ ^e,µ ^r) follows the reciprocal property of ρ(µ ^e,µ ^r); δ(µ ^e,µ ^r) is largest when ρ(µ ^e,µ ^r) is smallest and vice versa.

3.2. Application to experimental data from general tissue phantoms

With numerical results at hand, we studied the influence of the frequency-dependent SNR on the quality of reconstructions performed with experimental data. For this purpose we reconstructed the frequency-domain data obtained for two small-volume phantoms whose optical properties are nearly identical to those used in the numerical studies. We first give a brief description of our frequency-domain instrumentation and the phantom used.

3.2.1. Frequency-domain imaging system

The frequency-domain system [33] is designed for fast two-dimensional imaging of modulated light diffusely transmitted through small-tissue volume. As shown in Fig. 7, the main components of the system are the illumination part, the detection system, and the modulation sources for the light source and detector. The master signal generator provides a sinusoidal AC input to the laser diode driver that supplies the laser diode (wavelength λ=670 nm) with a bias and AC current. The laser illuminates the surface of the object. The position of the laser spot is adjusted with translation stages. The modulated light transmitted through the object is imaged by a lens to an intensified CCD (ICCD) camera.

The system operates in homodyne mode, i.e. the gain of the ICCD is modulated by a slave signal generator at the same frequency as the laser. As a result a steady state image at the intensifier output is imaged to the CCD. The signal in every pixel depends on the phase between source and detector modulation. Master and slave signal generators are linked together and the phase delay is adjustable. To detect the complete oscillation of the modulation multiple images are taken at phase delays covering the range of 2π and are transferred to a computer.

 

Fig. 7. Experimental setup of the frequency-domain system; (a) photograph; (b) schematic: 1. Laser diode, 2. Laser diode drive, 3. & 7. Signal generator, 4. Sample, 5. Lens, 6. ICCD camera, 8. High rate imager, 9. Computer.

Download Full Size | PDF

 

Fig. 8. Schematic of an experimental square phantom; (a) dimensions; (b) & (c) camera orientation for measurments. The CCD camera is fixed to the left side from the phantom. The turning arrow indicates that measurements are taken counterclockwise over the three sides.

Download Full Size | PDF

The system is controlled via a graphical user interface. From the stack of images two-dimensional amplitude and phase images are calculated by fast Fourier transformation (FFT) in every pixel. More details concerning this setup can be found in the reference [33].

3.2.2. Square phantom

The phantom with a square base has three holes along the z-direction. The sides are numerated as I, II, III and IV and their dimensions are shown in Fig. 8(a). The holes are filled with whole milk (fat content 3.5%) as a purely scattering perturbation or with Evans Blue solution as purely absorbing perturbation.

The measurements were made with a setup in which the laser and the camera are arranged at a right-angled as shown in Fig. 8(b). The phantom is illuminated at the center of each side. The distance from the object plane to the lens (50 mm, set to f/2) is approximately 20 cm, which forms the aperture angle of 7 degree when viewed from the optical components. On the CCD, an 8x8 binning is used that yields an image scale at which 1 mm equals 2.3 pixels. To image one side of the phantom with right-angled illumination, the phantom, mirror, and camera are moved so that the distance between phantom and camera is kept constant and the laser beam hits the phantom surface directly and not via the mirror [see for example Fig. 8(b), where the source is located side I and measurements are taken on side IV]. To image the phantom sides with illumination at the opposite side, e.g., illumination at side II and detection at side IV, the laser beam is redirected by the mirror [Fig. 8(c)]. At right-angled illumination we inserted a neutral density filter (NG9, D1.6, transmission T=0.011) because the detected intensity was much higher than during the illumination at the opposite side. We used the three sides for the measurements, for example, with illumination at side I and measurements taken on sides II, III, IV. The illuminated side was excluded from measurements and reconstructions.

 

Fig. 9. Reconstruction images of µ _a and µ _s for a tissue phantom, with well-known optical properties. The source modulation frequency was varied from 0 to 800 MHz.

Download Full Size | PDF

 

Fig. 10. Correlation (r) and deviation (d) factors for the reconstructed images shown in Fig. 9.

Download Full Size | PDF

All images were calibrated with respect to photocathode non-uniformity and to amplitude and phase of the laser source, measured by introducing neutral density filters. Therefore, the amplitude is relative to the filter transmission. The different positions of the phantom for opposite and right-angled illumination cause a difference in the optical path length and therefore a phase delay. We measured the phase delay for this path and excluded this phase delay in the calibration procedure. A phase delay can be quite accurately calibrated, but the amplitude signal is still relative to its absolute value because of some unknown filter effects. Therefore, for the image reconstruction we chose to use the normalized data that can eliminate this ambiguity. We obtained the measurement data at four frequencies: 0 (steady state), 200, 400, 600 and 800 MHz.

The reconstruction results are shown in Fig. 9. As already observed in the numerical studies, the DC data reveals strong crosstalk between scattering and absorption inhomogeneities. Using data obtained at higher source modulation frequencies, the cross-talk can be greatly diminished. The accuracy measures of these images in terms of correlation and deviation factors are also similar to those obtained from numerical studies (see Fig. 10). For the µ _a reconstruction, the correlation factor ρ(µ ^e,µ ^r) increases with increasing frequency whereas the deviation factor δ(µ ^e,µ ^r) decreases with frequency. The µ _s reconstructions show a maximum ρ(µ ^e,µ ^r) and a minimum δ(µ ^e,µ ^r) at 600 MHz whereas µ _a reconstructions show largest correlation rho and smallest deviation delta at 800 MHz.

3.3. Finger phantom

To extend our analysis to some practical applications, we furthermore considered measurements on a tissue phantom that mimics a human finger. It has been shown that

Table 2. Optical properties of the finger phantom shown in Fig. 11. The two sets of properties mimic the fingers of a normal healthy person and a patient affected by rheumatoid arthritis.

View Table | View all tables in this article

measurements of light transmission through finger joints can be used to obtain information on the physiologic state of the joint [17,18]. However all studies performed so far were limited to continuous wave data. As we have just shown, in this case, significant cross-talk between µ _a and µ ^′ _s can occur. In this study we explore the hypothesis that better images of the finger joint can be obtained using frequency-domain data.

 

Fig. 11. Structure of a finger phantom; (a) layout of major parts; (b) schematic of source-detector configuration.

Download Full Size | PDF

To this end we consider a finger phantom of well-characterized optical properties as introduced by Netz et al. [33], which is shown in Fig. 11. The optical properties of this phantom can be varied and were adjusted to mimic the fingers of a healthy person and a patient with rheumatoid arthritis (RA) (see Table 2). More details on this phantom can be found in the reference [33].

The experimental setup for this study is identical to the previously described system (see Fig. 7), except that the square phantom is replaced here by the finger phantom. For this experiment we illuminated the phantom with the laser source at 11 different locations, and measured the transmitted light intensity with the CCD camera [Fig. 11(b)] at 5 different source-modulation frequencies (0, 200, 400, 600 and 800 MHz). At each frequency we scanned the phantom several times to improve the SNR of the measured data. In concordance with our previous work [17,18], 31 detector readings were taken on the opposite surface of the side that was illuminated with the laser.

We start our analysis by comparing the difference in amplitude and phase shifts obtained at 600 MHz for the phantom representing a healthy person and a patient affected by RA. As can be shown in Fig. 12(a), the amplitude signal gradually increases towards the center of the finger, where the joint cavity is located. This increase is more pronounced in the healthy case when compared to the RA case. This can be explained by the increase of both µ _a and µ _s in the joint fluid and capsule. As for the phase delay, the RA phantom yields larger phase delays than the healthy phantom, which can be understood by the larger scattering coefficient of the joint fluid in the RA case.

 

Fig. 12. Profiles of amplitude and phase shift obtained along the line-of-measurement at 600MHz. The laser source is located at the joint center, on the opposite side of the phantom (see Fib. 11b).

Download Full Size | PDF

 

Fig. 13. SNR values of amplitude and phase shift obtained for the real finger phantom; (a) no RA; (b) RA

Download Full Size | PDF

We also examined the SNR values for the two cases. Figure 13 shows the SNR curves of amplitude and phase shifts obtained for a 10×10 pixel area opposite to the laser spot at the joint gap. It can be seen that the amplitude SNR weakly depends on the frequency, i.e., it decreases by only 2% per 100 MHz, whereas the phase SNR shows a strong dependence. It increases from 10 at 100 MHz to 80 at 800 MHz in the healthy case and 70 at 900 MHz in the RA case. This behavior is consistent with the numerical studies of similar optical properties shown in Figs. 2 and 3.

Three dimensional image reconstructions were performed using a computation domain that was discretized with 10234 tetrahedron elements as shown in Fig. 14. For the reconstructions we only used data from the center line of the images captured by the CCD camera. The CCD image represents the projected image of the outgoing light distribution from the object’s surface onto the virtual imaging plane of the CCD camera. Thus it requires the knowledge of a projection operator to make use of the entire image for the detection geometry. However, this projection operator was not available at the time of measurement. All reconstructions were started with an initial estimate of µ _a=0.4 cm^-1 and µ ^′ _s=10 cm^-1.

 

Fig. 14. The schematic of the finger phantom: cylinder height H=3.2 cm, diameter D=2 cm; (a) source-detector configuration; (b) computation domain with 11023 tetrahedrons and simplified cross-section view of the orientation of the finger joint phantom (red: bone, blue: synovial fluid, green: skin tissue).

Download Full Size | PDF

 

Fig. 15. Example of reconstructed distributions of µ _a values for the healthy (left) and RA (right) finger phantom. In this case the source-modulation frequency was 600 MHz.

Download Full Size | PDF

 

Fig. 16. Interactive images using VolView software, obtained for Fig. 15. (a) Healthy (View 1). (b) Affected (View 2).

Download Full Size | PDF

Results of the reconstructions are shown in Figs. 15 through 18. Figures 16 and 18 contain complete 3D data sets that can be viewed using the interactive VolView software; while Figs. 15 and 17 show representative 2-D cross-sections extracted from the 3D data sets. Looking at Figs. 15 and 16 one notices that the spatial distribution of optical properties varies much stronger in the healthy case, when compared to the results obtained for the finger phantom that mimics RA. This can be quantified by computing the minimum and maximum values of µ _a and µ ^′ _s. Table 3 summarizes these results for 5 different source-modulation frequencies. Evaluating the table one can see that in most cases the contrast (defined as ratio of maximum and minimum of µ _a and µ ^′ _s respectively) between the joint center and the other areas is larger in the healthy joint, except at 0 and 200 MHz where the contrast in the scattering images is slightly larger in the RA case. Furthermore we observe that in most cases the contrast appears to be strongest at 0 MHz. The reason for that is not quite clear at this point, but may be related to the fact that we use the relative rather than absolute data.

Table 3. Min. and max. values of reconstructed µ_a and µ^′_s for the ROI that covers the joint center and the bone area while excluding 2mm from the boundary. The values in parentheses represent the ratios (max./min.).

View Table | View all tables in this article

 

Fig. 17. Reconstruction images of µ ^′ _s obtained for the healthy finger-joint phantom at different frequencies. Note that the cross-talk in the bone is greatly diminished with increasing frequency. Note that the cross-talk in the bone is greatly diminished with increasing frequency.

Download Full Size | PDF

 

Fig. 18. Interactive images using VolView software, obtained for Fig. 17. The images are shown here for (a) 0 MHz (View 3) and (b) 600 MHz (View 4) only.

Download Full Size | PDF

The loss in contrast at higher modulation frequencies is offset by the reduction in cross-talk between absorption and scattering effects at higher frequencies (see Fig. 17 and 18). We found that the reconstructed images generated with 600 or 800 MHz source-modulation frequencies more closely reflect the true distribution of optical properties given in Table 2. For example, in the case of the healthy phantom, the reconstructed-µ _a in the bone area is approximately 1.4~1.7 cm^-1 at lower frequencies, and decreases to the more accurate value of ~0.9 cm^-1 at higher frequencies (see Fig. 15, 16, and Table 3). Furthermore, according to Table 2, µ ^′ _s’s of the bone and the skin-tissue complex are the same. However, in the reconstructed µ ^′ _s-image the lower frequency data reveals absorbing effect of the bone, while this effect is greatly diminished at higher frequencies (see Figs. 17 and 18). All these results comply well with the SNR characteristics shown in Fig. 13.

3.4 Small animal studies

In addition to the forgoing applications, we performed studies to explore how the choice of different source-modulation frequencies affects reconstruction results in small animal imaging. For these studies we use an anatomically accurate 3D model of a mouse shown in Fig. 19. The model was created from 90 axial magnetic resonance images with a slice separation of 0.6 mm. Each image was segmented into a whole mouse section by thresholding in combination with some minimal manual smoothing. Furthermore, to perform numerical simulation that mimic the “real world” as close as possible, we used the noise-characteristics of the measurement system described in section 3.2.1 [33]. With the mouse model and the noise characteristic at hand, we simulated measurements at 5 source-modulation frequencies (0, 200, 400, 600, 800 MHz) at several cross sections including the brain, lung, liver, and kidneys. These simulated measurements were input to our frequency-domain code and the reconstruction results were compared to the known target medium. As measures of the accuracy we calculated again the correlation factor ρ(µ ^e,µ ^r) and deviation factor δ(µ ^e,µ ^r) as defined in Eq. (15).

A reconstruction example is shown in Fig. 20. This figure shows a cross section through the model of the mouse’s head and the reconstruction results for µ _a (top row) and µ ^′ _s (bottom row) at source modulation frequencies of f=0 MHz (steady state) and 600 MHz. Clearly visible is the much higher accuracy of the 600 MHz reconstructions. At f=0 MHz, a strong cross-talk between absorption and scattering effects can be observed, with the µ _a reconstruction being too low and the µ ^′ _s reconstruction too high. Figure 21 shows the frequency-dependent values of ρ(µ ^e,µ ^r) and δ(µ ^e,µ ^r) respectively for both µ _a and µ ^′ _s. It can be seen that the correlation coefficient ρ(µ ^e,µ ^r) takes on a maximal value at 600 MHz for both µ _a and µ ^′ _s reconstruction. The deviation factor δ(µ ^e,µ ^r) is smallest at 600 MHz for the µ ^′ _s reconstruction, while for the µ _a reconstruction is smallest at 800 MHz.

 

Fig. 19. Numerical model of a mouse anatomy. Black arrows indicate cross-section for which reconstructions results were produced.

Download Full Size | PDF

 

Fig. 20. Cross-sectional images of µ _a and µ ^′ _s in the brain of the virtual mouse, obtained with source-modulation frequencies of f=0 MHz (steady state) and f=600 MHz.

Download Full Size | PDF

Similarly Figs. 22 and 23, show the frequency dependence of ρ(µ ^e,µ ^r) and δ(µ ^e,µ ^r) for cross section through the lungs and liver. It can be seen that dependent on the particular cross section and dependent on the optical properties (µ _a or µ ^′ _s reconstruction), ρ(µ ^e,µ ^r) and δ(µ ^e,µ ^r) take on optimal values (largest in the case of ρ(µ ^e,µ ^r) and smallest in the case of δ(µ ^e,µ ^r) at different frequencies. For example µ _a reconstructions in the lung (Fig. 22) show a weak maximum of the value between 400–600 MHz. The deviation factor δ(µ ^e,µ ^r) has a clear minimum at f=600 MHz. Interestingly the µ _a reconstruction yields a δ(µ ^e,µ ^r) that is worse (higher) at 200 MHz than in the steady state case (0 MHz). However, at 600 MHz, δ(µ ^e,µ ^r) is clearly smaller better than at 0 MHz.

Overall these studies show that a frequency-domain system provides superior reconstruction results for the tissue intrinsic background optical properties compared to a steady-state.

4. Conclusion

In this work we have studied the frequency dependence of amplitude and phase shift on optical properties typically encountered in small-tissue volumes. Employing an optical tomographic image reconstructions code that is based on the frequency-domain equation of radiative transfer, we have performed numerical and experimental studies to determine a range of frequencies for which the reconstruction results are best.

 

Fig. 21. Correlation (ρ)(left) and deviation (δ)(right) factors for the reconstructed µ _a and µ ^′ _s images in a cross-section that includes the brain of the virtual mouse (see Figs. 19 and 20), obtained with the noise-added synthetic data at 0 (steady-state), 200, 400, 600, and 800 MHz.

Download Full Size | PDF

 

Fig. 22. Same as Fig. 21, however for images of a cross-section that includes the lungs of the virtual mouse in Fig. 19.

Download Full Size | PDF

 

Fig. 23. Same as Fig. 22, however for images of a cross-section that includes the liver of the virtual mouse in Fig. 19.

Download Full Size | PDF

First we examined the signal-to-noise ratios (SNR) of measurement data using the transport theory forward code and altering the optical properties of the medium. We found that an increase of µ _s in the medium makes the frequency dependence more dominant, whereas an increase of µ _a weakens the frequency dependence of the amplitude SNR value. We also showed that the source modulation frequency for which the highest SNR values are found increases when µ _a is increased and µ _s is decreased. Next we performed image reconstructions using numerical and experimental data obtained with two different two tissue phantoms: a simple square phantom with 3 inhomogeneities and an anatomically correct finger phantom. Furthermore we use numerical studies involving an anatomically correct mouse model to determine optimal frequencies for small-animal imaging. In general we found that best image qualities can be achieved if source-modulation frequencies in the 400–800 MHz range are used.