|
1.IntroductionOne major goal of photoacoustic imaging is the recovery of the optical properties of an imaged subject. Of particular interest is optical absorption, which is directly proportional to the received pressures. However, the initial pressure distribution is determined by , where is the Grüneisen parameter, is the optical absorption, and is the fluence. All of these may vary spatially, and to complicate matters, depends on the propagation of light through turbid media with spatially varying absorption and scattering coefficients. For this manuscript, we assume that is recovered faithfully, which is itself nontrivial. Some approaches to that problem include time-reversal,1 filtered backprojection,2 and model-based inversion.3 Analytical solutions for have been attempted,4–6 but the inverse problem is difficult, and the solutions are limited in scope. Iterative approaches show some promise7,8 as a more general solution, but suffer from the potential for overiteration9 and nonuniqueness10 problems. We have introduced iterative techniques based on multiple illumination photoacoustic tomography (MIPAT) to tackle both of these problems.10,11 Fixed-point iterative techniques, like that proposed by Cox et al.7 and experimentally tested by Jetzfellner et al.,9 are quite sensitive to noise and may diverge with overiteration. The difficulty in applying these iterative techniques may be exacerbated by an MIPAT setup using shutters, or where fluence is already at the ANSI safety limit. The regularization parameter common to these techniques can help ensure convergence under noisy conditions, but negatively impacts reconstruction accuracy. An increase in the signal-to-noise ratio (SNR) of the initial pressure images provides more robust reconstruction of the optical absorption parameter by allowing a lower regularization parameter to be used. If increasing laser energy is not an option (due to system setup or safety concerns), then one possible solution is to use patterned illumination. We have recently shown that S-sequences can provide nearly the same SNR gain in synthetic transmit aperture ultrasound imaging as Hadamard encoding,12 but without the requirement of pulse inversion. This means that the technique can readily be adapted to photoacoustic imaging, which can only create an initial pressure distribution. It is possible to approximate Hadamard encoding by using separate pulse sequences for positive and negative elements in the Hadamard matrix and subtracting those intermediate encodings. However, this would require twice as many pulses. An alternative approach to MIPAT techniques uses diffuse optical tomography (DOT) to ameliorate the nonuniqueness problem.8 Potentially, this approach could be combined with our MIPAT approach, and S-sequence encoding could even be applied to DOT. This is, however, beyond the scope of this manuscript. In this work, we apply S-sequences to fixed-point iterative MIPAT and show that the resilience to noise is greatly improved over unencoded imaging. The basic idea is that multiple sources can be used simultaneously and then after a complete set of illumination patterns has been applied, multiplication by an inverse of the encoding matrix can recover images effectively due to a single source, but with enhanced SNR. The enhanced SNR images can then be used for improved quantitative reconstruction. 2.Theory2.1.Fixed-Point Iterative MIPATFixed-point iteration for MIPAT is a fairly straightforward technique which we have previously proposed11 based on the original work by Cox et al.7 Those works may be consulted for a more thorough treatment of the technique. Briefly, the goal is to iteratively reconstruct an estimate of the optical absorption, . To do this, we use sources (where the distribution and type of the sources may be arbitrary, but are known), resulting in reconstructed initial pressure estimates, which we call . It is assumed that these pressure estimates (or images) are proportional to the optical absorption by , where is the Grüneisen parameter (taken to be uniform) and is an estimate of the fluence. We then iteratively update and at each location over several iterations . The algorithm is described as follows:
The above algorithm was previously shown to provide excellent convergence properties over the single-illumination case but has not yet been tested for patterned illumination, which is the topic of this paper. 2.2.S-Sequence EncodingNow we consider the images . The previous assumption that the images are properly reconstructed is greatly complicated by concerns of SNR—our previous study showed that to keep a low , noise levels must be extremely low.11 The simplest way to increase SNR in photoacoustic images is to increase the fluence used for imaging. This may not always be possible due to hardware or safety limitations. The use of multiple illumination sources opens up an opportunity to easily apply spatially encoded illuminations. To accomplish this, we borrow from previous work using encoded optical spectroscopy14 and from our work on synthetic transmit aperture ultrasound.12 In Ref. 12, we introduced an encoding scheme based on S-sequences, which has practical advantages over other encoding schemes. Most importantly, an S-matrix only contains binary digits (0 and 1) and, thus, does not depend on an inverted signal, which is impossible to generate in photoacoustic imaging since we are constrained to positive fluence values. To adapt ultrasound encoding techniques to photoacoustic imaging, we can consider the same linear system representing encoding: . The ’th column of consists of photoacoustic data (for example, raw sinogram data or reconstructed initial pressure data) collected from the ’th illumination pattern and ordered into a column vector by rasterization. The matrix has elements , which are the data element (for example, the ‘th pixel in the reconstructed PAT image) due to individual source . represents the complete collection of image data from each respective source location and is the starting point for our previous work on multiple-illumination PAT. is a source weighting matrix. has columns with elements that select the sources used for the ’th illumination pattern. Thus, received pressures from effective single sources can be recovered by multiplication by the inverse of the encoding matrix . If were the identity matrix, would represent the collection of photoacoustic data from single-illumination sources. When has more than one nonzero element in a column, it means that multiple sources are simultaneously activated within a pulsed excitation. By judicious choice of , it may be possible to activate more than one source per illumination event and, hence, deliver more energy than is possible with a single source, while being able to recover the effective data from effective single illuminations (but with enhanced SNR). The recovered can then be used in the iterative least-squares multiple-illumination tomography algorithm described above.12 While many choices exist for the matrix, previous work has shown that the S-matrix is a good choice because it provides a best possible balance between inversion stability and energy delivery to provide improved SNR.12,14 The S-matrix is derived from a Hadamard matrix by replacing all 1’s with 0’s and all ’s with 1’s and then removing the first row and column. S-sequences are rows or columns of the S-matrix. Using the classic Hadamard construction, this constrains us to an matrix, where , . We, thus, use successive illumination patterns with sources to reconstruct a single image. The S-matrix is easily inverted and can be applied to encode and then decode the sources, resulting in a with potentially up to a improvement in intensity.12 Thus, we add an encoding step by selecting sources according to the S-sequence before imaging, and a decoding step via a matrix multiplication after imaging to produce higher-SNR before using the iterative technique. The encoding technique is similar to how S-sequences are applied in ultrasound imaging, but the apodization is used to select multiple sources. The images can then be formed into a matrix as described above and multiplied by to recover single-source images. Since all illumination sources and patterns are known, it is possible to directly use the nondecoded images (i.e., instead of ) in the finite element solver to save some computation time. However, the problem becomes much less well-posed since each image will use sources. This results in an illumination that appears much more uniform, especially deep in scattering tissue. Previous analytic work by Bal and Ren similarly stresses the importance of illumination choice through constraints on a vector field comparing two illuminations.15 We, therefore, expect that the iterative technique will underperform if the decoding step is not included. 3.SimulationsWe begin by defining a phantom identical to that used in our previous work: an ellipse of by 1.4 cm, where , with an inclusion of by 0.44 cm, where in a nonabsorbing background media. Over the entire field, . We consider a pixel area, representing a field with incident ring illumination approximated by 315 point sources evenly distributed and located one transport mean free path inside the phantom. These sources are simulated individually using the TOAST forward solver (using with Dirichlet boundary conditions) and recombined to form a portion of a ring illumination according to , 7, 15, 63. These intermediate illuminations (single-source illuminations) are used to produce an initial pressure (i.e., simulated photoacoustic image) by multiplying by a constant and the known . In the case of S-sequence encoding, single source images are combined into S-sequence images by adding images selected according to the S-matrix. Noise is added to these images at a level equivalent to the SNR for the ring illumination case (i.e., using the maximum of the sum of image data from all sources to scale Gaussian noise to effectively reach the SNR when all sources are used). This is equivalent to selectively blocking unused sources in an experimental setup, while keeping all other factors the same. This situation might be appropriate when considering ANSI-limited exposure for each illumination. We then apply the previously described iterative technique and track the normalized root-mean-squared error defined as , where . 3.1.Image QualityThe first important question to answer is whether or not S-sequence coding confers an SNR advantage over unencoded imaging on a per-image basis. We look at three different scenarios: unencoded images (using multiple single-source illuminations), nondecoded images (using S-sequence patterned illumination), and images that have been decoded after encoding. Since we are using simulations, we can exactly model the reconstructed image free of noise. This is most important for the decoded signal, where the noise signal will be some linear combination of noise signals from all images and, thus, the noise signal is not directly available at the time of image formation. SNR can then be characterized for each initial pressure distribution as follows: . Figure 1 illustrates the alternative MIPAT techniques and shows a qualitative SNR increase. Unencoded MIPAT uses illuminations from a single source as in the iterative technique and provides an absorption estimate that is sensitive to SNR. As a general rule, when more sources are used, the unencoded images lose SNR since an increasing amount of the incident fluence is discarded. The nondecoded images always use about half the available energy and, thus, the SNR improvement approaches the added SNR level. One might expect that SNR should be slightly less than the target level, which is added based on a uniform illumination, but our definition of SNR using the maximum value is easily achieved using only a few neighboring sources since sources on opposite sides of the phantom contribute little to their respective largest pixel value. Using these images directly as is possible, but results in poor reconstruction, motivating the decoding step. The nondecoded images result in poor reconstruction because they are similar to the uniform illumination case previously shown to diverge with overiteration. The decoded images exhibit very good agreement with the SNR of the nondecoded images over a wide range of included noise variances, while the unencoded images show a reduction in SNR appropriate for the number of sources. For example, in the 15 illumination case, we see an average SNR increase of , which is in very good agreement with the predicted increase of 10.5 dB. 3.2.Fixed-Point Iterative MIPATWe next consider the algorithm performance in the same three cases: using unencoded, nondecoded, or decoded images in the iterative technique. Figure 2 illustrates the behavior of the technique for the three examined scenarios over 30 iterations. While there is some similarity between the unencoded and decoded scenarios, it is clear that simply using the patterned illumination does not yield good results. Figure 3 illustrates the performance over many SNR conditions. In particular, the case is quite interesting. To ensure convergence, the unencoded images require an SNR that increases somewhat with , requiring added SNR in the case. Running the algorithm with the nondecoded images actually provides poor reconstruction even in the case, which is not surprising given the similarity in terms of incident fluence for the incident patterned illuminations. It may be difficult to tell from the figures, but the case performs quite a bit worse in the nondecoded case than the other two. Finally, when the decoded images are used, the required SNR for convergence in the case appears to be independent of at . Figure 4 shows a closer view of the unencoded and decoded views. The case does appear to favor a less accurate reconstruction, though the other two cases have similar behavior. An example of the different reconstructions can be seen in Fig. 5 for the different situations. While all figures saturate beyond the maximum true , the most faithful reconstruction is clearly the decoded image in Fig. 5(d). 4.Discussion and ConclusionsBalancing regularization with performance is a difficult task for iterative techniques where single points in the image may introduce numerical instability. Automated approaches to select an optimal regularization parameter have been proposed for the estimation of initial pressure distribution and may prove important for quantitative techniques.16,17 Starting with good initial pressure distributions with excellent SNR is one important factor for reducing in our iterative multiple-illumination technique. By applying S-sequence coding, we can recover higher-quality images through decoding, as shown in Fig. 1. We find that the nondecoded and decoded images have similar SNRs and exhibit close to the expected theoretical SNR increase. Figure 2 gives a sense of the speed and general trends of convergence, with the decoded case appearing to converge a bit faster, and possibly forcing convergence in the high case. Figure 3 shows that we can improve convergence in the case by using the decoded images. As expected, the nondecoded images result in poor performance of the iterative technique. When sources are selected according to an S-sequence, the illumination patterns are similar to uniform illumination, thus, we see the same overiteration problem from the experimental work of Jetzfellner et al.9 Figure 4 shows the danger of a that is too high in the case, where a relatively poor solution appears to be favored. This raises the very important question of the selection of . Unfortunately, it is not one that can be answered simply, especially based solely on simulations. Certainly, for a given set of images, there is likely to be a value of that will result in a convergent result, but its value may vary depending on the specifics of the object being imaged, and it may not be clear if the reconstruction is accurate. It is possible, though computationally inefficient, to start with a very low and increase it until numerical instability subsides. Experimental validation and exploration of in phantom and in vivo studies will be required to put this method into practice. The illumination patterns could be experimentally realized by either selectively blocking beams using shutters (which would have the disadvantage of losing light) or by using a steerable mirror array perhaps in combination with a fiber bundle, where each branch of a fiber-bundle output could be directed to one of two locations. Finally, Fig. 5 shows an example reconstructed compared to the true . It is quite clear that the decoded images provide the most faithful reconstruction of the phantom. One thing that is not immediately clear from viewing the figure is that the unencoded and nondecoded reconstructions have small areas in the inclusion that exceed the true by orders of magnitude, whereas the use of decoded images results in a more faithful reconstruction. The major caveat with this simple iterative technique remains SNR: if there is insufficient SNR deep in tissue, then additional illuminations will provide no additional constraint on the linear system, thus, the technique will not provide a faithful reconstruction. We have demonstrated that patterned illumination provides a powerful tool to boost SNR in tomographic systems to a level where fixed-point iterative MIPAT should be practical with little averaging required. While S-sequences are used here, any sort of binary coding scheme could be used, though other coding schemes may have a worse condition number, increasing the reconstruction error. Patterned illumination is applicable to any multiple-illumination technique and provides all the advantages of averaging without increasing imaging time. AcknowledgmentsWe gratefully acknowledge funding from NSERC (355544-2008, 375340-2009, STPGP 396444, student scholarships), Terry-Fox Foundation and the Canadian Cancer Society (TFF 019237, TFF 019240, CCS 2011-700718), Alberta Advanced Education & Technology (student scholarships), and China Scholarship Council scholarships for graduate student Peng Shao. ReferencesB. E. TreebyE. Z. ZhangB. T. Cox,
“Photoacoustic tomography in absorbing acoustic media using time reversal,”
Inverse Probl., 26
(11), 115003
(2010). http://dx.doi.org/10.1088/0266-5611/26/11/115003 INPEEY 0266-5611 Google Scholar
Y. ZhangY. Wang,
“An improved filtered back-projection algorithm for photoacoustic tomography,”
in 5th Int. Conf. on Bioinformatics and Biomedical Engineering,
1
–4
(2011). Google Scholar
A. RosenthalD. RazanskyV. Ntziachristos,
“Fast semi-analytical model-based acoustic inversion for quantitative optoacoustic tomography,”
IEEE Trans. Med. Imaging, 29
(6), 1275
–1285
(2010). http://dx.doi.org/10.1109/TMI.2010.2044584 ITMID4 0278-0062 Google Scholar
J. RipollV. Ntziachristos,
“Quantitative point source photoacoustic inversion formulas for scattering and absorbing media,”
Phys. Rev. E, 71
(3), 031912
(2005). http://dx.doi.org/10.1103/PhysRevE.71.031912 PLEEE8 1063-651X Google Scholar
Z. YuanH. Jiang,
“Quantitative photoacoustic tomography: recovery of optical absorption coefficient maps of heterogeneous media,”
Appl. Phys. Lett., 88
(23), 231101
(2006). http://dx.doi.org/10.1063/1.2209883 APPLAB 0003-6951 Google Scholar
B. Banerjeeet al.,
“Quantitative photoacoustic tomography from boundary pressure measurements: noniterative recovery of optical absorption coefficient from the reconstructed absorbed energy map,”
J. Opt. Soc. Am. A, 25
(9), 2347
–2356
(2008). http://dx.doi.org/10.1364/JOSAA.25.002347 JOAOD6 0740-3232 Google Scholar
B. T. Coxet al.,
“Two-dimensional quantitative photoacoustic image reconstruction of absorption distributions in scattering media by use of a simple iterative method,”
Appl. Opt., 45
(8), 1866
–1875
(2006). http://dx.doi.org/10.1364/AO.45.001866 APOPAI 0003-6935 Google Scholar
L. Yinet al.,
“Tomographic imaging of absolute optical absorption coefficient in turbid media using combined photoacoustic and diffusing light measurements,”
Opt. Lett., 32
(17), 2556
–2558
(2007). http://dx.doi.org/10.1364/OL.32.002556 OPLEDP 0146-9592 Google Scholar
T. Jetzfellneret al.,
“Performance of iterative optoacoustic tomography with experimental data,”
Appl. Phys. Lett., 95
(1), 013703
(2009). http://dx.doi.org/10.1063/1.3167280 APPLAB 0003-6951 Google Scholar
P. ShaoB. CoxR. J. Zemp,
“Estimating optical absorption, scattering, and Grueneisen distributions with multiple-illumination photoacoustic tomography,”
Appl. Opt., 50
(19), 3145
–3154
(2011). http://dx.doi.org/10.1364/AO.50.003145 APOPAI 0003-6935 Google Scholar
T. HarrisonP. ShaoR. J. Zemp,
“A least-squares fixed-point iterative algorithm for multiple illumination photoacoustic tomography,”
Biomed. Opt. Express, 4
(10), 2224
–2230
(2013). http://dx.doi.org/10.1364/BOE.4.002224 BOEICL 2156-7085 Google Scholar
T. HarrisonA. SampaleanuR. Zemp,
“S-sequence spatially-encoded synthetic aperture ultrasound imaging,”
IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 61
(5), 886
–890
(2014). http://dx.doi.org/10.1109/TUFFC.2014.2979 ITUCER 0885-3010 Google Scholar
S. R. Arridge,
“Optical tomography in medical imaging,”
Inverse Probl., 15
(2), R41
(1999). http://dx.doi.org/10.1088/0266-5611/15/2/022 INPEEY 0266-5611 Google Scholar
N. J. A. Sloaneet al.,
“Codes for multiplex spectrometry,”
Appl. Opt., 8
(10), 2103
–2106
(1969). http://dx.doi.org/10.1364/AO.8.002103 APOPAI 0003-6935 Google Scholar
G. BalK. Ren,
“Multi-source quantitative photoacoustic tomography in a diffusive regime,”
Inverse Probl., 27
(7), 075003
(2011). http://dx.doi.org/10.1088/0266-5611/27/7/075003 INPEEY 0266-5611 Google Scholar
T. Correiaet al.,
“Selection of regularization parameter for optical topography,”
J. Biomed. Opt., 14
(3), 034044
(2009). http://dx.doi.org/10.1117/1.3156839 JBOPFO 1083-3668 Google Scholar
C. B. Shawet al.,
“Least squares QR-based decomposition provides an efficient way of computing optimal regularization parameter in photoacoustic tomography,”
J. Biomed. Opt., 18
(8), 080501
(2013). http://dx.doi.org/10.1117/1.JBO.18.8.080501 JBOPFO 1083-3668 Google Scholar
BiographyTyler Harrison is a PhD candidate under Dr. Roger Zemp at the University of Alberta. He received his BSc from the same institution in 2009. His work focuses on photoacoustic imaging: systems development, image reconstruction, and quantitative imaging. His other research interests include transducer fabrication (largely CMUTs) and ultrasound imaging. Peng Shao earned his PhD in electrical and computer engineering with a biomedical engineering specialization from the University of Alberta in 2014. His thesis work focused on quantitative reconstruction of optical properties, OR-PAM imaging, and photoacoustic microscopy. He has recently accepted a postdoctoral appointment at Harvard Medical School. Roger J. Zemp is an associate professor in the Department of Electrical and Computer Engineering at the University of Alberta. He received his BSc from the University of Alberta in 1998, his MSc from the University of Toronto in 2000, and his PhD in biomedical engineering from the University of California, Davis, in 2004. His research focuses on novel methods of biomedical imaging, including optical and ultrasound approaches. |