A low memory cost model based reconstruction algorithm exploiting translational symmetry for photoacustic microscopy.

A low memory cost model based reconstruction algorithm exploiting translational symmetry for photoacustic microscopy Juan Aguirre,1,* Alexia Giannoula,1 Taisuke Minagawa,1 Lutz Funk,2 Pau Turon,2 and Turgut Durduran1 1

ICFO-Institut de Ciènces Fotòniques, 08860 Castelldefels, Barcelona, Spain 2 B.Braun Surgical S.A., Rubí, Barcelona, Spain *[email protected]

Abstract: A model based reconstruction algorithm that exploits translational symmetries for photoacoustic microscopy to drastically reduce the memory cost is presented. The memory size needed to store the model matrix is independent of the number of acquisitions at different positions. This helps us to overcome one of the main limitations of previous algorithms. Furthermore, using the algebraic reconstruction technique and building the model matrix “on the fly”, we have obtained fast reconstructions of simulated and experimental data on both two- and threedimensional grids using a traditional dark field photoacoustic microscope and a standard personal computer. ©2013 Optical Society of America OCIS codes: (110.5120) Photoacoustic imaging; (100.3010) Image reconstruction techniques.

References and links 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15.

L. V. Wang and S. Hu, “Photoacoustic tomography: in vivo imaging from organelles to organs,” Science 335(6075), 1458–1462 (2012). L. V. Wang, “Multiscale photoacoustic microscopy and computed tomography,” Nat. Photonics 3(9), 503–509 (2009). V. Ntziachristos, “Going deeper than microscopy: the optical imaging frontier in biology,” Nat. Methods 7(8), 603–614 (2010). H. F. Zhang, K. Maslov, and L. V. Wang, “In vivo imaging of subcutaneous structures using functional photoacoustic microscopy,” Nat. Protoc. 2(4), 797–804 (2007). E. Z. Zhang, J. G. Laufer, R. B. Pedley, and P. C. Beard, “In vivo high-resolution 3D photoacoustic imaging of superficial vascular anatomy,” Phys. Med. Biol. 54(4), 1035–1046 (2009). V. Ntziachristos, “Clinical translation of optical and optoacoustic imaging,” Philos Trans A Math Phys. Eng. Sci. 369, 4666–4678 (2011). P. B. Garcia-Allende, J. Glatz, M. Koch, and V. Ntziachristos, “Enriching the Interventional Vision of Cancer with Fluorescence and Optoacoustic Imaging,” J. Nucl. Med. 54(5), 664–667 (2013). D. Razansky, A. Buehler, and V. Ntziachristos, “Volumetric real-time multispectral optoacoustic tomography of biomarkers,” Nat. Protoc. 6(8), 1121–1129 (2011). J. Yao, C. H. Huang, L. Wang, J. M. Yang, L. Gao, K. I. Maslov, J. Zou, and L. V. Wang, “Wide-field fastscanning photoacoustic microscopy based on a water-immersible MEMS scanning mirror,” J. Biomed. Opt. 17(8), 080505 (2012). R. Ma, S. Söntges, S. Shoham, V. Ntziachristos, and D. Razansky, “Fast scanning coaxial optoacoustic microscopy,” Biomed. Opt. Express 3(7), 1724–1731 (2012). L. Wang, K. Maslov, W. Xing, A. Garcia-Uribe, and L. V. Wang, “Video-rate functional photoacoustic microscopy at depths,” J. Biomed. Opt. 17(10), 106007 (2012). K. Maslov, G. Stoica, and L. V. Wang, “In vivo dark-field reflection-mode photoacoustic microscopy,” Opt. Lett. 30(6), 625–627 (2005). H. F. Zhang, K. Maslov, M. L. Li, G. Stoica, and L. V. Wang, “In vivo volumetric imaging of subcutaneous microvasculature by photoacoustic microscopy,” Opt. Express 14(20), 9317–9323 (2006). M. L. Li, H. E. Zhang, K. Maslov, G. Stoica, and L. V. Wang, “Improved in vivo photoacoustic microscopy based on a virtual-detector concept,” Opt. Lett. 31(4), 474–476 (2006). M. A. Araque Caballero, A. Rosenthal, J. Gateau, D. Razansky, and V. Ntziachristos, “Model-based optoacoustic imaging using focused detector scanning,” Opt. Lett. 37(19), 4080–4082 (2012).

#195566 - $15.00 USD

Received 13 Aug 2013; revised 27 Oct 2013; accepted 1 Nov 2013; published 12 Nov 2013

(C) 2013 OSA 1 December 2013 | Vol. 4, No. 12 | DOI:10.1364/BOE.4.002813 | BIOMEDICAL OPTICS EXPRESS 2813

16. X. L. Deán-Ben, A. Buehler, V. Ntziachristos, and D. Razansky, “Accurate model-based reconstruction algorithm for three-dimensional optoacoustic tomography,” IEEE Trans. Med. Imaging 31(10), 1922–1928 (2012). 17. A. Rosenthal, D. Razansky, and V. Ntziachristos, “Fast semi-analytical model-based acoustic inversion for quantitative optoacoustic tomography,” IEEE Trans. Med. Imaging 29(6), 1275–1285 (2010). 18. A. Rosenthal, V. Ntziachristos, and D. Razansky, “Model-based optoacoustic inversion with arbitrary-shape detectors,” Med. Phys. 38(7), 4285–4295 (2011). 19. S. Bu, Z. Liu, T. Shiina, K. Kondo, M. Yamakawa, K. Fukutani, Y. Someda, and Y. Asao, “Model-based reconstruction integrated with fluence compensation for photoacoustic tomography,” IEEE Trans. Biomed. Eng. 59(5), 1354–1363 (2012). 20. X. L. Dean-Ben, R. Ma, D. Razansky, and V. Ntziachristos, “Statistical approach for optoacoustic image reconstruction in the presence of strong acoustic heterogeneities,” IEEE Trans. Med. Imaging 30(2), 401–408 (2011). 21. M. Xu and L. V. Wang, “Universal back-projection algorithm for photoacoustic computed tomography,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 71(1), 016706 (2005). 22. P. Burgholzer, G. J. Matt, M. Haltmeier, and G. Paltauf, “Exact and approximative imaging methods for photoacoustic tomography using an arbitrary detection surface,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 75(4), 046706 (2007). 23. B. E. Treeby and B. T. Cox, “k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields,” J. Biomed. Opt. 15(2), 021314 (2010). 24. J. Laufer, F. Norris, J. Cleary, E. Zhang, B. Treeby, B. Cox, P. Johnson, P. Scambler, M. Lythgoe, and P. Beard, “In vivo photoacoustic imaging of mouse embryos,” J. Biomed. Opt. 17(6), 061220 (2012). 25. G. Paltauf, J. A. Viator, S. A. Prahl, and S. L. Jacques, “Iterative reconstruction algorithm for optoacoustic imaging,” J. Acoust. Soc. Am. 112(4), 1536–1544 (2002). 26. P. Ephrat, L. Keenliside, A. Seabrook, F. S. Prato, and J. J. Carson, “Three-dimensional photoacoustic imaging by sparse-array detection and iterative image reconstruction,” J. Biomed. Opt. 13(5), 054052 (2008). 27. K. Wang, S. A. Ermilov, R. Su, H. P. Brecht, A. A. Oraevsky, and M. A. Anastasio, “An imaging model incorporating ultrasonic transducer properties for three-dimensional optoacoustic tomography,” IEEE Trans. Med. Imaging 30(2), 203–214 (2011). 28. D. Queirós, X. L. Déan-Ben, A. Buehler, D. Razansky, A. Rosenthal, and V. Ntziachristos, “Modeling the shape of cylindrically focused transducers in three-dimensional optoacoustic tomography,” J. Biomed. Opt. 18(7), 076014 (2013). 29. L. H. Wang, Biomedical Optics: Principles and Imaging (Wiley, Hoboken, New Jersey, 2007). 30. C. G. Hoelen and F. F. de Mul, “A new theoretical approach to photoacoustic signal generation,” J. Acoust. Soc. Am. 106, 11 (1999). 31. K. P. Köstli and P. C. Beard, “Two-dimensional photoacoustic imaging by use of Fourier-transform image reconstruction and a detector with an anisotropic response,” Appl. Opt. 42(10), 1899–1908 (2003). 32. X. Intes, V. Ntziachristos, J. P. Culver, A. Yodh, and B. Chance, “Projection access order in algebraic reconstruction technique for diffuse optical tomography,” Phys. Med. Biol. 47(1), N1–N10 (2002). 33. G. T. Herman and L. B. Meyer, “Algebraic reconstruction techniques can be made computationally efficient [positron emission tomography application],” IEEE Trans. Med. Imaging 12(3), 600–609 (1993). 34. D. Ros, C. Falcón, I. Juvells, and J. Pavía, “The influence of a relaxation parameter on SPECT iterative reconstruction algorithms,” Phys. Med. Biol. 41(5), 925–937 (1996). 35. T. Shin, J. F. Nielsen, and K. S. Nayak, “Accelerating dynamic spiral MRI by algebraic reconstruction from undersampled k--t space,” IEEE Trans. Med. Imaging 26(7), 917–924 (2007). 36. A. H. Andersen, “Algebraic reconstruction in CT from limited views,” IEEE Trans. Med. Imaging 8(1), 50–55 (1989). 37. T. Durduran, R. Choe, W. B. Baker, and A. Yodh, “Diffuse optics for tissue monitoring and tomography,” Rep. Prog. Phys. 73(7), 076701 (2010). 38. M. B. Taisuke, A. Pau Giannoula, and T. Durduran, “MBio: a comprehensive Monte-carlo package for Diffuse Correñation Spectroscopy/Tomography,” in European Conferences on iomedical Optics, 2013) 39. A. Rosenthal, T. Jetzfellner, D. Razansky, and V. Ntziachristos, “Efficient framework for model-based tomographic image reconstruction using wavelet packets,” IEEE Trans. Med. Imaging 31(7), 1346–1357 (2012). 40. T. Koichi, W. Katsuhiro, F. Kazuhiko, and S. Tsuyoshi, “Advanced model-based reconstruction algorithm for practical three-dimensional photoacoustic imaginbg,” in SPIE, 2011) 41. B. Shuhui, L. Zhenbao, S. Tsuyoshi, and F. Kazuhiko, “Matrix compression and Compressed sensing reconstruction for photoacoustic tomography,” Elektonika Ir Elektrotechnika, 18 (2012). 42. B. E. Treeby, T. K. Varslot, E. Z. Zhang, J. G. Laufer, and P. C. Beard, “Automatic sound speed selection in photoacoustic image reconstruction using an autofocus approach,” J. Biomed. Opt. 16(9), 090501 (2011). 43. X. L. Deán-Ben, V. Ntziachristos, and D. Razansky, “Artefact reduction in optoacoustic tomographic imaging by estimating the distribution of acoustic scatterers,” J. Biomed. Opt. 17(11), 110504 (2012). 44. J. Chamorro-Servent, J. Aguirre, J. Ripoll, J. J. Vaquero, and M. Desco, “Feasibility of U-curve method to select the regularization parameter for fluorescence diffuse optical tomography in phantom and small animal studies,” Opt. Express 19(12), 11490–11506 (2011).

#195566 - $15.00 USD



45. J. L. Herraiz, S. España, J. J. Vaquero, M. Desco, and J. M. Udías, “FIRST: Fast Iterative Reconstruction Software for (PET) tomography,” Phys. Med. Biol. 51(18), 4547–4565 (2006).

1. Introduction Photoacoustic imaging (PAI) techniques [1–3] are non-invasive biomedical imaging modalities that seek to retrieve the three-dimensional (3D) distribution of light energy deposition due to absorption that may, in turn, correspond to the distribution of chromophores such as hemoglobin species in vivo. They are based on the photoacoustic effect which is produced when tissues are irradiated with a short light pulse. This pulse leads to a thermoelastic expansion that induces the emission of an acoustic wave that can be detected by an appropriate detector. By employing appropriate data processing schemes, the distribution of the absorbed light energy can be obtained with ultrasonic- and/or optical scale spatial resolution even in deep, i.e. >>1 mm, tissues. Since light absorption serves as the source of contrast PAI provides functional and metabolic images of the vasculature, hemodynamics and oxygen metabolism [2, 4]. Many applications of PAI in preclinical research and clinical diagnostics have been envisioned and attempted [1, 5–8]. One class of PAI systems utilizes raster scanning or parallel detection at a large number of coplanar source points [9–11]. This paper is mainly related to this class of PAI systems. In particular, for example, the dark-field illumination photoacoustic microscopy (dark-field PAM) [12] has allowed for unprecedented in vivo imaging of the microvasculature as deep as 2-3 mm below the tissue surface. Dark-field PAM and most other raster scan systems utilize a spherically focused transducer directed toward the subject to be imaged. The tissue is illuminated by a ring of light and the ultrasound detector, often a transducer, is placed at the center of the field, hence the term “dark-field,” avoiding the strong photoacoustic signal that is generated due to absorption by the melanin on the skin surface. The transducer together with the illumination optics is placed on a scanning head that performs a raster scan on the x-y plane. At each position several “amplitude mode” lines (“A lines”) are obtained in the depth direction, i.e. z direction, that allows us to put together several two-dimensional (2D) planes or “brightness mode” lines (“B-lines”) yielding a 3D image. Relatively high lateral (~45 µm, x-y plane) and axial/depth resolution (~15 µm, z- direction) [12] is achieved by using transducers with a high numerical aperture, a high central frequency and a large bandwidth. There is, however, a tradeoff between the lateral resolution and the depth of focus (0.3 mm for the ultrasonic parameters commonly used [12]), hence, the spatial resolution of dark-field PAM images and other raster scan systems that are based on spherically focused transducers rapidly degrade above and below the focal zone of the ultrasound detector. This leads to systematic image artifacts and the lateral resolution is particularly affected. To minimize this limitation, contour scanning [13] and synthetic aperture focusing technique (SAFT) [14] have been proposed. In contour scanning, first, a preliminary image is obtained by performing a standard raster scan. Then, from this image, the skin contour of the subject under study is obtained which allows a second scan while dynamically adjusting the height (z coordinate) of the scanning head in order to keep the ultrasonic focus at a constant distance from the skin. By doing this, several improvements, such as, the uniform focusing of the subcutaneous blood vessels were reported [13]. However, the acquisition time is doubled and the depth of focus is not actually improved. The SAFT approach is a signal processing algorithm. When it is applied to spherically focused transducers, it assumes that a virtual detector is placed in its focal point. The acquired A-lines are delayed and summed in accordance with the position of this virtual detector at an attempt to correct for the realistic detector geometries. This algorithm improves the resolution of the out of focus regions of the PAM images while improving the signal-to-noise ratio. However, it has recently been shown [15] that it fails to account for some of the image artifacts that are due to the finite shape of the spherical transducers

#195566 - $15.00 USD



The next step was to introduce model based image reconstruction techniques [15, 16]. These techniques utilize a forward model for the acoustic propagation which leads to a linear system of equations whose factors form the so-called model matrix. The system can then be solved (inverse problem) leading to photoacoustic images. The main advantage of model based photoacoustic reconstructions is that several effects can be taken into account in the forward problem, such as the frequency response of the transducer [17], the spatial distribution of the transducer [18] and the light fluence corrections [19]. Furthermore, the algorithm can be modified to reduce artifacts due to acoustic scatterers [20] and/or take into account a priori information regarding the photoacoustic reconstruction such as nonnegativity. It has been suggested, for a general photoacoustic raster scanning configuration, that the results obtained using model based image reconstruction techniques clearly improve the results obtained by SAFT [15]. Other strategies based on closed-form back projection formulae [21, 22] do not, in general, allow the researchers to explicitly take the geometry of focused detectors or nonuniform illumination into account. Another option is the time-reversal algorithms which assume that the photoacoustic signals eventually decay to zero inside the reconstruction domain (Huygens principle) [22, 23]. These algorithms have been applied successfully in 3D [24]. By assuming this, the differential equation (wave equation) that describes the ultrasound propagation can be solved backwards in time using a zero initial condition. This class of algorithms lead to accurate reconstructions; however, they might not be computationally feasible for reconstructing the high frequency components of photoacoustic signals. Here we focus on the model based reconstructions. A common strategy for building the forward model consists on assuming that the light propagation is a delta function in time which allows us to express the photoacoustic pressure at a certain position and a certain time point as a poisson type integral defined over a spherical surface. By discretizing the integral and interpolating it in a cubic grid a system of linear equations is derived [16]. Other approach consists on deriving the model matrix by assuming that the measured photoacoustic signals are formed by the sum of the incoming signals generated by homogenous spherical emitters situated at center of each voxels of the reconstruction grid [25–27]. Unfortunately, for most problems (see Section 2.2 and 2.3), this requires very large matrices to be stored and manipulated making the problem intractable in standard computer platforms [28]. In this paper, we report a low memory cost time domain model based reconstruction algorithm for planar data acquisition systems and test the algorithm using simulations and experiments on phantoms using a dark field PAM. This approach and the use of the algebraic reconstruction technique (ART) for inversion allow us to carry out the image reconstruction on a standard personal computer at reasonable time-scales. At times, this approaches real-time imaging even without any specialized optimization or parallelization. The paper is organized as follows: In Section 2.1 we derive the theoretical model matrix from the forward model, in Section 2.2 we describe how the transducer shape and light fluence are taken into account in the model. In Section 2.3 we describe how the system symmetries allow us to reduce the matrix size. The implementation of the ART algorithm is described in Section 2.4. Sections 2.5 and 2.6 are devoted to describe the simulations and the phantom experiment. The results are shown in Section 3 and are discussed in Section 4. 2. Materials and methods 2.1 Forward modeling We build a forward model by assuming that the probed medium is acoustically homogeneous. A light pulse (whose profile is approximated by a Dirac delta in time) is shone on this   medium creating a photoacoustic pressure wave, p ( r , t ) , at position r and time t . Assuming that the thermal confinement condition is fulfilled, the production and propagation of  p ( r , t ) is described by the following expression [16]:

#195566 - $15.00 USD



 ∂ 2 p( r , t ) 2 2   ∂δ (t ) − cs ∇ p( r , t ) = ΓH ( r ) . 2 ∂t ∂t

(1)

Here Γ is the dimensionless Gruneisen parameter, cs is the speed of the sound in the medium,  and λ is the heating function. H ( r ) is defined as the thermal energy absorbed by the tissue per  unit volume. μa ( r ) is the optical absorption coefficient at the wavelength of the pulse and  φ ( r ) is the light fluence.  When a sphere located at rj of radius R is homogenously heated up at t = 0 , an initial    pressure p0 ( rj ) = Γμa ( rj )φ ( rj ) is generated inside the sphere. The resulting photoacoustic  pressure wave at position r situated outside the sphere can be analytically calculated from Eq. (1), as [29, 30]   r − ri − cs t            p ( r , rj , t ) = p0 ( rj ) U ( R − c s t − r − rj ) +   U ( r − rj − R − c s t )U ( R + c s t − r − rj )  2 r − r   j

Here U is the Heaviside step function. Equation (2) can be rewritten as      p( r , rj , t ) = p0 ( rj ) psp ( r , rj , t ),

(2)

(3)

  where psp ( r , rj , t ) represents the photoacoustic pressure generated by the sphere with an initial

pressure p0 = 1 N/m2 .  Our forward model assumes that the pressure p( r , t ) measured by a detector is a linear combination of the pressure generated by a set of spheres, each of them homogeneously  heated, located at the N points, {rm } ( m = 1.... N ) , of a regular, Cartesian grid. Using Eq. (3) the model can then be expressed as N     p( r , t ) ≈  po ( rj ) psp ( r , rj , t )

(4)

j =1

If the medium is homogeneously illuminated and several point like detectors are placed in a  set of locations, {rd } ( d = 1....S ) that record the pressure wave at several time points, {tk } ( k = 1....T ) then from Eq. (4) the following system of linear equations can be obtained: N d = 1...S     p( rd , tk ) ≈  po ( rj ) psp ( rd , rj , tk ), where k = 1...T j =1

(5)

  The known coefficients are psp ( rd , rj , tk ) and can be analytically calculated from Eq. (3). The

unknowns are the initial photoacoustic pressure at each voxel of the reconstruction grid, i.e.    p0 ( rj ) = Γμa ( rj )φ ( rj ) . In matrix notation, the set of linear equations described by Eq. (5) can be expressed as b = W  x

(6)

  where W is the model matrix [16] and its elements are w ij = psp ( rd , rj , tk ) where

i = ( d − 1)T + k .The first T rows of the matrix correspond to the detector d = 1, the next T rows to the detector d = 2 and so on. x is a vector that contains the unknowns

#195566 - $15.00 USD



   p0 ( rj ) = Γμa ( rj )φ ( rj ) , and, b is a vector that contains the measured photoacoustic signal for

every detector at every time point. 2.2 The modeling of the transducer, and the light fluence correction

According to [31] the response to an acoustic wave of a finite detector defined by a surface A  and a position rd is given by     pdet ( rd , t ) =  p( r , t )D ( r )dr

(7)

  1∀r ∈ A . where D ( r ) =   0∀r ∉ A Equation (7) does not take into account the convolution of the theoretical acoustic wave form with the impulse response function of the detector. This expression would be correct if the acquired signals are deconvolved with the impulse response function. Let us discretize the  3D space in an infinite set of homogeneous voxels of volume h3. We can then define {rd ,l }  (l = 1....B ) as the positions of the voxels for which D ( r ) = 1 corresponding to the dth detector. By discretizing the right hand side of Eq. (7) and using Eq. (4), we obtain N N   B      pdet ( rd , tk ) ≈ h 3  po ( rj ) psp ( rd ,l , rj , tk ) =  po ( rj ) w( rd , rj , tk ) j =1

l =1

(8)

j =1

Equation (8) is another linear system of equations that can be written as b = [W ] x.

(9)

B

  Here wij = h 3  psp ( rd ,l , rj , tk ) where again i = ( d − 1)T + k . l =1

Unfortunately, calculating W from expression 8 is computationally expensive. Some authors approximate Eq. (7) to a spatial impulse response function which is computationally more manageable [31]. Temporal impulse functions based on assuming point like sources have been also proposed in the past [18]. To give an example, for the phantom experiment described in Section 2.6, we have computed Eq. (8) directly by launching ten processes in parallel on an Intel Xeon (R) CPU ES-1660 @3.30 GHz x 12 work station. In our case it had to be calculated only once and for only one detector position, as discussed in Section 2.3, the computation time was of approximately one week. In a more traditional approach, it would scale roughly linearly with the number of measurement positions making it an unmanageable problem. To include the light fluence in the forward model, Eq. (8) has to be reorganized as follows N N         pdet ( rd , tk ) ≈  Γμ ( rj )φ ( rj )w( rd , rj , tk ) =  pˆ o ( rj )wˆ ( rd , rj , tk ). j =1

(10)

j =1

Equation (10) represents a new system of equations in which the unknowns are   pˆ 0 ( rj ) = Γμ ( rj ) and the coefficients are weighted by the fluence      as wˆ ( rd , rj , tk ) = φ ( rj ) w( rd , rj , tk ) . 2.3 Reduction of the matrix size by exploiting the translational symmetry

We imagine a PAM experiment where the data is acquired on the x-y plane covering a rectangular region of size Lx x Ly. The spacing of the scanning points is constant in the x- and y- directions and the homogeneous cubic grid used for reconstruction convers the same region #195566 - $15.00 USD



on the x-y plane as the scanning points. Let us assume that the number of the reconstruction grid points in the x-y plane is proportional to the number of scanned points and that the spacing of the grid points is equal or finer than the spacing of the scanning points. Therefore  for each detector position defined by the coordinates of its focal point, rd = ( xd , yd , zd ) there  is a subset of voxels with coordinates rvox = ( xvox , yvox , zvox ) that fulfill the conditions xvox = xd and yvox = yd . In order to model the light fluence, we further assume that the background medium has homogeneous optical properties. This assumption allows us to carry out a first order correction whose utility was previously demonstrated [19] albeit with some systematic errors in quantification.This is a reasonable compromise between complexity and speed. If we also assume that the surface of the imaged object is parallel to the x-y plane, and, that the time discretization of the photoacoustic signal is equal for all the detector positions, the following simple relationship can be fulfilled due to the translational symmetries of the system:       (11) w( rd 2 , rj , tk ) = w( rd 1 , rj − ( rd 1 − rd 2 ), tk ).     Equation (11) is true as long as rj − ( rd 1 − rd 2 ) ∈ {rm } (see Fig. 1).

Fig. 1. : (Top) Representation of a row of the weight matrix W corresponding to a specific time point and a detector characterized by the position of its focal point rd1. (Bottom) Representation of a row of the weight matrix for the same time point, corresponding to the same detector but now situated at position rd2 (focal point). A large area of the bottom image can be seen as a translation of the top image. The weights correspond to a circular transducer with a radius of 12.7 mm and a numerical aperture of 0.36 that is, situated in the positive part of the z axis.

 An alternative to Eq. (11) is to calculate a generator matrix, W , that would correspond to the model matrix of an experiment with a single scanning point that corresponds to the central scanning point of the experiment corresponding to the model matrix W. In this new imaginary

#195566 - $15.00 USD



experiment, the reconstruction grid density would be the same, however, it would cover a spatial region of size 2Lx x 2Ly. In this case, the following equation is always fulfilled.        w( rd 2 , rj , tk ) = w( rd 1 , rj − ( rd 1 − rd 2 ), tk ). (12)   where w are the elements of the generator matrix W . The size of W is equal to N columns and T x S rows where N is the number of grid points, T is the number of time points from the detected photoacoustic signal that are used for the  reconstruction and S is the number of the detection points, whereas the size of W is equal to 4N columns and T rows, being independent of the number of detection points.

2.4 The inverse problem

To solve the inverse problem we have used the ART reconstruction algorithm which can be expressed as [32] xi( l +1) = xi( l ) + λ

bi −  j wij x (jl )



j

wij2

wi .

(13)

Where xi( l ) is the lth estimate of the ith row contribution to the unknown reconstructed image. wi is the ith row vector of the weight matrix and bi the ith measurement. ART works by projecting and initial guess of the reconstructed images onto one hyperplane defined by a row of the system of equations to be solved. This process is done iteratively through all the hyperplanes defined by all the rows of the system of equations. Once the algorithm has gone through all the rows, an iteration of the algorithm is completed. λ is the relaxation parameter which adjusts the projection step for each iteration. The value of λ is commonly calculated empirically [33, 34] and it determines the number of iterations needed to obtain a satisfactory reconstruction. It plays a similar role as the regularization parameter when performing the commonly employed Tikonov regularization for reconstruction purposes [32]. In general, care must be taken when choosing λ and the number of iterations. Our results showed that using a relaxation parameter of 0.01 and 3 iterations gave the most suitable reconstructions for the purposes of this study. The solution of the inverse problem was programmed in C, using Eq. (13), and accessing each element of each row of the matrix on the fly from the pre-calculated generator matrix using Eq. (12). Assuming that the full matrix is arranged as described in Section 2.2, i.e. made of S different sub-matrices, each of them corresponding to a detector d, Eq. (12) has to be applied only S times for every iteration of the ART algorithm. This adds only a minor delay in its performance. The ART algorithm has been widely used for different imaging modalities, like diffuse optical tomography, single photon emission computed tomography, X-ray CT or MRI [33– 36]. It is well suited for photoacoustic tomography since it easily allows applying positive constrain for the imaged object. 2.5 Numerical experiment

2.5.1 Data generation To show the performance of the algorithm, a single B-scan was simulated and reconstructed on a 2D grid (x-z plane). The simulation consisted of five homogeneously heated spheres situated along the same line in the z-direction, i.e. at a different depth, each of them with a diameter of 300 μm with their centers separated by 500 μm. The geometry is illustrated in Fig. 2. For simplicity a circular transducer was defined with a radius of 12.7 mm and a numerical aperture of 0.36. The z-coordinates of the two spheres was situated below z-coordinate of the

#195566 - $15.00 USD



focal point of the transducer, two others were above, and the z- coordinate of the center of the remaining one was coincident with the z coordinate of the acoustic focal point. The B-scan was simulated in the x direction covering 5 mm with a step size of 50 μm. The signal corresponding to each sphere was calculated using Eq. (2) and Eq. (7) by assuming p0 = 1 N/m2. The discretization of Eq. (7) was done by taking 300 points in the detector line. 2.5.2 Reconstruction The reconstruction is carried out on a grid covering 3 mm in the z-direction and 5 mm in the x-direction with a homogeneous spacing of 50 μm in both directions. The generator matrix was calculated assuming that each voxel had a sphere of radius also of 50 μm, using Eq. (8). The detector was discretized as was done in Section 2.5.1. The time resolution for sampling the photoacoustic signal was taken roughly as two times finer than the spatial resolution of the grid ( Δt ⋅ cs = Δx / 2 ) i.e. Δt = 0.017 μs [17].

Fig. 2. A schematic of the simulated objects for the numerical simulation. Five spheres of 300 μm diameter separated by 500 μm each are shown. In the simulation each sphere is homogeneously heated.

2.6 The phantom experiment

2.6.1 The ICFO-PAM system The imaging system we have used is based on the one presented in Ref [12]. Briefly, the laser that is used to generate the photoacoustic excitation is a tunable near-infrared Ti:sapphire laser (LT-2211A, LOTIS TII, Minsk, Belarus) pumped by a Q-switched Nd:YAG laser (LS2137/2, LOTIS TII, Minsk, Belarus). The laser pulse width is ~10 ns and, for this experiment, we have tuned it at 700 nm with a repetition frequency of 10 Hz. The light is delivered through an optical fiber to the PAM scanner head. Dark-field confocal (optical and transducer) illumination is achieved by allowing the laser beam to pass from the fiber through an axicon towards a concave mirror (optical condenser) which “focuses” (it is focused in case there is no turbid medium present) the beam into the tissue with the focal region coaxially overlapping the ultrasonic focus inside the tissue. The system is illustrated in Fig. 3. The optical illumination on the surface of the excited sample is ring shaped with a dark center.

#195566 - $15.00 USD



Fig. 3. The diagram of the ICFO-PAM scanning head.

To detect the generated photoacoustic waves, a high-frequency (20 MHz) spherically focused (f = 13mm) ultrasonic transducer (V214-BB-RM, Olympus, Essex, UK) was used. The electrical signals obtained from the transducer are then amplified, digitized and stored by an oscilloscope with an embedded personal computer (Wavesurfer 24-MXsA, Lecroy, Geneva, Switzerland). The oscilloscope also controls an XY-linear translation stage (NRT100/M, Thorlabs, Dachau, Germany) allowing the PAM head (lens, mirror and transducer) to perform a two-dimensional (2-D) scan. For our experiments, the scanning head is introduced in a water container with a window sealed with an acoustically and optically transparent plastic membrane (Fig. 3) to allow imaging liquid phantoms. 2.6.2 Phantom, acquisition and reconstruction parameters We have prepared a phantom where a strongly absorbing thread (poly 4-hydroxybutyrate, MonoMax®, B. Braun Surgical, S.A., Barcelona, Spain) of ~300 µm diameter was immersed ~3mm deep from the surface in a solution of distilled water and Lipofundina (B. Braun Medical, S,A,, Barcelona, Spain). The optical properties at 700 nm were µ’s ≈10 cm−1 and µa ≈0.0057 cm−1. The optical properties of the phantom were confirmed by time resolved spectroscopy [37]. The thread was held in place tightly laid parallel to the liquid surface. The raster scan covered a 0.4x0.4 cm2 region in the x-y plane with a scanning step of 50 µm in both directions. For the reconstruction, a grid of 160 x 160 x 81 (x,y,z) elements was selected using a voxel separation of 25 µm. To calculate the generator matrix, from Eq. (8), the surface of the transducer at each position was discretized into 7669 volume elements using an extension of the mid-point circle algorithm [23]. The time sampling was of the photoacoustic signal was of Δt = 0.017 μs. 2.6.3 The simulation of the light fluence Dark field illumination produces a ring shaped image. From the geometry and properties of the optics involved in the illumination, the angle of incidence of the light was obtained. The inner and outer diameter of the ring was measured by adjusting the optical fiber to a continuous source of white light, and measuring them in situ with a caliber. We ran a MonteCarlo based photon propagation code which employs the method which is described in [29] but has the capability of defining an inhomogeneous medium and a user-defined source with an arbitrary shape [38]. The simulation was done in a grid of 100x100x50 voxels, the size of each voxel was 0.1 mm in every direction. The calculated fluence was interpolated onto a grid of the same size of the one used for calculating the matrix (Eq. (10). #195566 - $15.00 USD



3. Results

The ART algorithm, as described in Section 2.4, was used to reconstruct the numerically simulated data on a 2D grid. The reconstruction was performed using an Intel Core i5-2400 CPU operating at 3.10 GHz with 4 GB of RAM. The reconstruction time was 0.9 seconds and the memory needed to store the generator matrix was 5.49 Mb which is roughly 50 times less than the memory that would have been needed to store the equivalent full model matrix. This data is shown in Table 1. Table 1. Reconstruction time and size of the matrix corresponding to the numerical experiment Reconstruction grid 2D (simulated experiment) 3D (phantom experiment)

Reconstruction time 0.9 s 6hrs

Size of the generator matrix 5.49 Mb 3.65 Gb

Size of the equivalent full matrix 277.3 Mb 5980 Gb

Figure 4 shows the simulated acquired PAM raw data together with its corresponding reconstructed image. The five spheres are resolved and the quality of the image remains constant along the z direction. This indicates that the model based algorithm has allowed us to correctly model the effect of signals below, at and above the ultrasound transducer’s focal depth. However, all the reconstructed objects are slightly smeared (compare to Fig. 1) in the x-direction leading to an ellipsoidal shape. There are also systematic image artifacts due to the fact that the detector does not rotate 360 ° around the imaged objects, and, therefore, the information content of the data set is not complete. A larger scanned area and/or an optimized regularization scheme would improve the quality of the image.

Fig. 4. (Left) Raw, simulated PAM data of the imaged objects obtained after the linear scan. (Right) Model based reconstruction after three iterations of the ART algorithm using the low memory cost generator matrix. All five objects (shown in Fig. 1) are clearly resolved in both dimensions. All data are normalized to the maximum.

The reconstruction of the phantom experiment on a 3D grid took six hours of computation using the same Intel Core i5-2400 CPU operating at 3.10 GHz with a 4 GB of RAM. The memory needed to store the generator matrix was of 3.65 GB which is well within the range of the current RAM of standard desktop PCs. The equivalent full matrix would have been roughly 1600 times bigger. The number of scanning points is 80x80 (0.4 cm x 0.4 cm, 50 µm step size), therefore expected benefit is 80x80/(2x2) = 1600 as illustrated in Table 1. In Fig. 5 we show the reconstructed 3D image from the experimental data. The thread is clearly resolved appearing on the bottom right hand side of the reconstructed image. More advanced modeling including the frequency response function of the transducer and an optimized regularization scheme would improve the reconstruction especially in the zdirection. Other effects related to the acoustic impedance of the thread, the light absorption profile in the thread itself and the limited angle of the acquired data further reduce the image quality.

#195566 - $15.00 USD



Fig. 5. The 3D reconstruction of the thread after three iterations of the ART algorithm using the low memory cost generator matrix. The reconstructed image is represented with 30 equally spaced axial slices. The image is normalized to the maximum.

The effects of the simple model appear more obvious if a single slice in the z-y plane or the z-x plane is observed (Figs. 6 and 7). The thread is reconstructed with an elliptical shape, and some artifacts arise due to reflections. .

Fig. 6. Representation of a slice in the z-y plane of the 3D reconstruction corresponding to x = 2.9 mm.

Fig. 7. Representation of a slice in the z-x plane of the 3D reconstruction corresponding to y = 4.5 mm

4. Discussion and conclusions

In this paper, we have presented a model based reconstruction algorithm which has been applied for planar photoacoustic microscopy. Its main contribution is that it successfully

#195566 - $15.00 USD



exploits the translational symmetry of the data acquisition system in order to greatly reduce the memory size needed to store the model matrix. For this approach, the model matrix is independent of the number of acquisition points as long as the conditions described in Section 2.3 are fulfilled, and, is generated by a single, generator matrix that can be managed by current standard personal computers. By using the algebraic reconstruction technique, fast reconstructions from single B-scans in a 2D plane were obtained in a numerical simulation covering a plane of 5 x 3 mm2 with a grid spacing of 50 µm and 100 scan points. Finally a reconstruction on a 3D grid from experimental data, obtained in acceptable practical time using a standard personal computer, is presented. In general, the size of the model matrix is equal to N columns and T x S rows, where N is the number of grid points, T is the number of time points from the detected photoacoustic signal, that are used for the reconstruction and S is the number of the detection points. Typically, a dark-field PAM raster scan covers areas ranging from 0.5 x0.5 cm2 to 1x1 cm2 and the scan step is of 100 μm (this quantity depends on the lateral resolution of the transducer). If we assume that the reconstruction grid spacing in the x,y and z direction is of 50 μm that the time resolution for sampling the acoustic field measurements is two times finer than the spatial resolution of the grid and that the reconstruction field of view (FOV) should cover 3 mm in the Z direction, this would lead to a matrix of the order of 1011 elements. For this reason, even though the rows of the matrix are sparse, its whole size is too large to be handled by a standard PC. To overcome this problem, the reconstructions are often done separately for each B-Scan at the expense of neglecting the out of plane signals, e.g. see Ref [15]. In this work, we have been able to significantly reduce the size of the model matrix allowing it to be handled by the random-access memory of current standard PC’s even in the case of full 3D reconstructions. Furthermore, the generator matrix contains all the information needed to build the whole matrix for the acquired data with different number and positions of scanning points, hence it has to be calculated once and can be re-used for different experiments. 4.1 Comparison to previous work

Other authors, e.g. Rosenthal et al [39], have also proposed an efficient method to reduce the size of the matrix for general model based tomographic reconstruction. Under this scheme the imaged object and the acquired data are represented in a wavelet packed base, leading approximately to a block-diagonal matrix where each block corresponds to a single spatialfrequency band in the imaged object. Each frequency is reconstructed separately, leading to a significant reduction of the system matrix and fast reconstructions. However, first order corrections are needed to account for effect that arise from neglecting the coupling between different spatial-frequency bands and, even in this case, the memory needed to reconstruct on a 3D grid is too high. On the other hand, in reference [25], a strategy based on the compression of the model matrix is presented, allowing to reduce the matrix size by a factor of two. However, this reduction is not enough to be handled in a standard personal computer for the case of true 3D reconstructions. Tanji et al [40] have previously suggested the utilization of translational symmetries and a generator matrix in the context of planar imaging and demonstrated its utility on simulations. However, no experimental data was presented and it may not be suitable for raster scanning. Following this work, other authors have used the symmetries of the sound wave propagation to propose matrix compression methods that can achieve a reduction of the matrix size of a factor of 1/250 [41]. However, with the method presented in this article, the size of the generator matrix becomes independent of the number of detection positions. We have built the model matrix assuming that every point of the reconstruction grid is equivalent to a uniform spherical source using a simple and intuitive derivation, a more rigorous derivation can be found in [27]. The low memory cost strategy could be equally

#195566 - $15.00 USD



implemented if the matrix would be calculated from interpolation of the surface where the Poisson like integral is defined as described in the introduction section. Other authors [25] have proposed to utilize a forward problem assuming that the pressure profile measured at each detector could be described as a linear combination of individual photoacoustic pressure profiles generated by spherical sources located on a 3-D grid of voxels. Particularly, in [25] the photoacustic pressure is integrated to recover the velocity potential from the heating function. 4.2 Limitations of the forward modeling and the quality of the reconstructions

In order for Eq. (2) to hold, the stress confinement condition must be fulfilled [30]. This condition states that the stress relaxation time, defined as τ = d/cs, where d is the characteristic dimension of the heated region must be much shorter than the laser pulse width. Under this assumption, stress propagation is negligible during the pulse. If this condition is not fulfilled Eq. (2) must be convolved with the temporal shape of the heating function. If we set d to be the diameter of the spheres, we obtained for the experimental case that τ = 16 ns which is of the same order of the laser duration (10 ns). All these specific modeling issues are beyond the objectives of this paper. The symmetries of the system are preserved only if homogeneous acoustic properties are assumed. Regarding the homogeneity of the speed of sound, high quality reconstructed images can be obtained by choosing an appropriate average speed of sound, which can be determined with the aid of focusing algorithms [42]. The issue of the presence of acoustic scatterers can be dealt by discriminating the scattered signals from the acquired data sets using probalisitic approaches [43], whereas further work should focus on the homogeneity of the acoustic attenuation and the Grunesesin coefficient for the regions assessed in this work (3mm depth subcutaneous imaging). 4.3 Reconstruction times

The ART algorithm has been widely used for diffuse optical tomography problems [32]. In this context, ART is not cataloged as a fast algorithm, since many iterations are needed to reach a meaningful reconstructed image. This is due to the strong ill-posed nature of this kind of imaging problems [32]. The photoacoustic reconstruction problem, when distributing the detectors covering 360° around the imaged object, and using a sufficiently accurate forward problem is a well posed problem [17]. However, the photoacoustic co-planar raster scan problem lacks in the number of geometrical projections, i.e. the data is acquired at a limited number of angles and some assumptions are taken in the forward modeling which makes it a relatively low ill posed problem. As an example, the condition number, defined as the ratio between the highest and the lowest eigen-value for the matrix of the 2D reconstruction described in Section 2.5 is 52.7, whereas the condition number of a typical fluorescence diffuse optical tomographic problem is ~106 [44]. In ART, the algorithm starts by projecting an initial estimate of the solution onto the hyperplane defined by one row of the linear system. This projection is used as the next estimation of the solution and is projected onto the hyperplane defined by another row of the linear system. This process is done iteratively until a satisfactory reconstructed image is obtained. Since at each step the algorithm only uses one row of the matrix, the accessing order to the elements of the generator matrix has to be changed only S times every time the algorithm goes through the model matrix (see section 2.4) introducing a minor delay in the computational time of the algorithm. For the single B-scan case presented in this work, the equivalent matrix would have had 10201 x 6161 elements, as opposed to the generator matrix for our algorithm that had only 101x2x6161 elements. This, in turn, implied a reconstruction in 0.9 s which is significantly better than the usual speeds achieved by common model based algorithms which generally require many minutes or hours to perform a similar task [17]. This is not too far from the real

#195566 - $15.00 USD



time imaging (albeit without a model and quantification). We note that this can readily be improved by parallelization using simultaneous algebraic reconstruction algorithm 4.4 Other symmetries and applications to other PA systems

The spherical symmetry of the weights around the axis of the transducer can be also taken into account leading to a higher reduction of the model matrix. Furthermore, compression techniques [25] can further reduce the size of the matrix, however, the reduction in the speed of the reconstruction would be affected by this further symmetry and the compression. We note that this reduction was critical for the 3D experimental case. We note that other reconstruction strategies that exploit the symmetry of the acquisition system similarly to the one shown in this paper have been presented in the past for other tomographic techniques like PET [45] that use rotating gantrys for data acquisition. Following similar reasoning, this strategy could be implemented for the wide range of photoacoustic systems that use acquisition set ups based on rotational or cylindrical geometries. This kind of systems can greatly benefit by current model based reconstruction algorithms since they exhibit superior modeling capabilities. However, current implementations in a standard personal computer are restricted to reconstructions on 2D planes, neglecting the full modeling capabilities of this kind of algorithms [28]. Acknowledgments

The authors would like to thank B. Braun Surgical, S.A., Fundació Cellex Barcelona, and the Spanish Ministerio de Economía y Competitividad for their funding support. We acknowledge useful discussions with Dr Romain Quidant, Dr Jan Laufer and Dr Konstantin Maslov.

#195566 - $15.00 USD



The Index-based Subgraph Matching Algorithm with General Symmetries (ISMAGS): exploiting symmetry for faster subgraph enumeration.

Beyond biologics: absorbable mesh as a low-cost, low-complication sling for implant-based breast reconstruction.

Ultra-low-dose CT for left atrium and pulmonary veins imaging using new model-based iterative reconstruction algorithm.

Short-term memory for symmetry.

Exploiting parameter sparsity in model-based reconstruction to accelerate proton density and T(2) mapping.

Low-cost model for teaching ophthalmoscopy.

Image reconstruction for structured-illumination microscopy with low signal level.

Microscopy illumination engineering using a low-cost liquid crystal display.

A validated low-cost training model for suprapubic catheter insertion.

A Limited-Memory BFGS Algorithm Based on a Trust-Region Quadratic Model for Large-Scale Nonlinear Equations.

2 Regularization.

A Novel Iterative CT Reconstruction Approach Based on FBP Algorithm.

Improving pulmonary vessel image quality with a full model-based iterative reconstruction algorithm in 80kVp low-dose chest CT for pediatric patients aged 0-6 years.

Explicit Filtering Based Low-Dose Differential Phase Reconstruction Algorithm with the Grating Interferometry.

A hippocampal indexing model of memory retrieval based on state trajectory reconstruction.

An improved phase shift reconstruction algorithm of fringe scanning technique for X-ray microscopy.

A blueprint for cost-efficient localization microscopy.

Implementation of the Hungarian algorithm to account for ligand symmetry and similarity in structure-based design.

Exploiting protein symmetry to design light-controllable enzyme inhibitors.

Memory and Perception-based Facial Image Reconstruction.

Template Deformation-Based 3-D Reconstruction of Full Human Body Scans From Low-Cost Depth Cameras.

A scalable and memory-efficient algorithm for de novo transcriptome assembly of non-model organisms.

Online learning algorithm for time series forecasting suitable for low cost wireless sensor networks nodes.

Focal-series reconstruction in low-energy electron microscopy.