Two step imaging reconstruction using truncated pseudoinverse as a preliminary estimate in ultrasound guided diffuse optical tomography

: Due to the correlated nature of diffused light, the problem of reconstructing optical properties using diffuse optical tomography (DOT) is ill-posed. US-, MRI- or x-ray-guided DOT approaches can reduce the total number of parameters to be estimated and improve optical reconstruction accuracy. However, when the target volume is large, the number of parameters to estimate can exceed the number of measurements, resulting in an underdetermined imaging model. In such cases, accurate image reconstruction is difficult and regularization methods should be employed to obtain a useful solution. In this manuscript, a simple two-step reconstruction method that can produce useful image estimates in DOT is proposed and investigated. In the first step, a truncated Moore-Penrose Pseudoinverse solution is computed to obtain a preliminary estimate of the image that can be reliably determined from the measured data; subsequently, this preliminary estimate is incorporated into the design of a penalized least squares estimator that is employed to compute the final image estimate. By use of phantom data, the proposed method was demonstrated to yield more accurate images than those produced by conventional reconstruction methods. The method was also evaluated with clinical data that included 10 benign and 10 malignant cases. The capability of reconstructing high contrast malignant lesions was demonstrated to be improved by use of the proposed method.

A variety of image reconstruction methods have been employed to improve target reconstruction accuracy in DOT. This includes the algebraic reconstruction technique (ART) [14], nonlinear iterative gradient based optimization methods [15][16][17], and Newton-like methods that requires the direct calculation and inversion of the Jacobian or weight matrix [18,19]. Reconstruction methods using a prior information determined from co-registered high resolution MRI [9, 20-24], x-ray imaging [10,11], and ultrasound imaging (US) [5,6,12,21] have been investigated extensively. These methods segment the lesion and background regions or different tissue types seen by a high resolution modality, and therefore reduce voxels with unknown optical properties and improve the ill-posed and underdetermined DOT reconstruction problem.
Another means of incorporating a prior information into an iterative image reconstruction method is through the initial image estimate. Behnoosh and Zhu proposed a two step reconstruction using Genetic Algorithm (GA) to find a suitable initial image estimate that was subsequently refined by use of a conjugate gradient (CG) method, which showed improved target quantification as compared to CG with zero initial estimate [25]. However, GAs are time-consuming and an optimal CG stopping criterion for use with experimental data is difficult to specify.
In this manuscript, a simple two-step reconstruction method that can produce useful image estimates in DOT is proposed and investigated. In the first step, a truncated Moore-Penrose Pseudoinverse (MPP) solution is computed to obtain a preliminary estimate of the target image that can be reliably determined from the measured data; subsequently, this preliminary estimate is incorporated into the design of a penalized least squares estimator that is employed to compute the final image estimate. The MPP was employed to compute the initial estimate of the target image for several reasons. Firstly, the MPP pseudoinverse, by definition, produces the least-squares estimate of the image that possess the minimum norm. This yields a solution that can be interpreted as an orthogonal projection of the true target image onto a subspace that is the orthogonal complement to the null space of the imaging operator. Therefore, the MPP solution describes the estimate of the target that is closest to the true target but contains no component in the null space. This is a reliable strategy for image reconstruction when no reliable a priori information about the target is available. Secondly, the MPP pseudoinverse solution can be easily regularized by excluding contributions that correspond to small values. Therefore, there is little ambiguity in how to choose the regularization parameter. Thirdly, the MPP pseudoinverse operator for our system can be explicitly stored in memory, which leads to near real-time image reconstruction. By use of phantom data, the proposed method was demonstrated to yield more accurate images than those produced by conventional reconstruction methods. The method was also evaluated with clinical data that included 10 benign and 10 malignant cases. The capability of reconstructing high contrast malignant lesions was demonstrated to be improved by use of the proposed method.

Constrained optimization
The propagation of diffused light through tissue can be described by photon diffusion approximation as [26]: where S(r) is the equivalent isotropic source and, U(r) is the photon density wave, D is the diffusion coefficient, a μ and ' s μ are the absorption and reduced scattering coefficients, respectively. The inverse problem is typically linearized by Born approximation [26]. By digitizing the imaging space into N voxels, the resulted integral equations are formulated as following: where sc U is the measured scattered photon density wave, M is the number of measurements, denotes the unknown changes of absorption coefficient at each voxel. The weight matrix, W, describes the distribution of diffused wave in the homogenous medium and characterizes the measurement sensitivity to the absorption and scattering changes. At the end, the inverse problem can be formulated as an unregularized optimization problem as: In our ultrasound-guided DOT image reconstruction, a dual-zone mesh scheme is used to segment the imaging volume into a lesion region identified by co-registered US and a background region with fine and coarse voxel sizes, respectively [27]. This scheme effectively reduces the total number of voxels with unknown optical properties. The conjugate gradient (CG) method is utilized to iteratively solve the inverse problem. As a result, the target quantification accuracy can be significantly improved. However, when the lesions are larger, the total number of finer voxels and coarse voxels, N, can be much larger than the total measurements, M, which is the number of sources × the number of detectors × 2 = 14 × 9 × 2 = 252 counting for both amplitude and phase data. Due to the correlated nature of diffused light measured at closely spaced source and detector positions and also measurement noise, increasing the number of sources and detectors does not effectively mitigate the ill-conditioned nature of the DOT inversion problem.
In this manuscript, we formulate the inverse problem as: where 0 X is a preliminary estimate of the optical properties that can be reliably determined from the measured data and λ is a regularization parameter. A Newton-like or conjugate gradient optimization method will be employed to approximately solve Eq. (4). No spatial or temporal filters were used on solution of ( ).

Truncated pseudoinverse as an initial estimate
We propose to employ a truncated pseudoinverse (PINV) operator MPP of W is, From system of linear equations, Eq. (2), Since our measurement contains noise, we assume additive noise n, .
sc noiseless U U n = + Then the reconstructed absorption X  is given as: For very small singular values 0 n σ → , noise X may contain image artifacts.
In the truncated MPP approach, a threshold th σ is set for singular values and the initial solution using MPP is: In the phantom and clinical data, we have chosen th σ as 10% of 1 σ as a cut off value.
From the truncated pseudoinverse, a preliminary estimate of unknown optical properties can be obtained. A simple projection operation is used to suppress pixels outside the region of interest identified by a sphere B obtained from measurements of co-registered ultrasound image. This projected absorption map is used as an initial solution for Newton or Conjugate gradient search method.

Newton method
The Newton method uses 2 nd derivative of objective function (known as hessian) to calculate a 2 nd order search direction resulting in quadratic convergence rate. We reformulate our penalized least square problem as a quadratic optimization problem, Clearly, the hessian is positive definite when 0 λ > . Our solution is iteratively updated using following equations, The iteration process is terminated when change of objective function between successive iterations become smaller than a preset tolerance. Choice of the regularization parameter is a critical part of solution design. Based on tumor size measured from ultrasound image and largest singular value of weight matrix, 1 σ we empirically chose our regularization parameter where p is proportional to tumor size.

Conjugate gradient method
The Conjugate gradient (CG) method is a well-known iterative technique for solving symmetric positive definite linear systems of equations. We investigated this method both with regularization and without regularization. For the unregularized optimization formulation as given in Eq. (3)

Comparison of five reconstruction methods
Five reconstruction methods have been compared using phantom and clinical data. Using zero as an initial estimate of target optical properties and regularized Newton optimization (Newton Zero initial) and regularized CG optimization (CG Zero initial); using PINV as an initial estimate of target optical properties and regularized Newton optimization (Newton PINV initial), regularized CG optimization (CG PINV initial), and using zero initial estimate and unregularized CG. Additionally, target centroid error i.e. the absolute difference between the center of a phantom target measured by co-registered US and the centroid of corresponding reconstructed target absorption map, is calculated as a measure of reconstruction quality.

US-guided DOT system, data acquisition
Our US-guided DOT system has been reported earlier [

Phantom experiments
Phantom experiments were performed with solid ball phantoms of different sizes and different optical contrasts emulating tumors. These targets were submerged in 0.  Table 1 and shown in Fig. 1. Errors of both high and low contrast targets reconstructed using different methods are given in Table 2. As seen from the Table, Newton and CG with PINV as an initial accurately estimate absorption coefficient while Newton and CG with a zero initial produce larger errors. Unregularized CG gives a better estimate for high contrast phantoms but results in under reconstruction for low contrast ones.

Patient data
Performance of proposed method is demonstrated using clinical data obtained from 20 patients [12]. Based on biopsy results, 10 patients had benign lesions and 10 patients had cancers. The study was approved by the local Institution Review Boards (IRB) and was compliant with the Health Insurance Portability and Accountability Act (HIPPA). Informed consent was given by each patient. Data used in this study have been de-identified.
An example of a cancer case is shown in Fig. 2. Figure 2(a) is the co-registered US image with the suspicious lesion marked by a circle. Absorption maps of PINV initial image from truncated PINV (b) and reconstructed images using Newton with zero initial (c), PINV initial (d), CG with zero initial (e), CG with PINV initial (f) and unregularized CG (g) have shown similar lesion position and shape, however, the Newton's method with PINV initial yields highest reconstructed µ a = 0.268 cm −1 . An example of a benign lesion is shown in Fig. 3. Figure 3(a) is the co-registered US image with the suspicious lesion marked by a circle. Absorption map of MPP estimated image is shown in Fig. 3(b), reconstructed images using five corresponding reconstruction methods are given in Fig. 3(c)-3(g) and reconstructed maximum µ a are quiet similar. 11.6 ± 13.8% 0.04 ± 9.1% 11.8 ± 13.2% 0.1 ± 9.0% 3.5 ± 9.9%

Error (high contrast)
12.0 ± 16.1% 9.6 ± 14.6% 15.6 ± 10.9% 8.8 ± 15.8% 26.5 ± 8.5% Fig. 1. Box plot of phantom data obtained from 1 to 3 cm size absorbers of high contrast (red) and low contrast (blue) located at different depths (1.5-3.5 cm center depth) using zero and PINV as an initial guess and Newton as optimization, respectively (first and second columns), zero and PINV as initial guess and CG, respectively (third and fourth columns), and unregularized CG (last column).
Box plot of maximum total hemoglobin (tHb) of all clinical cases is shown in Fig. 4. The two sample t test was performed between malignant and benign groups of each method. Newton's method and CG with PINV as an initial provide highest statistical significance.
Additionally, the malignant to benign contrast ratios are 1.61, 2.11, 1.61, 2.07, 1.93, for Newton's with zero initial, PINV initial, CG zero initial, PINV initial and unregularized CG respectively. The average and standard deviation of maximum tHb concentration obtained from each method is given in Table 3. For benign cases, reconstructed tHb are comparable using five methods, however, for malignant cases the total Hb contrast is much higher when Newton's and CG are used with PINV initial estimate.  and each sub-image shows spatial x and y distribution of absorption coefficients reconstructed from 0.5 cm to 3.5 cm depth range from the skin surface. The spacing between the sub-images in depth is 0.5 cm. The color bar is absorption coefficients in cm −1 . We chose the µ a display range from 0 to 0.2 cm −1 because most of the reconstructed absorption values fall within this range. Each subfigure dimension is 8cm x 8cm with scales from −4 cm to 4 cm in both X and Y axis.

Convergence speed analysis
It is important for an iterative image reconstruction algorithm to converge quickly and therefore to provide images needed for on-site diagnosis by physicians. To compare the convergence of differenrt reconstruction methods, we normalize the Least Square Error . Newton and CG with PINV as an initial estimate converge in 1 and 2 iterations, respectively. Newton and CG with zero initial converges in 1 and 3 iterations, respectively, and the residual LSE of CG is slightly higher than that with PINV as an initial. Unregularized CG converges in 3 ietrations. Note that for our early studies using unregularized CG, the iteration was stopped at 3 ietrations because it provided optimal performance for phantom data.

Target centroid error analysis
To compare different reconstruction methods, the target centroid error i.e. the absolute difference between the center of a phantom target measured by co-registered US and the centroid of corresponding reconstructed target absorption map, is calculated as a measure of reconstruction quality. Phantom data of both low and high contrast targets of 1 cm diameter located at different depths and measured at 780 nm were used to estimate the centroid error and results are shown in Table 4. MATLAB function 'regionprop' is used to estimate the centroid of target absorption map and the difference between the estimated centroid and the measured target center from corresponding co-registered US is calculated. As seen from Table 4, the target centroid error which is less than one voxel size of 0.25 cm does not depend on reconstruction method. Thus, all reconstruction methods provide the same target centroid.

Discussion and summary
Choice of regularization parameter, λ, is an important part of reconstruction. If λ is too small, then the penalty may not have any effect on reconstruction, however, a larger λ heavily penalizes data fidelity term and solution may not converge near true minimum of unregularized objective function. In our approach, λ is chosen as μ , the λ regulates more to improve the conditioning of the Q matrix. Additionally, because the huge difference between the first and the rest of the eigenvalues, the / n λ σ increases with n and therefore λ regulates more for smaller eigenvalues and further improves the conditioning of Q matrix.
Choice of regularization parameter is always a difficult problem and mathematical techniques like L-curve and U-curve are not often useful. We have determined regularization parameter by trial and error using phantom data to ensure convergence, reconstruction accuracy and lower image artifacts. In further study, we will apply machine learning techniques to automatically select regularization parameter to minimize the reconstruction error. The ultimate clinical use of ultrasound-guided diffused light imaging is to maximize the separation of benign and malignant lesions.
In the past two decades, researchers in DOT community have tried to simultaneously reconstruct target absorption coefficient, a μ , and diffusion coefficient, D ( , see Eq. (1)). However, since the lesion diffusion coefficient is much smaller than the absorption coefficient, correctly reconstructing the scattering coefficient is a challenge. Also, simultaneously reconstructing the absorption and diffusion coefficients doubles the number of unknown optical parameters to estimate. Therefore, the reconstruction becomes more illposed and under-determined. However, with a better initial estimate and an appropriate choice of regularization parameter λ, it is possible to explore simultaneous reconstruction of both parameters. This has been demonstrated in reference 25 using GA as an initial estimate and unregulated CG to iteratively reconstruct target absorption and scattering maps. In this manuscript, our objectives were to 1) evaluate the performance of the proposed simple, robust, two-step reconstruction algorithm; and 2) compare this algorithm with a group of four algorithms including the unregulated CG algorithm that we have used in the past. Therefore, we did not attempt to simultaneously reconstruct both parameters but focused on absorption coefficient, a μ , which is the most important parameter to reveal tumor angiogenesis. Thus, our phantoms have similar reduced scattering coefficient as the background medium. In future study, we will evaluate the performance of the proposed novel algrothm in simultaneously recovering both target absorption and scattering maps.
In summary, a simple, robust, two-step reconstruction algorithm has been proposed and its performance has demonstrated using phantom and clinical data. Using a truncated pseudoinverse as a preliminary estimate of target optical properties and regularized Newton and CG optimization search methods to iteratively reconstruct target optical properties within region of interest identified by co-registered US gave best results. The truncated pseudoinverse as a preliminary estimate and regularized Newton optimization converges in one iteration. This two-step reconstruction technique is generally applicable to x-ray-guided and MRI guided DOT imaging reconstruction.

Funding
The authors thank the funding support of this work from the National Institute of Health (R01EB002136).

Disclosures
QZ is the inventor of the patents related to ultrasound-guided near-infrared tomography technologies and patents owned by the University of Connecticut. The other authors declare there are no conflicts of interest related to this article.