Introduction
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 micrometertomillimeter 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 modelbased 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. Modelbased 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, modelbased methods require large numbers of repeated sparse matrixvector multiplications in an iterative manner, resulting in significant computational cost (12, 13). Accelerated modelbased 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 realtime, modelbased reconstruction (15, 16).
A particular challenge in optoacoustic tomography is the implementation of limitedview projections, that is, cases where 360° projections are not available. This could be the case, for example, in imaging large volumes (wholeanimals or humans), whereby access is afforded only from 1 side of the tissue, analogous to ultrasound imaging. Limitedview implementations typically result in lower image fidelity and a larger number of artifacts compared with 360° view data sets (12). Nevertheless, sparsitybased algorithms showed better performance with limitedview data sets (12, 17, 18) compared with Tikhonovbased reconstructions, albeit at a higher computational cost. Moreover, sparse recoverybased methods may amplify noise in limiteddata scenarios (19).
In this work, we propose improvements to sparsitybased 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.
Theory
Acoustic Forward Problem.
The generation and propagation of the acoustic wave is given by the following wave equation (20, 21):
where
p(
r,
t) indicates the acoustic pressure at a position
r and time
t, and
H(
r,
t) indicates the heating function, which is obtained as a product of absorption coefficient and light fluence. β is the thermal expansion coefficient of the tissue, and
C_{p} is the specific heat at constant pressure.
k_{a} represents the acoustic wave number given as
k_{a} = ω/
v_{}, where ω is the temporal frequency and
v_{} is the speed of sound. Note that the heating function is independent, both spatially and temporally, that is,
H(
r,
t) =
H_{r}(
r)
H_{t}(
t). The solution for the above equation is given as follows:
whereby R = r  r′ and spherical integration are performed over surface element
dA′. This integral is discretized to form a model matrix using an interpolated model matrix method to result in the following matrix equation:
where
b is the recorded data and
x is the reconstruction image.
A is obtained by the linear interpolation of the heating function over the image grid.
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:
whereby ∥ ∥
_{2} represents the L2 norm. The above objective function can be solved using normal equations (
22), as follows:
However, Equation 5 is computationally expensive because of the timeconsuming matrix inversion. Alternatively, an LSQR approach can be used (23) as follows:
whereby
B_{k} represents a bidiagonal matrix,
Vk is the right orthogonal matrix resulting from Lanczos bidiagonalization (
23,
24), and β
_{0} is defined as ∥b∥
^{2}_{2}. e
_{1} 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 limiteddata scenarios. Equation 7 is minimized using the SALSA scheme, which has demonstrated the fastest convergence among existing sparsity normbased 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:
where α represents the regularization parameter (depends on the noise). Equation 8 is solved using a
maximum a posterioribased algorithm to obtain an estimate for initial pressure rise (
x). Equation 9 is minimized to obtain an estimate for
v, using a soft thresholding operation (which acts as a derivative for sparsity minimization). The update for the alternated direction method of multiplier parameter is given as
dk+1=dk−(xk+1−vk+1). The minimization in Equations 8 and 9 and the alternated direction method of multiplier parameter update is repeated until convergence.
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 L1normbased 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 limitedview 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:
where
N represents the total number of pixels in the reconstructed domain.
x_{recon} is the reconstructed image obtained using ASALSA and X
_{back} is the backprojection reconstruction. The numerator in Equation 10 represents the energy of the coherent sum of the signals, and the denominator is the total energy sum. The CF values can be interpreted as a focusing quality index estimated from the measured optoacoustic data, ranging from 0 to 1. It is maximal when all signals emitted by an optoacoustic absorber at position r' arrive in the same phase at different detector positions r. After being projected, real signals will constructively superimpose on their point of origin. In this way, good focusing properties can be achieved and consequently a sharp reconstruction. In contrast, incoherent signals will not superimpose on their point of origin after summation, but rather smear out, overall resulting into a degradation of the image quality. Therefore, weighting the amplitude of each image pixel with the corresponding CF can suppress contributions from incoherent signals, enabling identification of noise/artifacts in the reconstructed image and consequently thresholding them. Therefore, the CF is further used for weighting the reconstructed image given as follows:
Table 1.
Aim: Estimation of x by solving Equation 7

Input: A, b, λ, α, max_iter 
Output: x

Initialize ADMM parameter d = 0. 
1. Calculate backprojection solution (x = A^{T}b), CF = A^{T} b^{2}

for k = 1, 2, … max_iter

2. Optimize Equation 9: v = soft(x + d, (0.5 × λ)/α) 
3. Optimize Equation 8: x_{1} = Sk((CkTCk + λIk)−1β1CkTe1), where matrix C is obtained by Lanczos bidiagonalization of A^{T} with k iterations and β1 is the euclidean norm of (A^{T}b + α* (v − d)). 
x_{recon} = Mk((HkTHk + λIk)−1β1HkTe1), where β2 is the Euclidean norm of x_{1}. 
4. Update ADMM parameter: d = d + x_{recon} − v 
end

5. CF = (x^{2}/CF)

6. x = x_{recon}. × 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.
Experimental Methods
Optoacoustic Imaging System
A multispectral optoacoustic tomography small animal scanner (MSOT256TF, iThera Medical GmbH, Munich, Germany) was used to experimentally examine the performance of the proposed reconstruction method. In brief, a custommade, 256element 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 wavelengthtunable 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 signaltonoise 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.9cmdiameter 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 limitedview scenarios, we assumed 2 downsampled data sets, one using 128 positions over 270° coverage angle and one using 128 positions over 135° coverage angle.
Figure 1.
Reference U.S. Air Force phantom printed on white paper with black ink, which was embedded in scattering agar (A). The reconstructed image by L2norm (B). The ASALSA method (C) and the proposed method (D). The subsects in (B–D) are the zoomedin regions marked in a red rectangle of (A). The line profiles in the horizontal and vertical directions marked in (B) are represented in (E) and (F), respectively.
A tissuemimicking agar phantom was also prepared to examine the ASALSACF performance. The phantom included areas containing 0.016% India Ink to impart higher optical absorption than the background. The absorption coefficient of the inkcontaining 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.
Image Reconstruction
The performance of reconstructions based on Equation 6 (L2norm), Equations 8 and 9 (ASALSA), and Equation 11 (ASALSACF) was compared using the experimental measurements collected from phantoms and mouse kidneys. Before reconstruction, the optoacoustic signals were bandpass filtered, with cutoff frequencies between 0.2 and 7 MHz to remove lowfrequency offsets and highfrequency 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 L2normbased scheme was obtained using an Lcurve 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).
Results
The reconstruction results corresponding to the printedpaper 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 L2norm 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 ASALSACF method achieves sharper structure and lesser background artifacts compared with other methods. The artifact reduction is apparent from the zoomedin areas shown in the insets of Figure 1BD, and the zoomedin 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 reddashed 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 L2norm 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 L2normbased 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 L2normbased approach. Figure 2C displays the result obtained by the proposed method (ASALSACF), where the line features are better resolved as observed in the zoomedin 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.
Figure 2.
Images reconstructed using 128 transducer positions over 270 degrees. The reconstructed image by L2norm (A). The ASALSA method (B) and the proposed method (D). The subsects in (A–C) are the zoomedin regions marked in a red rectangle of Figure 1A. The line profiles in the horizontal and vertical directions marked in Figure 1B are represented in (D) and (E), respectively.
Underdamped data with limitedview condition (128 transducer positions over 135°) are reconstructed, and the corresponding results are shown in Figure 3. The L2normbased 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 zoomedin 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 ASALSACF method achieves better image contrast than other methods.
Figure 3.
Images reconstructed using 128 transducer positions over 135 degrees. The reconstructed image by L2norm (A). ASALSA (B) and the proposed method (D). The subsects in (A–C) are the zoomedin regions marked in a red rectangle of Figure 1A. The line profiles in the horizontal and vertical directions marked in Figure 1A are represented in (D) and (E), respectively.
Table 2.
Contrast (CNR) Comparison
Methods 
Object 1 
Object 2 
Object 3 
Object 4 
L2norm 
0.1(D^{a}) 
0.4(D) 
0.2(D) 
1.4 
SALSA 
0.9 
0.7 
1.2 
2.1 
ASALSACF 
1.2 
1.4 
2.4 
3.7 
The reconstruction results pertaining to the tissuemimicking 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 tissuemimicking agar); thus, more reconstruction artifacts are produced. The L2norm 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.
Figure 4.
Reconstructed images of the tissuemimicking agar phantom, which includes a hollow cavity filled with air and 2 high absorbing areas. Reference image of the phantom (A). The reconstructed image by L2norm (B). The ASALSA method (C) and the proposed method (D). Arrows indicate artifacts caused by reflections or scattering of the acoustic waves, which are significantly reduced with the proposed method.
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 L2norm 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 ASALSACF approach compared with the results obtained using other schemes (insets of Figure 5AC). 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 ASALSACF reconstructions.
Figure 5.
Reconstructed images of the mouse kidney from 256 transducer positions over 270 degrees. The reconstructed image by L2norm (A). The ASALSA method (B) and the proposed method (C). The subsects in (A–C) are the zoomedin regions marked in a red rectangle of Figure 5A. The line profiles marked by the red line in Figure 5A are represented in (D).
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 L2norm 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 L2normbased scheme. The reconstruction result (Figure 6C) obtained by the proposed method significantly reduces artifacts. Insets in Figure 6AC, 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 ASALSACF image than those in the images obtained using other methods.
Figure 6.
Reconstructed images of the mouse kidney from 128 signal positions over 135 degrees. The reconstructed image by L2norm (A). The ASALSA method (B) and the proposed method (C). The subsects in (A–C) are the zoomedin regions marked in a red rectangle of Figure 6A. The line profiles marked by the red line in Figure 6A are represented in (D).
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 i53470 @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 L2norm 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.
Table 3.
Computation Time(s) and Memory (GB)
Methods 
Figure 5

Figure 6

L2norm 
79/2.5 
35/1.3 
SALSA 
Out of memory 
3554/14.7 
ASALSA_CF 
294/4.7 
137/2.2 
Discussion
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 2step procedure. Further, many L1normbased algorithms are present in the literature, namely, 2step iterative shrinkage thresholding, fast iterative shrinkage thresholding algorithm, optimization based on majorizationminimization, 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 realtime implementation.
Table 3 also shows that the proposed scheme is slower compared with the traditional L2normbased approach by about 4 times. The reason is that the ASALSAbased 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 L2based 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 sparsitybased 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 limitedview 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 limitedview scenarios were retrospectively studied, and it was observed that the ASALSA method outperformed the L2normbased 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 SNRbased 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 limitedview data sets compared with conventional modelbased methods and much faster reconstruction than original sparse methods.
Acknowledgments
Hailong He and Jaya Prakash contributed equally to this work.
Andreas Buehler acknowledges funding from the European Union project FAMOS (FP7 ICT, Contract No. 317744). Jaya Prakash acknowledges Alexander von Humboldt Postdoctoral Fellowship program. Jaya Prakash also acknowledges initial discussion with Dr. Calvin B. Shaw. The authors acknowledge the second phantom data from Dr. X. Luís DeánBen.
Disclosure: No disclosures to report.
Conflict of Interest: None reported.
References

Taruttis A, van Dam GM, Ntziachristos V. Mesoscopic and macroscopic optoacoustic imaging of cancer. Cancer Res. 2015;75(8):1548–1559.

Taruttis A, Ntziachristos V. Advances in realtime multispectral optoacoustic imaging and its applications. Nat Photonics. 2015;9(4):219–227.

Beard P. Biomedical photoacoustic imaging. Interface Focus. 2011;1(4):602–631.

Ntziachristos V, Razansky D. Molecular imaging by means of multispectral optoacoustic tomography (MSOT). Chem Rev. 2010;110(5):2783–2794.

Xu MH, Wang LHV. Photoacoustic imaging in biomedicine. Rev Sci Instrum. 2006;77(4):041101.

DeanBen XL, Razansky D, Ntziachristos V. The effects of acoustic attenuation in optoacoustic signals. Phys Med Biol. 2011;56(18):6129–6148.

Xu M, Wang L. Universal backprojection algorithm for photoacoustic computed tomography. Phys Rev E Stat Nonlin Soft Matter Phys. 2005;71(1 Pt 2):016706.

Paltauf G, Viator JA, Prahl SA, Jacques SL. Iterative reconstruction algorithm for optoacoustic imaging. J Acoust Soc Am. 2002;112(4):1536–1544.

Rosenthal A, Ntziachristos V, Razansky D. Acoustic inversion in optoacoustic tomography: A review. Curr Med Imaging Rev. 2013;9(4):318–336.

Caballero MAA, Rosenthal A, Buehler A, Razansky D, Ntziachristos V. Optoacoustic determination of spatiotemporal responses of ultrasound sensors. IEEE Trans Ultrason Ferroelectr Freq Control. 2013;60(6):1234–1244.

Rosenthal A, Ntziachristos V, Razansky D. Optoacoustic methods for frequency calibration of ultrasonic sensors. IEEE Trans Ultrason Ferroelectr Freq Control. 2011;58(2):316–326.

Buehler A, Rosenthal A, Jetzfellner T, Dima A, Razansky D, Ntziachristos V. Modelbased optoacoustic inversions with incomplete projection data. Med Phys. 2011;38(3):1694–1704.

DeanBen XL, Ntziachristos V, Razansky D. Acceleration of optoacoustic modelbased reconstruction using angular image discretization. IEEE Trans Med Imaging. 2012;31(5):1154–1162.

Wang K, Xia J, Li C, Wang LV, Anastasio MA. Fast spatiotemporal image reconstruction based on lowrank matrix estimation for dynamic photoacoustic computed tomography. J Biomed Opt. 2014;19(5):056007.

Ding L, DeanBen XL, Razansky D. Realtime Modelbased inversion in crosssectional optoacoustic tomography. IEEE Trans Med Imaging. 2016. [Epub ahead of print].

DeanBen XL, Ozbek A, Razansky D. Volumetric realtime tracking of peripheral human vasculature with GPUaccelerated threedimensional optoacoustic tomography. IEEE Trans Med Imaging. 2013;32(11):2050–2055.

Han Y, Tzoumas S, Nunes A, Ntziachristos V, Rosenthal A. Sparsitybased acoustic inversion in crosssectional multiscale optoacoustic imaging. Med Phys. 2015;42(9):5444–5452.

Provost J, Lesage F. The application of compressed sensing for photoacoustic tomography. IEEE Trans Med Imaging. 2009;28(4):585–594.

Lee K, Bresler Y, Junge M. Subspace Methods for Joint Sparse Recovery. IEEE Trans Inf Theory. 2012;58(6):3613–41.

Xu MH, Xu Y, Wang LHV. Timedomain reconstructionalgorithms and numerical simulations for thermoacoustic tomography in various geometries. IEEE Trans Biomed Eng. 2003;50(9):1086–1099.

Xu M, Wang L. Analytic explanation of spatial resolution related to bandwidth and detector aperture size in thermoacoustic or photoacoustic reconstruction. Phys Rev E Stat Nonlin Soft Matter Phys. 2003;67(5 Pt 2):056605.

Calvetti D, Morigi S, Reichel L, Sgallari F. Tikhonov regularization and the Lcurve for large discrete illposed problems. J Comput Appl Math. 2000;123(12):423–446.

Saunders CCPaMA. LSQR: An algorithm for sparse linearequations and sparse leastsquares. ACM Trans Math Softw. 1982;8(1):43–71.

Prakash J, Yalavarthy PK. A LSQRtype method provides a computationally efficient automated optimal choice of regularization parameter in diffuse optical tomography. Med Phys. 2013;40(3):033101.

Lustig M, Donoho D, Pauly JM. Sparse MRI: The application of compressed sensing for rapid MR imaging. Magn Reson Med. 2007;58(6):1182–1195.

Prakash J, Shaw CB, Manjappa R, Kanhirodan R, alavarthy PK. Sparse recovery methods hold promise for diffuse optical tomographic image reconstruction. IEEE J Sel Top Quantum Electron. 2014;20(2):74–82.

Afonso MV, BioucasDias JM, Figueiredo MAT. Fast image recovery using variable splitting and constrained optimization. IEEE Trans Image Process. 2010;19(9):2345–2356.

Wang YH, Li PC. SNRdependent coherencebased adaptive imaging for highframerate ultrasonic and photoacoustic imaging. IEEE Trans Ultrason Ferroelectr Freq Control. 2014;61(8):1419–32.

Li PC, Li ML. Adaptive imaging using the generalized coherence factor. IEEE Trans Ultrason Ferroelectr Freq Control. 2003;50(2):128–141.

Razansky D, Buehler A, Ntziachristos V. Volumetric realtime multispectral optoacoustic tomography of biomarkers. Nat Protoc. 2011;6(8):1121–1129.

Buehler A, Herzog E, Razansky D, Ntziachristos V. Video rate optoacoustic tomography of mouse kidney perfusion. Opt Lett. 2010;35(14):2475–2477.

Dima A, Burton NC, Ntziachristos V. Multispectral optoacoustic tomography at 64, 128, and 256 channels. J Biomed Opt. 2014;19(3):36021.

Mandal S, Nasonova E, DeanBen XL, Razansky D. Optimal selfcalibration of tomographic reconstruction parameters in wholebody small animal optoacoustic imaging. Photoacoustics. 2014;2(3):128–136.

Prakash J, Raju AS, Shaw CB, Pramanik M, Yalavarthy PK. Basis pursuit deconvolution for improving modelbased reconstructed images in photoacoustic tomography. Biomed Opt Express. 2014;5(5):1363–1377.
Research Articles
Download PDF (4.56 MB)
TOMOGRAPHY, June 2016, Volume 2, Issue 2:138145
DOI: 10.18383/j.tom.2016.00148
Optoacoustic Tomography Using Accelerated Sparse Recovery and Coherence Factor Weighting
Hailong He^{1}, Jaya Prakash^{1}, Andreas Buehler^{1}, Vasilis Ntziachristos^{1}
Abstract
Sparse recovery algorithms have shown great potential to accurately reconstruct images using limitedview optoacoustic (photoacoustic) tomography data sets, but these are computationally expensive. In this paper, we propose an improvement of the fast converging Split Augmented Lagrangian Shrinkage Algorithm method based on least square QR inversion for improving the reconstruction speed. We further show image fidelity improvement when using a coherence factor to weight the reconstruction result. Phantom and in vivo measurements show that the accelerated Split Augmented Lagrangian Shrinkage Algorithm method with coherence factor weighting offers images with reduced artifacts and provides faster convergence compared with existing sparse recovery algorithms.
Introduction
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 micrometertomillimeter 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 modelbased 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. Modelbased 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, modelbased methods require large numbers of repeated sparse matrixvector multiplications in an iterative manner, resulting in significant computational cost (12, 13). Accelerated modelbased 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 realtime, modelbased reconstruction (15, 16).
A particular challenge in optoacoustic tomography is the implementation of limitedview projections, that is, cases where 360° projections are not available. This could be the case, for example, in imaging large volumes (wholeanimals or humans), whereby access is afforded only from 1 side of the tissue, analogous to ultrasound imaging. Limitedview implementations typically result in lower image fidelity and a larger number of artifacts compared with 360° view data sets (12). Nevertheless, sparsitybased algorithms showed better performance with limitedview data sets (12, 17, 18) compared with Tikhonovbased reconstructions, albeit at a higher computational cost. Moreover, sparse recoverybased methods may amplify noise in limiteddata scenarios (19).
In this work, we propose improvements to sparsitybased 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.
Theory
Acoustic Forward Problem.
The generation and propagation of the acoustic wave is given by the following wave equation (20, 21):
(1)
(2)
(3)
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:
(4)
(5)
However, Equation 5 is computationally expensive because of the timeconsuming matrix inversion. Alternatively, an LSQR approach can be used (23) as follows:
(6)
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:
(7)
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 limiteddata scenarios. Equation 7 is minimized using the SALSA scheme, which has demonstrated the fastest convergence among existing sparsity normbased 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:
(8)
(9)
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 L1normbased 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 limitedview 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:
(10)
(11)
Table 1.
ASALSACF Algorithm
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.
Experimental Methods
Optoacoustic Imaging System
A multispectral optoacoustic tomography small animal scanner (MSOT256TF, iThera Medical GmbH, Munich, Germany) was used to experimentally examine the performance of the proposed reconstruction method. In brief, a custommade, 256element 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 wavelengthtunable 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 signaltonoise 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.9cmdiameter 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 limitedview scenarios, we assumed 2 downsampled data sets, one using 128 positions over 270° coverage angle and one using 128 positions over 135° coverage angle.
Figure 1.
Reference U.S. Air Force phantom printed on white paper with black ink, which was embedded in scattering agar (A). The reconstructed image by L2norm (B). The ASALSA method (C) and the proposed method (D). The subsects in (B–D) are the zoomedin regions marked in a red rectangle of (A). The line profiles in the horizontal and vertical directions marked in (B) are represented in (E) and (F), respectively.
A tissuemimicking agar phantom was also prepared to examine the ASALSACF performance. The phantom included areas containing 0.016% India Ink to impart higher optical absorption than the background. The absorption coefficient of the inkcontaining 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.
Image Reconstruction
The performance of reconstructions based on Equation 6 (L2norm), Equations 8 and 9 (ASALSA), and Equation 11 (ASALSACF) was compared using the experimental measurements collected from phantoms and mouse kidneys. Before reconstruction, the optoacoustic signals were bandpass filtered, with cutoff frequencies between 0.2 and 7 MHz to remove lowfrequency offsets and highfrequency 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 L2normbased scheme was obtained using an Lcurve 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).
Results
The reconstruction results corresponding to the printedpaper 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 L2norm 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 ASALSACF method achieves sharper structure and lesser background artifacts compared with other methods. The artifact reduction is apparent from the zoomedin areas shown in the insets of Figure 1BD, and the zoomedin 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 reddashed 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 L2norm 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 L2normbased 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 L2normbased approach. Figure 2C displays the result obtained by the proposed method (ASALSACF), where the line features are better resolved as observed in the zoomedin 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.
Figure 2.
Images reconstructed using 128 transducer positions over 270 degrees. The reconstructed image by L2norm (A). The ASALSA method (B) and the proposed method (D). The subsects in (A–C) are the zoomedin regions marked in a red rectangle of Figure 1A. The line profiles in the horizontal and vertical directions marked in Figure 1B are represented in (D) and (E), respectively.
Underdamped data with limitedview condition (128 transducer positions over 135°) are reconstructed, and the corresponding results are shown in Figure 3. The L2normbased 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 zoomedin 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 ASALSACF method achieves better image contrast than other methods.
Figure 3.
Images reconstructed using 128 transducer positions over 135 degrees. The reconstructed image by L2norm (A). ASALSA (B) and the proposed method (D). The subsects in (A–C) are the zoomedin regions marked in a red rectangle of Figure 1A. The line profiles in the horizontal and vertical directions marked in Figure 1A are represented in (D) and (E), respectively.
Table 2.
Contrast (CNR) Comparison
i] D^{a}: Distorted
The reconstruction results pertaining to the tissuemimicking 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 tissuemimicking agar); thus, more reconstruction artifacts are produced. The L2norm 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.
Figure 4.
Reconstructed images of the tissuemimicking agar phantom, which includes a hollow cavity filled with air and 2 high absorbing areas. Reference image of the phantom (A). The reconstructed image by L2norm (B). The ASALSA method (C) and the proposed method (D). Arrows indicate artifacts caused by reflections or scattering of the acoustic waves, which are significantly reduced with the proposed method.
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 L2norm 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 ASALSACF approach compared with the results obtained using other schemes (insets of Figure 5AC). 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 ASALSACF reconstructions.
Figure 5.
Reconstructed images of the mouse kidney from 256 transducer positions over 270 degrees. The reconstructed image by L2norm (A). The ASALSA method (B) and the proposed method (C). The subsects in (A–C) are the zoomedin regions marked in a red rectangle of Figure 5A. The line profiles marked by the red line in Figure 5A are represented in (D).
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 L2norm 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 L2normbased scheme. The reconstruction result (Figure 6C) obtained by the proposed method significantly reduces artifacts. Insets in Figure 6AC, 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 ASALSACF image than those in the images obtained using other methods.
Figure 6.
Reconstructed images of the mouse kidney from 128 signal positions over 135 degrees. The reconstructed image by L2norm (A). The ASALSA method (B) and the proposed method (C). The subsects in (A–C) are the zoomedin regions marked in a red rectangle of Figure 6A. The line profiles marked by the red line in Figure 6A are represented in (D).
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 i53470 @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 L2norm 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.
Table 3.
Computation Time(s) and Memory (GB)
Discussion
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 2step procedure. Further, many L1normbased algorithms are present in the literature, namely, 2step iterative shrinkage thresholding, fast iterative shrinkage thresholding algorithm, optimization based on majorizationminimization, 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 realtime implementation.
Table 3 also shows that the proposed scheme is slower compared with the traditional L2normbased approach by about 4 times. The reason is that the ASALSAbased 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 L2based 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 sparsitybased 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 limitedview 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 limitedview scenarios were retrospectively studied, and it was observed that the ASALSA method outperformed the L2normbased 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 SNRbased 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 limitedview data sets compared with conventional modelbased methods and much faster reconstruction than original sparse methods.
Notes
[2] Abbreviations:
SALSA
Split Augmented Lagrangian Shrinkage algorithm
ASALSACF
SALSA method with coherence factor
CF
coherence factor
GPU
graphics processing unit
LSQR
least square QR
ASALSA
accelerated SALSA
Acknowledgments
Hailong He and Jaya Prakash contributed equally to this work.
Andreas Buehler acknowledges funding from the European Union project FAMOS (FP7 ICT, Contract No. 317744). Jaya Prakash acknowledges Alexander von Humboldt Postdoctoral Fellowship program. Jaya Prakash also acknowledges initial discussion with Dr. Calvin B. Shaw. The authors acknowledge the second phantom data from Dr. X. Luís DeánBen.
Disclosure: No disclosures to report.
Conflict of Interest: None reported.
References
Journal Information
Journal ID (nlmta): tom
Journal ID (publisherid): TOMOG
Title: Tomography
Subtitle: A Journal for Imaging Research
Abbreviated Title: Tomography
ISSN (print): 23791381
ISSN (electronic): 2379139X
Publisher: Grapho Publications, LLC (Ann Abor, Michigan)
Article Information
Copyright statement: © 2016 The Authors. Published by Grapho Publications, LLC
Copyright: 2016, Grapho Publications, LLC
License (openaccess, http://creativecommons.org/licenses/byncnd/4.0/):
This is an open access article under the CC BYNCND license (http://creativecommons.org/licenses/byncnd/4.0/).
Publication date (print): June 2016
Volume: 2
Issue: 2
Pages: 138145
Publisher ID: TOM0014816
DOI: 10.18383/j.tom.2016.00148
PDF
Download the article PDF (4.56 MB)
Download the full issue PDF (76.23 MB)
Mobileready Flipbook
View the full issue as a flipbook (Desktop and Mobileready)