Two-step reconstruction method using global optimization and conjugate gradient for ultrasound-guided diffuse optical tomography

Behnoosh Tavakoli; Quing Zhu

doi:10.1117/1.JBO.18.1.016006

7 January 2013 Two-step reconstruction method using global optimization and conjugate gradient for ultrasound-guided diffuse optical tomography

Behnoosh Tavakoli, Quing Zhu

Author Affiliations +

Journal of Biomedical Optics, Vol. 18, Issue 1, 016006 (January 2013). https://doi.org/10.1117/1.JBO.18.1.016006

Abstract

Ultrasound-guided diffuse optical tomography (DOT) is a promising method for characterizing malignant and benign lesions in the female breast. We introduce a new two-step algorithm for DOT inversion in which the optical parameters are estimated with the global optimization method, genetic algorithm. The estimation result is applied as an initial guess to the conjugate gradient (CG) optimization method to obtain the absorption and scattering distributions simultaneously. Simulations and phantom experiments have shown that the maximum absorption and reduced scattering coefficients are reconstructed with less than 10% and 25% errors, respectively. This is in contrast with the CG method alone, which generates about 20% error for the absorption coefficient and does not accurately recover the scattering distribution. A new measure of scattering contrast has been introduced to characterize benign and malignant breast lesions. The results of 16 clinical cases reconstructed with the two-step method demonstrates that, on average, the absorption coefficient and scattering contrast of malignant lesions are about 1.8 and 3.32 times higher than the benign cases, respectively.

1. Introduction

Diffuse optical tomography (DOT) is a noninvasive functional imaging modality in which tissue is illuminated by near-infrared (NIR) diffused light, and the reflected or transmitted light is measured at the tissue surface. The interior optical properties of the tissue are then estimated from these measurements by image reconstruction algorithms. DOT has demonstrated clinical potential for distinguishing benign from malignant breast lesions.¹^–¹⁵ Two of the estimated optical parameters used to characterize malignant and benign breast tissues are the absorption and scattering coefficients. By estimating these parameters that are obtained from multiwavelength imaging systems, it is possible to quantify tissue constituents such as oxygenated and deoxygenated hemoglobin, lipid and water concentrations, as well as structural parameters, which are linked to the dimensions and densities of the scattering centers. The complex blood vessel network found in malignant tumors provides the target absorption contrast. Clinical studies of female breasts imaged by different optical tomography systems have demonstrated increased total hemoglobin concentration in the malignant lesions.¹^–⁴ Scattering contrast, on the other hand, is caused by changes in the organelle population within tumor beds. Based on clinical data, a significant contrast ratio in the mean cellular size and the volume fraction between malignant and benign lesions has been found.⁵^,⁶ Accordingly, malignant tissues have shown increased mean values of the reduced scattering coefficients in comparison with the benign tissue.⁷^–¹⁰ Simultaneous reconstruction of absorption and scattering distributions of human breast has also been investigated. It was shown that cysts in the breast can be differentiated from solid tumors because cysts generally demonstrate lower absorption and scattering coefficients compared with the surrounding normal tissue, whereas solid tumors show both higher absorption and scattering relative to normal tissue.¹¹^,¹² In another study, the reconstructed distributions of hemoglobin, oxygen saturation, water fraction, and subcellular scattering in breast tissue showed that the fractions of both blood and water were higher in fibroglandular than in adipose tissue.¹³ In some studies, a synthetic optical index was derived from these intrinsic physiological properties and yielded an average of twofold contrast difference between malignant tumors and normal tissue.¹⁴^,¹⁵

Due to the intense light scattering in biology tissue, the inverse problem of DOT is ill-posed. In addition, the reconstructed optical properties are not unique unless constrained by additional prior information. Furthermore, the inverse problem is underdetermined because the number of measurements is one or two orders of magnitude smaller than the number of voxels with unknown absorption and scattering coefficients.¹⁶ A variety of regularization techniques such as Tikhonov regularization, spatially varying regularization, and Levenberg-Marquardt regularization have been applied to obtain a stable solution for the DOT inverse problem.¹⁷^,¹⁸ Different optimization methods such as the modified Newton-Raphson method, and the truncated Newton methods are then used to solve the regularized problems.¹⁹^,²⁰

In our approach, we have used ultrasound-guided DOT to overcome the limitation of inaccurate target quantification by optical tomography alone.²¹ Having prior anatomical information from co-registered ultrasound (US), we are able to solve the inverse problem with the conjugate gradient (CG) iterative optimization method without using the regularization method. Our initial results obtained from more than 200 patients have demonstrated twofold higher absorption or total hemoglobin contrast between malignant and benign lesions.⁴ The feasibility of simultaneous reconstruction of absorption and scattering distributions using our approach has been investigated; however, the crosstalk between these two parameters often causes the CG to converge to the suboptimal values, especially in cases where there is a high scattering contrast between the target and the background.²² It is well known that gradient-based search algorithms often find local minima, which are closer to an initial guess or do not converge if the gradient is very small. In contrast, global optimization techniques such as genetic algorithm (GA) can solve problems with multiple solutions and are not prone to be trapped in the local minima; additionally, they are applicable to multidimensional, nondifferential problems. These algorithms seek to find the global optimum and are less sensitive to an initial guesses and regions with small gradients. It has been shown that evolutionary strategies can reliably find solutions to highly ill-posed inverse problems. The use of evolution strategy algorithms has been introduced for estimating the optical properties of homogenous media.²³^–²⁵ The feasibility of reconstructing bioluminance source distribution in the scattering medium using principles of evolutionary strategies is also shown.²⁶ An improved performance in absorption reconstruction with GA as compared with Tikhonov regularization for simulated time-gated data is also reported.²⁷

In this study, we investigated a new two-step approach that takes advantages of both global optimization and gradient techniques to solve the DOT inverse problem. The method utilizes the fitted values of the unknown target parameters estimated from GA as the initial values for the CG reconstruction. This new approach is validated using simulations and phantom experiments, and demonstrated with 16 clinical benign and malignant cases.

2. Method

2.1.

US-Guided Image Reconstruction

The propagation of diffused light through tissue can be described by photon diffusion approximation as Eq. (1):¹⁷

Eq. (1)

[\nabla^{2} + k^{2}] U (r) = - \frac{1}{D} S (r), k^{2} = \frac{- υ μ_{a} + j w}{D}, D = \frac{1}{3 μ_{s}^{'}},

where

S (r)

is the equivalent isotropic source and,

U (r)

is the photon density,

D

is the diffusion coefficient and

μ_{a}

,

μ_{s}^{'}

are the absorptions and reduced scattering coefficients, respectively. The inverse problem is typically linearized by Born or Rytov approximation. In Born approximation, the photon density is divided into superposition of the homogenous and scattered parts denoted by

U_{0}

and

U_{sc}

assuming that

U_{sc} ≪ U_{0}

. By digitizing the imaging space into

N

voxels, the resulting integral equations are formulated as follows:²⁸

Eq. (2)

{[U_{sc}]}_{M \times 1} = {[\begin{matrix} W^{A} & W^{S} \end{matrix}]}_{M \times 2 N} {[\begin{matrix} \int_{v} δ μ_{a} d v \\ \int_{v} δ D d v \end{matrix}]}_{2 N \times 1} = W X,

where

M

is the number of measurements,

δ μ_{a}

and

δ D

denote the unknown changes of absorption and diffusion coefficient at each voxel, respectively.

\int_{v} δ μ_{a} d v

and

\int_{v} δ D d v

are integral changes at each voxel. The weight matrices of absorption and scattering,

W^{A}

and

W^{S}

, describe the distribution of diffused wave in the homogenous medium and characterize the measurement sensitivity to the absorption and scattering changes. In the end, the inverse problem is formulated as an unconstrained optimization problem as:

\min {‖ Usc - W X ‖}^{2}

, where

‖ . ‖

is the Euclidean norm. In our ultrasound-guided DOT image reconstruction, a dual-zone mesh scheme is used to segment the imaging volume into a lesion region and a background region with fine and coarse voxel sizes, respectively.²¹ The reconstructed integral changes

\int_{v} δ μ_{a} d v

and

\int_{v} δ D d v

are divided by the different voxel sizes after CG optimization search, to obtain

δ μ_{a}

and

δ D

distributions. For computing the weight matrix, we used a semi-infinite absorbing boundary condition at the surface and the corresponding analytic solution of the diffusion equation. The

i j

’th element of

W^{A}

and

W^{S}

corresponding to

i

’th source-detector pair with the source and the detector located at

r_{si}

and

r_{di}

, the

j

’th voxel located at

r_{j}

, is calculated as:²⁸

Eq. (3)

\begin{array}{l} W_{i j}^{A} = \frac{1}{D_{0}} U_{0} (r_{j}, r_{si}) . G (r_{di} - r_{j}) \\ W_{i j}^{S} = \frac{- 1}{D_{0}} \nabla U_{0} (r_{j}, r_{si}) . \nabla G (r_{di} - r_{j}) \end{array} .

In Eq. (3),

D_{0}

is the background diffusion coefficient,

U_{0}

is the homogenous part of photon density, and

G

is the green function of the semi-infinite geometry. The average

μ_{a}

and

μ_{s}^{'}

of the background medium needed for the weight matrix computation are obtained from the measurements made with the homogenous intralipid solution for phantom experiments and the normal contra-lateral breast for clinical study. Finally after computing the weight matrix, we solve the aforementioned optimization problem in order to obtain the desired absorption and scattering distributions.

2.2.

Genetic Algorithm

The genetic algorithm (GA) is an approach used to solve optimization problems based on the survival of the fittest individuals.²⁹ First, a starting population is generated, formed by $N_{P}$ individuals while each member is a possible solution of the optimization problem and can be represented by a vector of real numbers. These individuals then undergo a set of genetic operations—selection, crossover, and mutation—in order to promote the population evolution.

Generally, a DOT problem has thousands of unknowns, so the problem cannot be solved directly using GA. Since optical properties of biological tissues usually change gradually rather than dramatically, the Gaussian perturbation is used to represent the solution as introduced in Ref. 30, where Bayesian inversion was applied to estimate the unknown parameters. We have added the scattering term to that model as formulated in Eq. (4), and the optimization is performed with genetic algorithm.

Eq. (4)

[\begin{matrix} μ_{a} \\ μ_{s}^{'} \end{matrix}] = [\begin{matrix} μ_{am} \\ μ_{sm}^{'} \end{matrix}] \exp [\frac{- {(P - C)}^{2}}{2 Σ^{2}}] .

This model has eight unknown parameters that are three-dimensional (3-D) center position

C : (c_{x}, c_{y}, c_{z})

, 3-D variance

2 Σ^{2} : (σ_{x}, σ_{y}, σ_{z})

), the maximum absorption and reduced scattering coefficient denoted as

μ_{am}

and

μ_{sm}^{'}

. The 3-D variance is calculated as:

2 Σ^{2} = - R^{2} \ln (0.25)

, where

R : (r_{x}, r_{y}, r_{z})

is the 3-D radius in (

x, y, z

) dimensions. This means that the radius is defined at the point where the target absorption/scattering coefficient is less than 25% of its maximum value. In our imaging procedure, first we examine the desired area with the handheld probe and monitor the ultrasound scans such that the target appears in the center of the imaging area. After that, the co-registered ultrasound and optical data are acquired. The target center and its radius in the

z

direction, denoted as

c_{z}

and

σ_{z}

, are also readily available using the ultrasound B-scan. Therefore GA algorithm is applied to obtain the rest of the parameters. Table 1 shows a pseudo code for the GA algorithm. First, an initial population of

N_{p} = 500

individuals is generated while each individual is a vector of four elements, i.e.,

[\begin{matrix} μ_{am} & {μ^{'}}_{sm} & r_{x} & r_{y}] \end{matrix}

. For initialization, these unknown variables are sampled from a uniformly distributed random variables as:

μ_{am} ({cm}^{- 1}) \sim Uniform (0, 0.4)

,

{μ^{'}}_{sm} ({cm}^{- 1}) \sim Uniform (1, 40)

, and

r_{x} (cm), r_{y} (cm) \sim Uniform (r_{\min}, 4)

.

Table 1

Pseudo code for GA algorithm.

Initialization:

P^{0} = {Y i}^{0}, i = 1, \dots, N_{p}

and

Y i = [\begin{matrix} μ_{am} & {μ^{'}}_{sm} & r_{x} & r_{y}] \end{matrix}

Evaluate fitness of

P^{0}

While not terminated:

1. Select parents:

X_{p 1}

,

X_{p 2}

2. Form child via crossover:

X_{c} = X_{p 1} + R_{c} (X_{p 1} - X_{p 2})

and

R_{c} \sim Uniform

(0,1).

3. Mutate the child:

{\tilde{X}}_{c} = X_{c} + R_{m}

and

R_{m} \sim Normal

(0,

σ_{m}

)

4. Evalute the fitness of the generated child

5. Choose the fittest member among (

X_{p 1}

,

X_{p 2}

,

{\tilde{X}}_{c}

) as a member of the new population

6. Repeat 1 to 5 until

N_{p}

members are added to the new population

End

The ranges of the distributions for absorption and scattering are selected according to the typical values of lesion optical properties at NIR wavelength. The minimum radius, $r_{\min}$ is chosen according to the target size measured by the US, and maximum radius is limited to the imaging field of view, which is half of our probe size. The error of each individual is defined as a square error between the measurement and the model and is formulated in Eq. (5) as,

Eq. (5)

Er = {‖ U_{sc} - W X ‖}^{2}, X = [\begin{matrix} \int_{v} δ μ_{a} d v \\ \int_{v} δ D d v \end{matrix}] = [\begin{matrix} \int_{v} (μ_{a} - μ_{a 0}) d v \\ \int_{v} (μ_{s}^{' - 1} - μ_{s 0}^{' - 1}) / 3 d v \end{matrix}],

where

μ_{a 0}

,

{μ^{'}}_{s 0}

are the background absorption and reduced scattering coefficients, respectively. The fitness of each individual is then calculated by inversing this error i.e.,

1 / Er

. In Table 1, steps 1 to 6 show how a new population is generated with GA algorithm. In step 1, the parents are selected using the roulette-wheel selection. In this method, the probability of selecting an individual as a parent is proportional to its fitness value. This means that the members of the population with higher fitness are selected more often and are less likely to be eliminated. In step 2, in order to form a child, heuristic crossover is applied followed by mutation in step 3. Mutation is an operator that changes a piece of genetic information of the new individual, and it is necessary to prevent the algorithm from converging toward a local minimum condition. In this work, Gaussian mutation is performed by adding a random number taken from a normal distribution with zero mean to the generated child. The standard deviation of this normal distribution decreases by increasing the number of generations. The fitness of the child is evaluated in step 4 and compared with the fitness of its parents in step 5. After that, the most fit member among the parents and generated child becomes a member of next population. To generate

N_{p}

members for a new population, steps 1 to 5 are repeated

N_{p}

times. With this procedure, a new population is generated such that the members have generally higher fitness levels in comparison with the older ones. The generation of the populations is then repeated until the algorithm reaches its termination criteria; that is, when the improvement in the fitness of members is less than

10^{- 6}

, or the generation number exceeds 100. By performing this GA optimization, the most fit values for the unknown parameters of the model are obtained. Finally, in order to obtain the optical distributions, the result of the GA is applied as an initial guess of CG reconstruction method. The iterative CG optimization for solving DOT inverse problem is described in detail in Ref. 31. Briefly, CG starts the iterations with an initial guess of the solution,

X_{0}

. At iteration

k

, the solution is updated as

X^{k + 1} = X^{k} + α^{k} d^{k}

, where the step size,

α^{k}

, is obtained by line minimization technique, and the direction

d^{k}

is obtained by using the gradient of the objective function, Er. The stopping criterion for our CG method is defined such that the algorithm stops when the cumulative change of the objective function of three consecutive iterations is less than 15%.

2.3.

Phantom and Clinical Experiments

Phantom and clinical experiments were performed using a frequency domain NIR system co-registered with ultrasound. For this study, we used the laser diode of wavelength 780 nm, modulated at 140 MHz. The source was sequentially coupled to nine positions on a handheld probe via optical fibers, and the reflected light was detected by 14 light-guides of 3 mm diameter that were coupled to 14 parallel photomultipliers. The system is described in detail in Ref. 32.

3. Results

3.1.

Simulated Targets

The DOT fitting with GA followed by CG reconstruction was applied to estimate the optical properties of simulated targets of diameters 1, 1.5, and 2.5 cm in size, and indicated as small, medium, and big targets, respectively. The optical absorption coefficients of $μ_{a} = [0.02, 0.08, 0.16] {cm}^{- 1}$ and the scattering coefficients of $μ_{s}^{'} = [4, 7, 16, 34] {cm}^{- 1}$ were used for targets. This gives 36 cases of targets with different sizes and optical properties. The detected fluence, back reflected from the top boundary was calculated for all cases, while the targets were assumed to be embedded inside a turbid medium with $μ_{a} = 0.02 {cm}^{- 1}$ and $μ_{s}^{'} = 7 {cm}^{- 1}$ . The targets were located at different depths such that the distances of the medium boundary to the top of the targets were 0.5, 1.0, 1.5, and 2 cm.

The fitted maximum values of $μ_{a}$ and $μ_{s}^{'}$ of all targets obtained by GA are presented in Fig. 1(a) and 1(b), respectively, and the estimated lateral diameter of the targets is shown in Fig. 1(c). These parameters were then applied to the CG-based reconstruction method. As an example, for the medium-size target with $μ_{a} = 0.08 {cm}^{- 1}$ , $μ_{s}^{'} = 16 {cm}^{- 1}$ , the GA algorithm converged after 51 iterations as shown in Fig. 2, and the CG method initialized by fitted values converged after five iterations. The GA fitted values for this target were $μ_{a} = 0.09 {cm}^{- 1}$ , $μ_{s}^{'} = 14.3 {cm}^{- 1}$ , and the reconstructed values using CG after GA were $μ_{a} = 0.071 {cm}^{- 1}$ , $μ_{s}^{'} = 15.4 {cm}^{- 1}$ . The CG method that was also applied for simultaneous reconstruction of absorption and scattering started from zero initial value. It converged after 11 iterations to a local minimum with $μ_{a} = 0.047 {cm}^{- 1}$ , $μ_{s}^{'} = 8 {cm}^{- 1}$ . The CG method solved for the absorption coefficient only and started from a zero initial value and converged after six iterations. The maximum reconstructed absorption was $0.132 {cm}^{- 1}$ . For this example, the final error, Er, calculated after CG iterations using GA fitting was about 4.5 times smaller than the one that started from a zero initial value.

Fig. 1

The maximum fitted (a) absorption coefficient; (b) reduced scattering coefficient; and (c) diameter of simulated targets using GA. The targets are located at depths of 0.5 to 2 cm. The dotted line shows the true values.

Fig. 2

An example of the reduction of the mean error, Er, of the population calculated at each generation of GA.

The spatial absorption and scattering distributions of all cases were reconstructed simultaneously with CG method after GA fitting, and the result was compared with the one that started from a zero initial value. In addition, the spatial absorption distributions that were reconstructed with the CG method started from a zero initial value and assumed that there was no scattering contrast between the target and the background. The maximum $μ_{a}$ and $μ_{s}^{'}$ of each target were extracted from its corresponding reconstructed maps, and the results are shown in Fig. 3. The light-gray bars in Fig. 3(a) show the maximum $μ_{a}$ reconstructed with the CG method after applying the GA fitting, while the gray bars are results of the CG obtained from simultaneous reconstruction of $μ_{a}$ and $μ_{s}^{'}$ , starting from a zero value. The black bars are the maximum $μ_{a}$ values obtained from the CG reconstructions for $μ_{a}$ only, and which started from a zero initial value. The standard deviations in Fig. 3(a) are due to different target $μ_{s}^{'}$ , sizes and depths. The error of the reconstructed absorption for high- and low-contrast targets ( $μ_{a} = 0.16 {cm}^{- 1}$ and $μ_{a} = 0.08 {cm}^{- 1}$ ) using CG for $μ_{a}$ only and CG for simultaneous $μ_{a}$ and $μ_{s}^{'}$ started from a zero initial value are on average 18.5% and 52%. While with this method, the ratio of high contrast target of $μ_{a} = 0.16$ to the low contrast one of $μ_{a} = 0.08$ is 1.25 and 1.33, respectively. This error reduced to 14% and the contrast increased to 1.7 by applying the CG method after GA fitting. In addition, the standard deviation for all cases has also reduced by about five times after using the GA initial values. The light-gray and gray bars in Fig. 3(b) show the maximum $μ_{s}^{'}$ reconstructed with the CG method after applying the GA fitting, and CG started from a zero initial value, respectively. It was determined from calculation that the maximum scattering was reconstructed with 18% error on average using CG after GA fitting. It can be seen that if CG started from a zero initial value is used for the simultaneous reconstruction of absorption and scattering, it cannot get close to the true value and is trapped in the local minimum close to the background values. The final Er calculated using CG after GA fitting is about 3.7 times smaller than the one using CG started from a zero initial value. The Er obtained using zero-initialized CG for reconstructing $μ_{a}$ only is also 1.6 times smaller than the one for the simultaneous reconstruction of $μ_{a}$ and $μ_{s}^{'}$ .

Fig. 3

The maximum reconstructed (a) absorption coefficient and (b) reduced scattering coefficient of 1-, 1.5-, and 2.5-cm diameter simulated targets located at depths of 0.5 to 2 cm. The dotted line shows the true values.

Figure 4(a) and 4(b) shows the bar plots of the reconstructed maximum $μ_{a}$ and $μ_{s}^{'}$ versus the size of the targets, respectively. The absorption of targets are reconstructed with 8.3%, 14.5%, and 24% errors while scattering are reconstructed with 11.4%, 19.5%, and 25% errors for the small, medium, and big targets, respectively. It can be seen that the optical values of the small targets are estimated more accurately in comparison with the larger ones. Furthermore, the increase in the estimation error caused by increasing the absorption and scattering values is also seen in this figure. The reconstructed values of the big target with highest absorption of $0.16 {cm}^{- 1}$ and scattering of $32 {cm}^{- 1}$ have a maximum error of about 40%.

Fig. 4

The maximum reconstructed (a) absorption coefficient and (b) reduced scattering coefficient versus size calculated for simulated targets located at depths of 0.5 to 2 cm. The dotted line shows the true values.

The reconstructed maximum $μ_{a}$ by CG only versus the target scattering is demonstrated in Fig. 5(a). The crosstalk between the absorption and scattering parameters is clearly seen in this figure (i.e., absorption values have increased by increasing the target scattering). Figure 5(b) indicates that the absorption values are not affected by changes in the scattering when simultaneous reconstruction of both parameters is performed with CG method after GA fitting. The standard deviations shown in Fig. 5 are due to the different sizes and depths of the targets. In addition to reducing the cross-talk between these two parameters, it is seen that even for the case when there is no scattering contrast between target and background (i.e., $μ_{s}^{'} = 7$ ), the absorption values are reconstructed more accurately after GA fitting since the initial guess is closer to the true value.

Fig. 5

Maximum reconstructed absorption versus reduced scattering of simulated targets of different sizes located at depths of 0.5 to 2 cm: (a) reconstructed with CG method and (b) with CG after fitting with GA.

3.2.

Phantom Experiments

Phantom experiments were performed to validate the simulations. Phantom targets of diameter 1.5 and 2.5 cm were submerged in 0.8%-intralipid solution with the calibrated absorption coefficient $μ_{a} = 0.03 {cm}^{- 1}$ and reduced scattering coefficient $μ_{s}^{'} = 7.5 {cm}^{- 1}$ , to model the breast tissue optical properties. The targets used were thin glass balls filled with various intralipid concentrations of 0.4%, 0.16%, and 0.32%, which provided mean calibrated scattering coefficients of 3, 13, and, $24 {cm}^{- 1}$ , respectively. In another set of experiments, Indian ink was added to the different concentrations of the intralipid solutions to obtain a calibrated average absorption value of $0.085 {cm}^{- 1}$ . Measurements were made with our frequency domain NIR system co-registered with ultrasound while the targets were located at depths of 0.5 to 2 cm measured from the surface of the target to the imaging probe.

The absorption and scattering distributions were reconstructed with the new two-step method. The fitted lateral diameters of 1.5 and 2.5 cm of the phantom targets are 2.43 and 3.24 cm on average, with standard deviations of 0.22 and 0.14, respectively, as shown in Fig. 6. For the phantom target of diameter 1.5 cm located at depth of 1.5 cm, with calibrated $μ_{a} = 0.085 {cm}^{- 1}$ , $μ_{s}^{'} = 13 {cm}^{- 1}$ , the fitted values using GA were $μ_{a} = 0.088 {cm}^{- 1}$ , $μ_{s}^{'} = 14.48 {cm}^{- 1}$ . The maximum reconstructed coefficients using CG after GA were $μ_{a} = 0.091 {cm}^{- 1}$ , $μ_{s}^{'} = 15.9 {cm}^{- 1}$ , while the maximum absorption reconstructed with CG started from zero value was $0.1268 {cm}^{- 1}$ .

Fig. 6

Fitted diameter of phantom targets using GA. The targets are located at depths of 0.5 to 2 cm.

The maximum reconstructed absorption of all phantom targets is shown in Fig. 7(a). For the target phantoms with $μ_{a} = 0.085 {cm}^{- 1}$ , the error of maximum reconstructed absorption is reduced from 14% to 4% after using the GA fitted values. The error is reduced about 22% on average for all experimental cases. The standard deviation shown is due to the differences in depth and $μ_{s}^{'}$ of the targets and is reduced about 1.5 times by using the GA fitting. The bar plots in Fig. 7(b) show the reconstructed $μ_{s}^{'}$ for all phantom targets. The scattering coefficient of the medium-size targets is reconstructed with about 16% error, while the standard deviation is about 1.8. The scattering of the big-size target is reconstructed with less than 31% error for lower scattering coefficients of 3 and $13 {cm}^{- 1}$ while there is about 52% error for higher scattering coefficient of $24 {cm}^{- 1}$ . This result agrees with the simulated data that the algorithm is less accurate in reconstructing higher scattering coefficients of big targets.

Fig. 7

Maximum reconstructed (a) absorption and (b) reduced scattering of medium (1.5 cm diameter) and big (2.5 cm diameter) phantom targets located at depths of 0.5 to 2 cm.

3.3.

Clinical Studies

The two-step method is demonstrated with clinical cases obtained with our ultrasound-guided NIR system. The study protocol was approved by the local Institution Review Board (IRB), and all patients signed the informed consent. Patients were scanned in a supine position while multiple sets of optical reflectance measurements were made with co-registered ultrasound images at the lesion location and the normal contra-lateral breast of the same quadrant as the lesion.

Figure 8(a) shows the ultrasound B-scan of a small suspicious mass located at the 11 o’clock position of the left breast of patient #2. The location of the lesion is marked with the white dashed line. The estimated background tissue optical properties were $μ_{a} = 0.017 {cm}^{- 1}$ and $μ_{s}^{'} = 3.59 {cm}^{- 1}$ . The core needle biopsy revealed that the lesion was benign. The reconstructed absorption distribution with CG starting from zero initial values and CG using the GA fitted values are shown in Fig. 8(b) and 8(c), with the maximum values of 0.136 and $0.142 {cm}^{- 1}$ , respectively. Each subfigure is the spatial $x - y$ image of 8 by 8 cm, reconstructed at the depth marked on its title and with 0.5 cm increment in depth. Figure 8(d) shows the reconstructed reduced scattering map of this lesion with maximum value of $4.81 {cm}^{- 1}$ . It can be seen that there is a small contrast between the $μ_{s}^{'}$ of the lesion and the background region. Similarly, Fig. 9 shows a co-registered ultrasound and the optical maps of a suspicious lesion located in the left breast of patient #11. The estimated optical properties of background tissue were $μ_{a} = 0.048 {cm}^{- 1}$ and $μ_{s}^{'} = 4.09 {cm}^{- 1}$ . The maximum absorption coefficients were 0.182 and $0.19 {cm}^{- 1}$ reconstructed with CG method only, and after applying the GA fitting, respectively. The spatial reduced scattering map is also shown in Fig. 9(d), with the maximum value of $10.38 {cm}^{- 1}$ , which has a high contrast as compared with the background tissue. The core biopsy revealed that the lesion was malignant. The optical maps of benign and malignant lesions reconstructed with and without applying the GA fitting were evaluated for seven malignant and nine benign cases. Table 2 provides patient information, lesion size and center depth as measured by co-registered ultrasound. All lesion diameters are smaller than 2 cm in the $z$ -direction and located at center depths of 1 to 3 cm. The lateral diameters of the lesions estimated with the GA algorithm are also presented in the last column of Table 2. They are on average 2.6 times larger than the size measured with the ultrasound image. The measured background tissue scattering varies from 1.7 to $6.7 {cm}^{- 1}$ and the reconstructed lesion scattering changes from 1.9 to $23 {cm}^{- 1}$ . Therefore a scattering-contrast term has been introduced and is defined as the ratio of absolute difference of lesion and background scattering to the background scattering. It was also reported in Ref. 6 that only the relative scattering of lesion to background showed a significant difference between the benign and the malignant tumor cases. The maximum reconstructed absorptions are shown in Fig. 10, and the scattering-contrasts of all cases are presented in Fig. 11. These figures imply that the maximum absorption and scattering contrast are significant parameters for characterizing the lesions. Figure 12(a) shows the statistics of maximum reconstructed absorption of all lesions. The ratio of the maximum absorption coefficient of malignant to benign cases improved slightly from 1.75 to 1.82 by applying the GA fitting before CG reconstruction. The scattering contrast shown in Fig. 12(b) of the malignant cases is also found to be about 3.32 times higher than that of the benign cases on average.

Fig. 8

A benign lesion obtained from patient #2: (a) ultrasound B-scan, spatial absorption map reconstructed at depths from 0.5 to 3 cm with (b) CG method and (c) after GA fitting; (d) reconstructed scattering distribution after GA fitting.

Fig. 9

A malignant lesion obtained from patient #11: (a) ultrasound B-scan, spatial absorption map reconstructed at depths from 0.5 to 3 cm with (b) CG method and (c) after GA fitting; (d) reconstructed scattering distribution after GA fitting.

Fig. 10

Maximum absorption of nine benign (1 to 9) and seven malignant (10 to 16) cases reconstructed with CG method started from zero initial values compared with the CG after GA fitting.

Fig. 11

Scattering contrast calculated for nine benign (one to nine) and seven malignant (10 to 16) cases.

Fig. 12

(a) Maximum reconstructed absorption and (b) scattering contrast of the nine benign and seven malignant cases.

Table 2

Patient information and lesion size in lateral and depth dimensions and center depth measured by co-registered ultrasound. The lateral dimension estimated by GA algorithm is also given.

Patient age	Lesion type	US measurements ( $lateral \times depth$ in cm)	GA estimated lateral size (cm)	Target depth measured by US (cm)
82	Benign papillary	$1.16 \times 0.6$	5.02	1.4
33	Benign breast tissue	$0.53 \times 0.75$	2.5	1
49	Benign fibrocystic changes	$0.87 \times 0.8$	2.42	1.4
46	Benign fibrosis	$0.53 \times 1.6$	2.24	1.25
34	Fibroadenoma	$2.38 \times 1.2$	3.39	1.5
26	Fibroadenoma	$2.2 \times 1.5$	6.03	1.25
39	Fibrocystic changes and moderate degree of intraductal hyperplasia	$5.3 \times 2.4$	7.84	2
76	Complex cyst	$2.1 \times 1.2$	5.56	2
31	Fibroadenoma	$5.2 \times 3$	7.9	2
61	Invasive ductal Ca.	$1.1 \times 1.0$	2.77	1.7
51	Invasive lobular Ca.	$1.1 \times 0.7$	2.76	1.5
61	Invasive Ca. with ductal and lobular features	$0.9 \times 0.6$	2.88	1.25
60	Invasive ductal Ca.	$1.7 \times 1.1$	3.68	1.8
53	Invasive ductal Ca.	$2.2 \times 1.9$	4.17	1.7
78	Invasive ductal Ca.	$1.8 \times 2.0$	3.86	3
76	Invasive lobular Ca.	More than 3 cma	5.2	2
39	Inflammatory ductal Ca.	$3.1 \times 1.4$	7.1	2.0

^aRecurrent cancer in the previous surgical site. Tumor size cannot be estimated from US due to scar tissue. The size was estimated from x-ray CT.

4. Summary and Discussion

We have introduced a new two-step imaging reconstruction algorithm based on a global optimization approach for estimating target parameters and then reconstructing the absorption and scattering distributions with the conjugate-gradient method. The GA algorithm can converge to the optimum solution with less than 60 iterations and takes less than 45 s of processing using a personal computer with CPU of 2.40 GHz. Results of simulation and phantom experiments with targets located at a depth range of 0.5 to 2 cm have shown that the new algorithm recovers the target absorption with errors less than 10%, while the use of the traditional conjugate-gradient approach results in about 20% error. Therefore 10% more accuracy has been achieved with the use of global optimization. With this approach, it is possible to reconstruct the reduced scattering distribution of targets with error less than 25% for the simulated and phantom targets. We have applied the new method for reconstructing the optical maps of nine benign and seven malignant clinical cases. The result shows a slight improvement in the absorption contrast of malignant to benign lesions. However, the scattering contrast showed about 3.32 times higher contrast on average in malignant lesions compared with benign ones.

There are some malignant cases that may not have high absorption contrast but high scattering contrast. For example, the patient #17 had an advanced inflammatory breast cancer with no measurable vascular content. The maximum reconstructed absorption with CG method only was $0.028 {cm}^{- 1}$ , and it increased to $0.057 {cm}^{- 1}$ after GA fitting. Using either of these values, this lesion would be characterized as benign. However the scattering-contrast for this case, calculated as 2.44, is in the malignant lesion range. This example demonstrates that tumor scattering distribution can add significant diagnostic value to further improve the lesion characterization. This result is applicable for the accurate diagnosis of breast cancers once it is validated by more clinical cases.

Although the algorithm here is applied for estimating the optical properties of one target, it is straightforward to expand the method for general cases with multiple targets as follows: We assume that the optical map of each target follows a Gaussian distribution as in Eq. (4). Having the prior information of center location and radius in $z$ -direction for each target, we need to estimate four unknowns for each target in this model. So we generate the population with members of size four times of the number of targets and perform GA fitting to find the optimum values for these parameters. Finally we apply CG iteration starting from the fitted parameters to reconstruct the optical distribution of all targets.

Another point worth to note is that there could be a mismatch between the centers of the lesion optical map and ultrasound image because the contrast mechanism of the two modalities is different. Therefore the lateral center positions were added as unknown variables and the fitted values were estimated with the GA algorithm. It was seen that the final optical reconstruction results were not sensitive to these parameters especially for the large lesions.

We have used genetic algorithm as the first step since the promising result of applying evolutionary algorithms for reconstructing tissue optical properties was studied before. Other global optimization techniques can be investigated in future for solving DOT inverse problem.

In summary, we have introduced a two-step image reconstruction algorithm by estimating initial target optical properties using a global optimization procedure and then applying the initial values to the gradient search algorithm to solve DOT inverse problem. Simulations and phantom experiments have shown that the maximum absorption and reduced scattering coefficients are reconstructed with less than 10% and 25% errors, respectively. However, the CG method alone generates about 20% error for the absorption coefficient and does not accurately recover the scattering distribution. The scattering contrast introduced to characterize benign and malignant breast lesions demonstrated that, on average, the reconstructed target absorption and scattering contrast of malignant lesions were about 1.8 and 3.32 times higher than the benign cases, respectively.

Acknowledgments

The authors thank the funding support of this work from the National Institute of Health (R01EB002136) and the Donaghue Medical Research Foundation.

References

1.

B. Chanceet al., “Breast cancer detection based on incremental biochemical and physiological properties of breast cancers: a six-year, two-site study,” Acad. Radiol., 12 (8), 925 –933 (2005). http://dx.doi.org/10.1016/j.acra.2005.04.016 1076-6332 Google Scholar

2.

R. Choeet al., “Differentiation of benign and malignant breast tumors by in-vivo three-dimensional parallel-plate diffuse optical tomography,” J. Biomed. Opt., 14 (2), 024020 (2009). http://dx.doi.org/10.1117/1.3103325 JBOPFO 1083-3668 Google Scholar

3.

Q. Fanget al., “Combined optical and x-ray tomosynthesis breast imaging,” Radiology, 258 (1), 89 –97 (2011). http://dx.doi.org/10.1148/radiol.10082176 RADLAX 0033-8419 Google Scholar

4.

Q. Zhuet al., “Early-stage invasive breast cancers: potential role of optical tomography with US localization in assisting diagnosis,” Radiology, 256 (2), 367 –378 (2010). RADLAX 0033-8419 Google Scholar

5.

C. Liet al., “Noninvasive in vivo tomographic optical imaging of cellular morphology in the breast: possible convergence of microscopic pathology and macroscopic radiology,” Med. Phys., 35 (6), 2493 –2501 (2008). http://dx.doi.org/10.1118/1.2921129 MPHYA6 0094-2405 Google Scholar

6.

B. Pogueet al., “Near-infrared scattering spectrum differences between benign and malignant breast tumors measured in vivo with diffuse tomography,” in Biomedical Topical Meeting, (2004). Google Scholar

7.

H. Jianget al., “Near-infrared optical imaging of the breast with model-based reconstruction,” Acad. Radiol., 9 (2), 186 –194 (2002). http://dx.doi.org/10.1016/S1076-6332(03)80169-1 1076-6332 Google Scholar

8.

G. M. Palmeret al., “Monte Carlo-based inverse model for calculating tissue optical properties. Part II: application to breast cancer diagnosis,” Appl. Opt., 45 (5), 1072 –1078 (2006). http://dx.doi.org/10.1364/AO.45.001072 APOPAI 0003-6935 Google Scholar

9.

A. Garofalakiset al., “Optical characterization of thin female breast biopsies based on the reduced scattering coefficient,” Phys. Med. Biol., 50 (11), 2583 –2596 (2005). http://dx.doi.org/10.1088/0031-9155/50/11/010 PHMBA7 0031-9155 Google Scholar

10.

C. Zhuet al., “Diagnosis of breast cancer using diffuse reflectance spectroscopy: comparison of a Monte Carlo versus partial least squares analysis based feature extraction technique,” Lasers Surg. Med., 38 (7), 714 –724 (2006). http://dx.doi.org/10.1002/(ISSN)1096-9101 LSMEDI 0196-8092 Google Scholar

11.

L. Spinelliet al., “Characterization of female breast lesions from multiwavelength time-resolved optical mammography,” Phys. Med. Biol., 50 (11), 2489 –2502 (2005). http://dx.doi.org/10.1088/0031-9155/50/11/004 PHMBA7 0031-9155 Google Scholar

12.

X. Guet al., “Differentiation of cysts from solid tumors in the breast with diffuse optical tomography,” Acad. Radiol., 11 (1), 53 –60 (2004). http://dx.doi.org/10.1016/S1076-6332(03)00562-2 1076-6332 Google Scholar

13.

B. Brooksbyet al., “Imaging breast adipose and fibroglandular tissue molecular signatures by using MRI-guided near-infrared spectral tomography,” Proc. Natl. Acad. Sci. U. S. A., 103 (23), 8828 –8833 (2006). http://dx.doi.org/10.1073/pnas.0509636103 PNASA6 0027-8424 Google Scholar

14.

A. Cerussiet al., “In vivo absorption, scattering, and physiologic properties of 58 malignant breast tumors determined by broadband diffuse optical spectroscopy,” J. Biomed. Opt., 11 (4), 044005 (2006). http://dx.doi.org/10.1117/1.2337546 JBOPFO 1083-3668 Google Scholar

15.

P. Taroniet al., “Noninvasive assessment of breast cancer risk using time-resolved diffuse optical spectroscopy,” J. Biomed. Opt., 15 (6), 060501 (2010). http://dx.doi.org/10.1117/1.3506043 JBOPFO 1083-3668 Google Scholar

16.

D. A. Boaset al., “Imaging the body with diffuse optical tomography,” IEEE Signal Process. Mag., 18 (6), 57 –75 (2001). http://dx.doi.org/10.1109/79.962278 ISPRE6 1053-5888 Google Scholar

17.

P. K. Yalavarthyet al., “Weight matrix structured regularization provides optimal generalized least squares estimate in diffuse optical tomography,” Med. Phys., 34 (6), 2085 –2098 (2007). http://dx.doi.org/10.1118/1.2733803 MPHYA6 0094-2405 Google Scholar

18.

B. W. Pogueet al., “Spatially variant regularization improves diffuse optical tomography,” Appl. Opt., 38 (13), 2950 –2961 (1999). http://dx.doi.org/10.1364/AO.38.002950 APOPAI 0003-6935 Google Scholar

19.

R. RoyE. Sevick-Muraca, “A numerical study of gradient-based nonlinear optimization methods for contrast enhanced optical tomography,” Opt. Express, 9 (1), 49 –65 (2001). http://dx.doi.org/10.1364/OE.9.000049 OPEXFF 1094-4087 Google Scholar

20.

R. RoyE. Sevick-Muraca, “Truncated Newton’s optimization scheme for absorption and fluorescence optical tomography: Part I theory and formulation,” Opt. Express, 4 (10), 353 –371 (1999). http://dx.doi.org/10.1364/OE.4.000353 OPEXFF 1094-4087 Google Scholar

21.

N. G. Chenet al., “Simultaneous near-infrared diffusive light and ultrasound imaging,” Appl. Opt., 40 (34), 6367 –6380 (2001). http://dx.doi.org/10.1364/AO.40.006367 APOPAI 0003-6935 Google Scholar

22.

M. M. Huanget al., “Simultaneous reconstruction of absorption and scattering maps with ultrasound localization: feasibility study using transmission geometry,” Appl. Opt., 42 (19), 4102 –4114 (2003). http://dx.doi.org/10.1364/AO.42.004102 APOPAI 0003-6935 Google Scholar

23.

A. J. HielscherA. D. KloseJ. Beuthan, “Evolution strategies for optical tomographic characterization of homogeneous media,” Opt. Express, 7 (13), 507 –518 (2000). http://dx.doi.org/10.1364/OE.7.000507 OPEXFF 1094-4087 Google Scholar

24.

A. D. KloseA. H. Hielscher, “An approach for diffuse optical tomography combining evolution strategies and gradient techniques,” Proc. SPIE, 4250 11 –19 (2001). http://dx.doi.org/10.1117/12.434484 PSISDG 0277-786X Google Scholar

25.

M. L. Flexmanet al., “A wireless handheld probe with spectrally constrained evolution strategies for diffuse optical imaging of tissue,” Re. Sci. Instrum., 83 (3), 033108 (2012). http://dx.doi.org/10.1063/1.3694494 RSINAK 0034-6748 Google Scholar

26.

A. D. Klose, “Transport-theory-based stochastic image reconstruction of bioluminescent sources,” JOSA, 24 (6), 1601 –1608 (2000). http://dx.doi.org/10.1364/JOSAA.24.001601 JSDKD3 Google Scholar

27.

Q. Zhaoet al., “Reconstruction in diffuse optical tomography using genetic algorithm,” in Biomedical Optics, (2010). Google Scholar

28.

Q. ZhuN. ChenS. J. Kurtzman, “Imaging tumor angiogenesis by use of combined near-infrared diffusive light and ultrasound,” Opt Lett., 28 (5), 337 –339 (2003). http://dx.doi.org/10.1364/OL.28.000337 OPLEDP 0146-9592 Google Scholar

29.

M. Mitchell, An Introduction to Genetic Algorithms, MIT Press, Cambridge, MA (1996). Google Scholar

30.

J. Liuet al., “Parametric diffuse optical imaging in reflectance geometry,” IEEE J. Selected Topics in Quan. Electron., 16 (3), 555 –564 (2010). IJSQEN 1077-260X Google Scholar

31.

W. Zhuet al., “Iterative total least-squares image reconstruction algorithm for optical tomography by the conjugate gradient method,” J. Opt. Soc. Am. Optic. Image Sci. Vis., 14 (4), 799 –807 (1997). JSDKD3 Google Scholar

32.

Q. Zhuet al., “Optimal probing of optical contrast of breast lesions of different size located at different depths by U.S. localization,” Technol. Cancer Res. Treat., 5 (4), 365 –380 (2006). 1533-0346 Google Scholar

Citation Download Citation

Behnoosh Tavakoli and Quing Zhu "Two-step reconstruction method using global optimization and conjugate gradient for ultrasound-guided diffuse optical tomography," Journal of Biomedical Optics 18(1), 016006 (7 January 2013). https://doi.org/10.1117/1.JBO.18.1.016006

Published: 7 January 2013

Access the abstract

JOURNAL ARTICLE
11 PAGES

DOWNLOAD PAPER SAVE TO MY LIBRARY

GET CITATION

CITATIONS

Cited by 19 scholarly publications.

Explore citations on Lens.org

RIGHTS & PERMISSIONS

Get copyright permission Get copyright permission on Copyright Marketplace

KEYWORDS

Scattering

Absorption

Tissue optics

Breast

Reconstruction algorithms

Ultrasonography

Optical properties

1.

Introduction

2.

Method

2.1.

US-Guided Image Reconstruction

Eq. (1)

Eq. (2)

Eq. (3)

2.2.

Genetic Algorithm

Eq. (4)

Table 1

Eq. (5)

2.3.

Phantom and Clinical Experiments

3.

Results

3.1.

Simulated Targets

Fig. 1

Fig. 2

Fig. 3

Fig. 4

Fig. 5

3.2.

Phantom Experiments

Fig. 6

Fig. 7

3.3.

Clinical Studies

Fig. 8

Fig. 9

Fig. 10

Fig. 11

Fig. 12

Table 2

4.

Summary and Discussion

Acknowledgments

References

Show All Keywords

Keywords/Phrases

Search In:

Publication Years