Optoacoustic imaging combines the rich contrast of optical imaging and the high resolution of ultrasonic imaging (1–3). The method enables visualization of absorption chromophores that have high spatial resolution in the micrometer-to-millimeter range depending on the imaging depth. Therefore, optoacoustic imaging can provide anatomical, physiological, and molecular information at resolutions that are significantly improved over optical diffusion imaging techniques (1, 2, 4). The technique operates on illuminating the sample imaged by nanosecond laser pulses. Absorption of the light energy generates broadband ultrasonic waves via thermoelastic expansion; these waves have frequencies ranging from hundreds of kilohertz to many tens of megahertz (3, 5, 6). The recorded optoacoustic signals are used to reconstruct an image using analytical or model-based algorithms (7–9).
Different image reconstruction schemes can be used for image formation. The spherical radon transform is widely used because of its simple implementation and high efficiency. Model-based approaches are capable of incorporating information regarding detection geometry, acoustic attenuation, and transducer properties in the reconstruction process (10, 11), resulting in more accurate reconstructions. However, model-based methods require large numbers of repeated sparse matrix-vector multiplications in an iterative manner, resulting in significant computational cost (12, 13). Accelerated model-based methods were developed to reduce the computation cost (13–15). For example, the angular discretization method was used to generate the model matrix, which effectively reduced the computational cost and saved memory (13, 14). Other approaches performed inversion on parallel processing platforms based on graphics processing units, which enabled real-time, model-based reconstruction (15, 16).
A particular challenge in optoacoustic tomography is the implementation of limited-view projections, that is, cases where 360° projections are not available. This could be the case, for example, in imaging large volumes (whole-animals or humans), whereby access is afforded only from 1 side of the tissue, analogous to ultrasound imaging. Limited-view implementations typically result in lower image fidelity and a larger number of artifacts compared with 360° view data sets (12). Nevertheless, sparsity-based algorithms showed better performance with limited-view data sets (12, 17, 18) compared with Tikhonov-based reconstructions, albeit at a higher computational cost. Moreover, sparse recovery-based methods may amplify noise in limited-data scenarios (19).
In this work, we propose improvements to sparsity-based inversion for optoacoustic tomography. The sparse method proposed is implemented for accelerating the reconstruction process implemented within the Split Augmented Lagrangian Shrinkage Algorithm (SALSA) framework, using least square QR (LSQR) inversion and a novel coherence factor (CF) weighting scheme for suppressing noise and artifacts. We show the performance on several phantoms and biological tissue ex vivo to quantify performance improvements compared with previously described SALSA implementations.
Acoustic Forward Problem.
Regularization and Inversion.
Efficient inversion of Equation 3 requires regularization. We selected conventional Tikhonov regularization with parameter (λ), assuming an initial pressure rise distribution that is smoothly varying. In this case, the objective function to be minimized is given as follows:22), as follows:
However, Equation 5 is computationally expensive because of the time-consuming matrix inversion. Alternatively, an LSQR approach can be used (23) as follows:23, 24), and β0 is defined as ∥b∥22. e1 is [1 …]T. Equation 6 can be inverted in a faster fashion compared with Equation 5 because it involves inverting the bidiagonal matrix, which is computationally efficient.
Proposed SALSA Acceleration with CF.
The proposed method is based on applying a sparsity constraint and accelerating the reconstruction with the help of bidiagonal matrices. The accelerated SALSA (ASALSA) is proposed herein as an improved version of SALSA minimization implemented using Krylov subspace optimization. In this case, the objective function to be minimized is as follows:
Sparsity optimization schemes are expected to offer better performance over conventional Tikhonov regularization for limited projection data (25, 26). Equation 4 assumes a smooth solution, hence resulting in a large number of unknowns and edge smoothening. Equation 7 assumes that the number of unknowns are sparse (by considering only nonzero entries) and are known to perform well in limited-data scenarios. Equation 7 is minimized using the SALSA scheme, which has demonstrated the fastest convergence among existing sparsity norm-based optimization schemes (27). In this scheme, we use a variable splitting approach, wherein a new variable is introduced in the optimization procedure. The above objective function is now split into 2 quadratic minimizations with the help of the temporary variable (v), which is given as follows:
The original SALSA algorithm involved the inversion of a large matrix during the optimization procedure (27). To accelerate inversion, we recast the SALSA algorithm, as indicated in Table 1, by using the LSQR solver. Faster computations are achieved by using LSQR iterative inversion schemes for enabling accelerated SALSA (termed as ASALSA) reconstruction using the L1-norm-based approach. It can be seen that Equation 4 applies a smoothness constraint, (∥x∥2); hence, noise will be smoothed out during reconstruction. In contrast, because equation (7) applies sparsity constraint, it may amplify the weak signal and noise (19). To suppress noise amplification and artifacts that arise because of limited-view data and sparsity constraint, we introduce herein an additional operation using a CF (28, 29), defined as the ratio between the energy of the coherent sum of optoacoustic signals and the total incoherent energy, which is explained as follows:
|Aim: Estimation of x by solving Equation 7|
|Input: A, b, λ, α, max_iter|
|Initialize ADMM parameter d = 0.|
|1. Calculate backprojection solution (x = ATb), CF = AT b2|
|for k = 1, 2, … max_iter|
|2. Optimize Equation 9: v = soft(x + d, (0.5 × λ)/α)|
|3. Optimize Equation 8: x1 = , where matrix C is obtained by Lanczos bidiagonalization of AT with k iterations and is the euclidean norm of (ATb + α* (v − d)).|
|xrecon = , where is the Euclidean norm of x1.|
|4. Update ADMM parameter: d = d + xrecon − v|
|5. CF = (x2/CF)|
|6. x = xrecon. × CF|
This CF weighting enables amplification of the higher optoacoustic signal and suppresses the noise. The ASALSA algorithm along with the CF is indicated in Table 1.
Optoacoustic Imaging System
A multispectral optoacoustic tomography small animal scanner (MSOT256-TF, iThera Medical GmbH, Munich, Germany) was used to experimentally examine the performance of the proposed reconstruction method. In brief, a custom-made, 256-element cylindrically focused array was used to record the optoacoustic signals. The ultrasound array covered an angle of approximately 270° with a radius of 40 mm, allowing simultaneous acquisitions of the signals generated with each laser pulse. The central frequency of the array elements was 5 MHz with a bandwidth of 90%. The sample was illuminated with a wavelength-tunable optical parametric oscillator laser with a repetition rate of 10 Hz. The detected optoacoustic signals were simultaneously digitized at 40 megasamples/s and were averaged 10 times to improve the signal-to-noise ratio (SNR) of the signal. Detailed information about the imaging setup can be found in the literature (30–32).
Phantom and Tissue Experiments
To test the performance of the proposed method, a printed paper (U.S. Air Force resolution target, standard inkjet printer with black ink) embedded in a 1.9-cm-diameter diffuse agar cylinder (6% by volume intralipid in the agar solution) was imaged. The absorbing features of the phantom are shown in Figure 1A. The phantom consisted of several groups of line elements of different sizes, which can be used for resolution and image quality characterization at different levels. To mimic limited-view scenarios, we assumed 2 down-sampled data sets, one using 128 positions over 270° coverage angle and one using 128 positions over 135° coverage angle.
A tissue-mimicking agar phantom was also prepared to examine the ASALSA-CF performance. The phantom included areas containing 0.016% India Ink to impart higher optical absorption than the background. The absorption coefficient of the ink-containing inclusions was μa = 1.6 cm−1. In addition, a hollow cylindrical cavity was introduced for illustrating the nature of artifacts produced because of acoustic mismatches. The phantom was used to examine the performance of the proposed algorithm in relation to artifacts arising because of the reflecting material.
Furthermore, a mouse kidney was imaged ex vivo to examine the performance of the proposed method with the biological tissue. The kidney sample was extracted postmortem (nonperfused) according to institutional regulations regarding animal handling protocols and subsequently embedded in a diffuse agar block (6% by volume intralipid in the agar solution) for ensuring uniform illumination of the sample.
The performance of reconstructions based on Equation 6 (L2-norm), Equations 8 and 9 (ASALSA), and Equation 11 (ASALSA-CF) was compared using the experimental measurements collected from phantoms and mouse kidneys. Before reconstruction, the optoacoustic signals were band-pass filtered, with cutoff frequencies between 0.2 and 7 MHz to remove low-frequency offsets and high-frequency noise. A uniform speed of sound of 1510 m/s was used for all the reconstructions (33). For all phantom measurements, images were reconstructed with a pixel size of 100 μm [200 × 100 pixels (2)], and in the case of tissue data, a pixel size of 100 μm [200 × 180 pixels (2)] was used. The regularization parameter for the L2-norm-based scheme was obtained using an L-curve approach, while in the case of the ASALSA algorithm, it was chosen heuristically. The parameters α and λ were set as 100 and 1500, respectively, for the ASALSA algorithm. Note that in the case of the ASALSA algorithm, we have multiple parameters (α and λ, which are sensitive to noise; therefore, they can be adjusted based on the image quality of the reconstructed image.
To further quantify the reconstruction improvement, the average intensity (S) of the resolution lines was calculated for the paper phantom, as well as the standard deviation of background (B) signals in the spaces between the lines. The ratio of the average signal to the standard deviation of background noise provided a contrast ratio of CNR = 20 ∗ log(S/B).
The reconstruction results corresponding to the printed-paper U.S. Air Force resolution phantom using 256 detector elements over 270° are depicted in Figure 1. Figure 1A shows the structure of the paper phantom. Figure 1B is the image reconstructed by the L2-norm scheme, whereas Figure 1C indicates the reconstructed image obtained by the ASALSA method. Both reconstructions result in similar initial pressure rise distribution. In contrast, the proposed ASALSA-CF method achieves sharper structure and lesser background artifacts compared with other methods. The artifact reduction is apparent from the zoomed-in areas shown in the insets of Figure 1B-D, and the zoomed-in region is indicated by a red rectangle in Figure 1A. Although the image intensity of line features (label 1 marked in Figure 1B) is partially distorted, the line profiles along the red-dashed line indicated in Figure 1B (shown in Figure 1, E and F) suggest that line features are better resolved in the image reconstructed by the proposed method. Also the first line in label 1 was reconstructed very well using the L2-norm and SALSA method, but it was not accurately reconstructed using the proposed method.
Results of the limited data situation (128 positions over 270°) are depicted in Figure 2. Figure 2A shows the image reconstructed with L2-norm-based algorithm. Because of limited data, the line features (labels 1 and 2) are heavily blurred. However, the image obtained by the ASALSA method recovers the absorbing features much better than that obtained by the L2-norm-based approach. Figure 2C displays the result obtained by the proposed method (ASALSA-CF), where the line features are better resolved as observed in the zoomed-in areas. Analogous to Figure 1, Figure 2C shows fewer artifacts compared with the results reconstructed by other methods. The lateral and axial line profiles (Figure 2, D and E) marked by the red lines in Figure 1B also suggest significant resolution improvement achieved by the proposed reconstruction method on highly undersampled data.
Underdamped data with limited-view condition (128 transducer positions over 135°) are reconstructed, and the corresponding results are shown in Figure 3. The L2-norm-based reconstruction is fully distorted and blurred. Line features (labels 1 and 2) in Figure 3A cannot be identified. In contrast, the ASALSA method shows better performance in resolving the line pattern. Clearly, both images contain artifacts and blurry regions. However, Figure 3C shows that there are fewer artifacts and line features are much better resolved compared with the results of the other reconstruction methods, indicating the superiority of the proposed scheme. Line profiles and zoomed-in images present similar resolution improvement as in previous cases. Meanwhile, the CNR values of line features (yellow labels 1, 2, 3, and 4) are calculated and shown in Table 2. It can be seen that the ASALSA-CF method achieves better image contrast than other methods.
|Methods||Object 1||Object 2||Object 3||Object 4|
The reconstruction results pertaining to the tissue-mimicking agar phantom containing background optical absorption and scattering are shown in Figure 4. Optoacoustic signals from 128 detector positions over 180° are used for reconstruction. The light absorption takes place throughout the phantom, resulting in reflected waves (due to the mismatch between air and tissue-mimicking agar); thus, more reconstruction artifacts are produced. The L2-norm and ASALSA results contain obvious artifacts (white arrows). However, the CF weighting method removes background artifacts, and the 2 absorbing areas are recovered with higher contrast compared with other reconstruction approaches.
The results pertaining to the ex vivo kidney experiment reconstructed from 256 elements over 270° are presented in Figure 5. Figure 5, A and B shows images obtained with the L2-norm and ASALSA method. Analogous to the paper phantom, these 2 images display similar image quality. However, the CF method further improves the reconstruction performance of the SALSA scheme, as illustrated in Figure 5C, showing improved reconstruction quality. Specifically, blood vessel structures marked with the box indicated in Figure 5A are better distinguishable and less blurry with the ASALSA-CF approach compared with the results obtained using other schemes (insets of Figure 5A-C). The visual evaluation is further corroborated by the line profile (Figure 5D) drawn over a given image segment [indicated by the dash line in Figure 5A], indicating that blood vessels marked by the red line are better resolved in the ASALSA-CF reconstructions.
The results for the ex vivo kidney data reconstructed from 128 elements over 135° are presented in Figure 6. Figure 6, A and B show images obtained with the L2-norm and ASALSA method, respectively. It can be seen that the ASALSA method is able to reconstruct more details while producing more artifacts compared with the L2-norm-based scheme. The reconstruction result (Figure 6C) obtained by the proposed method significantly reduces artifacts. Insets in Figure 6A-C, marked by the red box and the line profile (Figure 6D), clearly show that the vessel structures are better distinguishable and have higher contrast in the ASALSA-CF image than those in the images obtained using other methods.
The comparison of different reconstruction schemes with respect to the computational time and memory requirements is presented in Table 3. We calculated the reconstruction time and memory usage for Figures 5 and 6 using a normal PC (Intel Core i5-3470 @2.3GHz and 16 GB memory). It can be seen from Table 3 that the proposed method takes more time and memory compared with the L2-norm approach. However, the conventional SALSA method cannot reconstruct the 256 signals because of computer memory limitation. For 128 signals, the original SALSA method is over 20 times slower and takes 7 times more memory compared with the proposed method.
In this work, we proposed a fast sparse recovery method along with CF weighting for optoacoustic tomographic image reconstruction. The interpolated model matrix method uses a sparse matrix; hence, it is beneficial to use mathematical tools pertaining to sparse algebra. Therefore, the original Basis Pursuit (solved using Augmented Lagrangian method) is rewritten using iterative Krylov subspace solvers framework (LSQR inversion), which tends to converge in fewer iterations. It has been proved that the accelerated SALSA approach can save enormous memory and significantly accelerate the computation time compared with the original SALSA approach.
Previous work used the SALSA algorithm in the deconvolution framework to remove the blurring caused by the regularization parameter (34). In this work, we directly used the SALSA algorithm for performing the image reconstruction and hence avoiding the 2-step procedure. Further, many L1-norm-based algorithms are present in the literature, namely, 2-step iterative shrinkage thresholding, fast iterative shrinkage thresholding algorithm, optimization based on majorization-minimization, and greedy algorithms like orthogonal matching pursuit. SALSA is known as the fastest converging algorithm among all these. Hence, in this work, the SALSA algorithm has been applied and rewritten with the LSQR inversion, foreseeing its utility for real-time implementation.
Table 3 also shows that the proposed scheme is slower compared with the traditional L2-norm-based approach by about 4 times. The reason is that the ASALSA-based approach needs 2 iterative inversion (LSQR) operations (as shown in Table 1). In terms of the order of computation, the ASALSA approach is O(4(M + N)N), while the L2-based reconstruction is O(2(M + N)N), where M is the number of measurements and N is the number of pixels to be reconstructed. This drawback can be overcome by the using graphics processing units to accelerate the reconstruction procedure.
We further hypothesized that the application of a CF weighting will reduce the noise and artifacts that arise during sparsity-based reconstruction, as noise and artifacts are incoherent, whereas optoacoustic signals are coherent. This hypothesis was motivated by the use of CF in ultrasonography for similar reasons (28, 29). CF weighting was integrated into the SALSA algorithm for reducing artifacts arising in limited-view cases. In contrast, artifacts also arise because of acoustically reflecting materials such as bone and air; previous works have reduced these artifacts by using CF weighting (28, 29). From Figure 4 it is evident that CF weighting clearly reduces artifacts arising because of acoustic mismatches. In the resolution phantom measurements, undersampled limited data and limited-view scenarios were retrospectively studied, and it was observed that the ASALSA method outperformed the L2-norm-based method in reconstructing the line features. The CF method further helps in reducing artifacts and improving resolution and contrast. However, the distorted line feature (label 1) in Figure 1 suggests that the CF method may underestimate the image intensity when SNR is low. This can be improved by using the SNR-based CF calculation method (28). Analogous to the phantom measurements, the reconstruction results of the mouse kidney data also prove that the proposed method can better recover vessel structure compared with conventional methods. Overall, the implication of the proposed reconstruction method could be in its utility for limited-view data sets compared with conventional model-based methods and much faster reconstruction than original sparse methods.