Selected mechanical parameters of human brain and skull
1. Introduction
Thermoacoustic tomography (TAT) [1][3] is a novel, noninvasive medical imaging technique that detects the large differences in microwave absorption between pathological and normal tissue [4], [5].
It applies the principle of the thermoacoustic effect [6], an effect that arises from the combination of the pressure oscillations of a sound wave with the accompanying adiabatic temperature oscillations [7]. Within TAT an input microwave impulse stimulates thermoexpansion in biotissue and consequently generates acoustic waves (P wave) to be recorded by transducers arranged outside the tissue. When the tissue is relatively uniform, the initial local acoustic amplitude is approximately proportional to the absorption ratio of the microwave [1], [8], [9]. Consequently, the TAT imaging problem is to retrieve the distribution of the initial acoustic amplitude based on recorded acoustic wave energy, and thus it can be characterized as a problem of “multiple sources localization”.
In recent years, imaging by using thermoacoustic effect has raise increasing concerns [10], [11]. In which TAT has also been well studied for brain imaging [8], [9], [12]. The reasons include the following: first and far more important, TAT takes advantage of deep penetration of the electromagnetic impulse and the high resolution of the ultrasonic wave for deep imaging. Second, brain tissue is fundamentally uniform and isotropic. For example, the acoustic velocity of human brain is narrowly ranged within 14831521 m/s [13]. Acoustic wave propagation inside the brain is very close to oneway lineofsight transmission without much aberration. The most commonly used imaging algorithm of TAT is to backproject the wave energy recorded by each transducer along the ray paths to all possible locations inside the imaging domain [1], [9]. The backprojection method [15][18] has been reported to be well performed on its simple scheme with high costeffectiveness [19]. However, the approximation of oneway sight transmission used in the backprojection algorithm is no longer valid when high velocity contrast exists, for example, when the skull with an acoustic speed of 25002900 m/s [20] is included in the imaging domain. For compensating skullrelated aberration Jin et al. [13] have developed a strategy based on the approximation of raytracing; however, it may still suffer from the difficulty of accurately calculating Green's function when the velocity structure becomes more irregular.
Currently reversetime migration (RTM) has emerged as a more precise and powerful imaging tool in the exploration geophysics community [21][24]. RTM takes full advantage of the wave equation that includes all the dynamic features of a propagating wave field. Different from backprojection, RTM is based on the insensitivity of the wave equation's solution to the directionality of time. During RTM, by solving the wave equation with either the finitedifferences time domain (FDTD) method [25] or the pseudospectral time domain method [26], all transducers act as a virtual source by broadcasting their own records back to the domain in a timereversed manner. If the velocity model is precise, the reversedtime wave field should converge and be enhanced at the origins of the tobeimaged structures. Previous studies [19], [27], [28] have compared backprojection method and RTM by using complex velocity structure models, and have confirmed better results by using RTM. However, with the advantage of RTM outlined previously, the barriers to RTM adoption are also evident. Compared with backprojection, RTM appears to be more sensitive to accuracy of given velocity model. Once the velocity model is inaccurate, the quality of imaging results from RTM may be insignificant in comparison to backprojection [19], [27]. Consequently the solutions for creating accurate velocity model become extremely vital to guarantee the imaging quality of RTM.
For estimate velocity model correctly before RTM, to the best of author’s horizon there are two types of approaches. One is tomography, such as ultrasonic transmission tomography. Although it has been reported to be effectiveness for velocity correction from Jin and Wang in 2006 [29], prior of RTM detection, additional laboratory experiment has to be conducted for tomography data collection. This setup obviously increases the complexity of imaging procedure and financial cost. The other approach is an iterative imaging procedure developed by Whitmore [23]. In his approach for estimate the velocity correctly, an initial guess of the velocity model is setup by incorporating all known external velocity information as the first step. After that the TAT image result is derived by RTM, and a comparison between the reconstructed image and the velocity model is made. The differences derived by this comparison are attempted to be eliminated by generating the new velocity model. This procedure is repeated until both RTM results and velocity models are unchangeable. This iterative RTM procedure will be applied in this work.
Based on our previous studies [28], [30], [31] in this paper we closely concentrate on two topics: 1) comparison between backprojection and RTM on transcranial TAT; 2) Application of iterative RTM procedure on velocity correction. The synthetic dataset derived from a twodimensional (2D) brain model and real datasets acquired from laboratory experiments by Xu and Wang [9] (hereinafter referred as XW06) on rhesus monkey heads will be used for test case. The comparison of imaging results derived by both methods supports the finding of previous studies that RTM is superior to backprojection in imaging quality and accuracy. Meanwhile by analysis the imaging results derived by this iterative RTM procedure, conclusion shows that iterative RTM procedure should be a promising approach for achieving transcranial TAT.
This paper is organized as follows. Firstly backprojection and RTM are briefly introduced in Section 2. Section 3 illustrates the validation of RTM and backprojection methods through both synthetic data and laboratory data. The iterative RTM procedure will be described in detail and demonstrated with the applications to synthetic data and real laboratory data in Section 4. Section 5 is detailed discussion and analysis on the results coming from Section 3 and Section 4. Finally, in Section 6 we restate the major findings as conclusion.
2. Backprojection and reverse time migration
In this section we briefly introduce backprojection (or named Kirchhoff migration in exploration seismic engineering) and RTM, two different migration approaches used for transcranial imaging. Generally speaking migration is an inversion operation involving rearrangement of seismic or acoustic information elements from time to depth domain so that reflections and diffractions are plotted at their true locations [32]. It plays a function of moving dipping events to its original location and increasing the accuracy of the image. From seismic exploration the seismic wave is generated by the source on the ground to propagate downward the earth. Once it hits the noncontinuity along subsurface, reflection wave is generated and detected by transducers located on the ground [33]. The misfits exist between observed data and subsurface’s real location even sources and receivers are perfectly overlapped. The reasons include: 1) Due to the low frequency and noninfinitesimal width of the input beam, the reflection wave could propagate to multiple direction rather than following the single line; 2) Due to the possible nonhorizontal structure owned by subsurfaces, their locations could be “shifted” to wrong position. As Figure 1 shows, combinations of source and receiver are collocated (black solid dots) along the ground line AB. Suppose the velocity is uniform underground, the slope reflectors along CD is shifted incorrectly to C’D’ from the results collected by transducers with wave’s travel time unchanged. And therefore, the purpose of migration is to shift C’D’ back to CD, and hence the correct structure can be reached.
The migration problem can be transformed as a more straightforward case when it is applied on TAT imaging. The purpose of migration on TAT imaging is to shift the locations of different acoustic sources to their correct locations. In Figure 2, the shadowed area is the tobeimaged tissue, in which each point is an independent acoustic source after receiving microwave impulse. The dashed circle surrounded shows the location of receiver array. Suppose the velocity distribution in all 2D space is relatively uniform, and tissue’s acoustic wavelet is a single narrow pulse, the possible location of certain source S is constrained on a circle with the center as transducer T_{N}. Its radius is the distance calculated by the wave propagation time took from S to T_{N}. During the TAT imaging, for achieving the exact location of S, backprojection and RTM are applied by different principles. Backprojection can be summed as a kind of spatial integration method. It redistributes the received acoustic energies and makes summation in space. In Figure 2, the results observed by different transducers are redistributed along different circles in tobeimaged domain. If there is a large number of transducers, the amplitude at the location S will be enhanced due to interference of multi circles, otherwise the amplitudes should be minimum [33]. This approach is easily to achieve when the velocity model is uniform; However when the velocity distribution inside the tobeimaged domain is complex, the possible location of S will be distorted from circle to an irregular shape. Combined with multiple scattering, the assumption of oneway sight transmission in backprojection is no longer established. Different from backprojection, RTM works by running the wave equation backward for all transducers. It sets up all transducers as pseudo sources, and the records are broadcasted with a timereversed mannerfrom the end of the trace to time zero. If the velocity model is presumed correctly, in TAT imaging the wave field at time 0 in space will coincide with the original source distribution [19]. Compared with backprojection, by applying full wave equation, RTM is capable of involving all wave field phenomenon including diffraction, aberration and multiple scattering. However it requires much higher usage of CPU time and memory [27].
3. Validation of RTM and backprojection methods
3.1. Validation via synthetic data
To test the effectiveness of backprojection and RTM, we have built a 2D synthetic human brain model with an intact skull. In this model, the brain is made of gray matter and white matter [14], and the skull is made of three layers, namely the inner table, diploe, and outer table [20]. To mimic a real laboratory experiment similar to [8], [13], we modeled the space outside the skull as mineral oil. The reason for us using mineral oil is that it provides both background with uniform acoustic speed and tiny loss coupling between acoustic transducers and human head surface. Additionally when compared with water, mineral oil has much weaker microwave absorption ratio (~0), this will guarantee the minimum errors are introduced during modeling. Instead of mineral oil, a plenty kinds of fluid with weak microwave absorption ratio and comparable acoustic velocity to biotissue could also been used in further realistic test. To mimic pathological changes and build a benchmark for results analysis, in the synthetic model we replaced a small area of the brain with blood. This area was located at the left cerebral hemisphere and defined as ellipticallyshaped. The distribution of acoustic velocity in the synthetic model is shown in Figure 3a. From a review of several literatures, mechanical parameters of all related biotissues were collected and are listed in Table I, in which the velocities and densities of grey matter and white matter were measured from lamb brain using acoustic frequency of 1 MHz [14]. The skull’s velocities, densities and thicknesses of different layers are applied from datasets in research [20], [34] and [35]. The loss factor is defined in [35] to depict the amplitude of the propagating acoustic wave’s energy decay. Its values come from [14] for brain and [35] for skull, respectively.


Grey matter  1039  1483  0.0046   
White matter  1044  1521  0.0069   
Outer Table  1870  2900  0.1542  2.0 
Diploe  1740  2500  0.1234  2.5 
Inner Table  1910  2900  0.1985  1.7 
Blood  1057  1500  0   
Mineral oil  900  1437  0   
After establishing the mechanical properties, we calculated the initial acoustic amplitudes in the synthetic model. We assumed that the initial acoustic pressure (dyne/cm2) in mineral oil and the skull was 0, so the acoustic wave field was entirely generated by multiple acoustic sources in the brain at time zero. Their amplitudes were closely linked to the microwave absorption ratio. As expressed in [6], the relationship between the power intensity I (W/cm2) of absorbed microwaves and the generated peak acoustic pressure P0 (dyne/cm2) is shown as Eq. 1, where c is the sound velocity, β is the volumetric thermoexpansion coefficient, and C_{p} is heat capacity. The definition of absorbed power intensity I can be expressed as Eq. 2 [8], where E, ρ, and σ are the maximum amplitude of the radiated electromagnetic field, tissue’s density, and electrical conductivity, respectively. By combining Eq. 1 and Eq. 2, Eq. 3 is derived as the theoretical relationship between tissue’s electrical conductivity and initial acoustic pressure generated based on the thermoacoustic effect.
The volumetric thermoexpansion coefficient β and heat capacity Cp of the brain are nearly uniform: (β=12.3X10^{5} /℃ [36] and



White matter  0.665  1.081 
Grey matter  1.009  1.525 
Blood  1.868  2.283 
Once the synthetic model was established, we applied the FDTD method to forward modeling acoustic wave propagation. Using the initial acoustic pressure shown in Figure. 3b, a zerooffset Ricker wavelet with a central frequency of 0.15 MHz was applied to each point of the brain as the acoustic source. The model space was meshed as 512 x 512 grids with a spatial interval of 0.5 mm. A total of 1600 time steps with time intervals of 0.115 µs were judged long enough to allow the acoustic waves to propagate to each receiver from the most remote grid in the brain. The outgoing acoustic signal was recorded by 240 receivers located outside the skull (white circle shown in Figure. 3). The derived synthetic time traces are shown in Figure. 4, which provides the input dataset for later imaging.
The velocity model used in both migration algorithms is critical to successful imaging. In this study we applied two velocity models: one (abbreviated as V1) assumes the average acoustic velocity is uniformly 1540 m/s in the model space. Essentially it approximates a “bare brain” model with the effect of the skull excluded. The second model (abbreviated as V2) includes the effect of the high acoustic speed of the skull, which is almost double the speed of brain tissues. In our study, velocity models V1 and V2 were applied to both backprojection and RTM. Due to the velocity variance in V2, we applied V2 to two migration methods by different approaches. In backprojection, raytracing was applied from every transducer to all directions. This procedure was similar to the methods described in [13], but for simplicity we considered only the rays’ travel time caused by velocity variance, and aberration around the skull was ignored. Different from backprojection, as a kind of fullwave migration, RTM uses the same scheme as forward modeling methods such as FDTD. Consequently V2 can be applied to RTM in a straightforward manner. Comparisons of migration imaging results are shown in Figure. 57.
Figure 5a and Figure 6a are results derived by backprojection and RTM using the “bare brain” model V1. When compared with the distribution of initial acoustic pressure of the original model (Figure 5b or Figure 6b) we can clearly observe two kinds of imaging artifacts. First, the area with higher initial pressure at the left cerebral hemisphere, which is set up artificially in an ellipticalshape, is seriously enlarged by backprojection (Figure 5a) and falsely elongated along the major axis of the ellipse by RTM (Figure 6a); Second, delicate features such as the gyrus, located in the outer part of the brain, and the gap which separates the left and right cerebral hemispheres in our model are totally blurred in both results when using velocity model V1. In contrast, from Figure 5c and Figure 6c, which are results based on velocity model V2 with the skull’s velocity included, these two misfits are substantially reduced. From all of these comparisons we can see that exclusion of the skull leads to severe error and distortion in migration imaging for both backprojection and RTM.
To further examine the differences between backprojection and RTM, we reorganized the results using backprojection, RTM, along with the original velocity model V2 as shown in Figure 7. From it we can see that although both backprojection and RTM can transform most of the wave field back to its original location correctly, there are obvious differences in imaging quality between the two methods. Compared with the original model shown in Fig. 5b, the detail features are blurred in backprojection result (Figure 7a) but appear to be clean and sharp in RTM result (Figure 7c). These visual differences can be further amplified through 1D comparison along the xdirection along the horizontal white dashed line shown in Figure 7 chosen to cross the brain’s left boundary, artificially blooded area, interhemispheric fissure, and right boundary. The source amplitudes along this profile for backprojection RTM, and the original model are shown in Figure 8. It is obvious that the result of RTM is far more superior than that of backprojection. Compared with the backprojection result (the dotted line), result of RTM (the solid line) has larger variance for depicting structures such as interhemispheric fissure around 12.50 cm as well as the boundaries of the artificially blooded area around 7.5 cm and 10 cm. These differences can also be clearly observed on brain boundaries around 5 cm and 20 cm, where the amplitude from RTM decays as sharp as the original model, but backprojection’s result is obviously incorrect and decay much slower outside of the brain.
3.2. Validation through laboratory data
The backprojection and RTM algorithms were also tested by using the laboratory data acquired by XW06. In their experiment, the monkey’s head was decapitated and fixed by a clamp and completely immersed in mineral oil. During TAT detection, this specimen was stimulated by 3GHz microwave pulses, and the derived acoustic wave field was recorded by a transducer with a 1 MHz central frequency and about 0.8 MHz bandwidth. The transducer was positioned from 614 cm to the center of the monkey’s head, and the sampling frequency was 20 MHz. During the experiment, the clamp fixing the monkey’s head was mounted on a rotary table driven by a stepper motor with a step size of 2.25 degrees. Accordingly in this laboratory application, the outgoing acoustic wave was observed by 160 receivers surrounding the head in a 2D circle. With data processing performed through the procedure of [38] for high frequency enhancement, only the segment with a spectrum of 0.31 MHz of the observed data was picked up and enhanced for imaging by backprojection and RTM. We applied the estimated average acoustic velocity to both image approaches, since any velocity information on velocity distribution of earlier experiment in XW06 was unknown, which may introduce some error and reduce the image quality.
Figure 9 shows the results based on a dataset collected from a onemonthold monkey head with a skull thickness of less than 1 mm. The velocity model is assumed to be uniform as 1437 m/s, as the same as acoustic velocity of mineral oil. The region shown is 53mm by 51mm along the coronal cross section. From the experiment of XW06, three steel needles with diameters of 0.9 mm were inserted in the approximate locations as shown in Figure. 9a (XW06). The results derived from backprojection (Figure. 9b) and RTM (Figure. 9c) are shown sidebyside for comparison. Both imaging algorithms show the three needles, and the black dot located at the center is believed to be an air bubble introduced by inserting the needles (XW06). Compared with the result derived from backprojection, the needle A in RTM result owns sharper edges. Meanwhile, backprojection provides less visibility for needle C than RTM. From the plots along the x crosssection shown in Figure. 9d, although both needle A and B can be detected by using backprojection and RTM, the backprojection image is much nosier. The existence of this noise causes seriously reduction of signal noise ratio and fussy in whole image by backprojection. However needle B is seriously distorted from results of both backprojection and RTM. This shows that due to the technical limitations of the coarsely estimated average acoustic velocity for TAT reconstruction, neither of these two methods can provide satisfying image quality.
4. Iterative RTM procedure
The previous result shows that the image quality of RTM is restricted by the precision of the velocity model. In some sense, the increased requirement of RTM sensitivity to indistinct velocity model brings a paradox: If the velocity model is well known by operation or any invaded techniques, RTM is not necessary; If the velocity model is not known, RTM lacks an essential input. For solving this paradox and improving the RTM result, rather than single RTM run, we developed iterative RTM procedure based on the theory of Whitmore [23]. The iterative RTM procedure works based on an assumption between velocity model and RTM result. Suppose there exists a function to map velocity from RTM result, complex velocity model could be renewed after RTM for one time. By repeating this procedure, velocity model will continue be updated after multiple RTM runs until the result becomes promising. This iterative RTM procedure is shown in Figure 10 as a schematic. An initial guess of a model is made, incorporating all known external velocity information. The image is derived through RTM. After that a new velocity model is made based on current RTM result. This procedure is repeated until both RTM results and velocity models are unchangeable.
As Figure 10 shows, the key procedure of iterative RTM is to update velocity model after each RTM. This requires a function for mapping velocity from image derived by RTM. Since for TAT imaging RTM image is the reconstruction of electrical conductivity distribution, the desired function of velocity and electrical conductivity could be built by referring to Table 1 and Table 2. Based on the given velocity and electrical conductivity values of white matter and gray matter, a function is established by using the second order polynomial fitting as shown in Figure 11. Obviously this function is coarsely approximated. It could be improved by involving more velocity–electrical conductivity pairs in further research. Meanwhile it’s noteworthy that the mapping function can’t derive correct velocity on blood, water and mineral oil, consequently a mask must be used to preserve the region outside sample from wrongly updating. At present the function shown in Figure 11 is used in our research to map the velocity from electrical conductivity derived by RTM; in this way the velocity model will be updated for the next RTM iteration.
An example of applying iterative RTM procedure is shown in Figure 12. The laboratory data of this case comes from XW06, which has been introduced in Section V. By referring to XW06, the velocity model would not be known during TAT. The only knowledge of the velocity field in this case is 1437 m/s of mineral oil, which submerges the sample during TAT detection. To demonstrate the iterative RTM procedure, we make an initial guess that velocity model is uniform as 1437 m/s. By using the initial velocity model, the first RTM output is derived as shown in Figure 12b, which is as the same as Figure 9c. Velocity model will then be updated as new velocity input for RTM. This procedure is repeated for several times until RTM output is unchangeable. The TAT results derived from iterative RTM procedure on step 3 and 5 are displayed on Figure 12b and Figure 12c.
From Figure 12, by using iterative RTM procedure the improvement on TAT image can be observed from several aspects: First, by comparing with Figure 12b (Figure 9c), the distortion of needle B has been well corrected after five times iteration. Second, the original RTM result shown in Figure 12b provides less visibility for needle C. After two iterations, Figure 12c shows that the boundary of needle C has been largely enhanced. Additionally, comparing among Figure 12 ac, the crosssection of needle A has been focused gradually with iteration times increasing. All of this shows that iterative RTM procedure has great potential of correcting coarsely estimated velocity model and therefore enhance imaging capability of RTM.
5. Discussion
By combining Section II and III, It is clear that RTM is superior to backprojection in terms of imaging quality and higher signal to noise ratio. Compared with backprojection, which makes a ray approximation, RTM bases its entire algorithm on solving a fullwave equation, without substantial approximation, and holds the original dynamic features of the wave field intact. The handle of velocity heterogeneity is fundamentally intrinsic. Usually, adapting ray tracing in backprojection is time consuming and has a limited improvement on image quality. Unlike backprojection, the quality of RTM’s results is independent of the complexity of the velocity model. This feature makes the wave propagation in the domain highly accurate compared with using raytracing. Figure 58 show that RTM is able to recover almost all features of the brain to their original potion when the skull’s velocity is included.
RTM can recover a complex structure’s boundary sharply. This has been reported by [19], [27] and is proved by our results in Figure 6. When both velocity models V1 and V2 are used in RTM, the boundaries of tiny features can be clearly seen. Especially, even though obvious distortions exist in the result using V1 (Figure 6a), with the exclusion of the skull, all features are still relatively unblurred in comparison with the backprojection results (Figure 5a). By comparison, looking at the crosssection in Figure 8, the sharp edges of brain are well recovered by RTM but seriously smeared by backprojection. Further, in Figure 9 when the skullexcluded model is applied, the contour of Needle C is well recovered by RTM but not by backprojection.
Nevertheless, it is noteworthy that the image quality of RTM is still limited by the precision of the velocity model. Consequently the key to capitalizing on the benefit of RTM is to build better velocity models before applying RTM. For achieving this goal, we developed iterative RTM procedure to update velocity model iteratively. Its principle is based on an assumption that acoustic velocity could be mapped from RTM image through a certain function. The results shown in Figure 12 demonstrated that currently even by using our coarsely approximated velocity function, the quality of RTM image can be well improved after several iterations within iterative RTM procedure. It could be improved by involving more velocity–electrical conductivity pairs in our further research.
6. Conclusion
In this paper we have compared backprojection and RTM for transcranial TAT imaging. Compared with backprojection, RTM offers better performance with regard to velocity variance, imaging quality, and noise suppression caused by spatial aliasing. The capability of RTM can be further improved by iterative RTM procedure, which aimed to provide velocity model with higher accuracy. Generally speaking RTM owns great potential for achieving acoustic localization on transcranial TAT imaging with high quality and accuracy.
References
 1.
M. Xu and L. V. Wang, “Timedomain reconstruction for thermoacoustic tomography in a spherical geometry,” IEEE Trans. Med. Imag. , vol. 21, pp. 814822, Jul. 2002.  2.
Y. Xu, D. Feng and L. V. Wang, “Exact frequencydomain reconstruction for thermoacoustic tomographyI: Planar Geometry,” IEEE Trans. Med. Imag. , vol. 21, pp. 823828, July 2002.  3.
Y. Xu, D. Feng and L. V. Wang, “Exact frequencydomain reconstruction for thermoacoustic tomographyII: Cylindrical Geometry,” IEEE Trans. Med. Imag. , vol. 21, pp. 829833, July 2002.  4.
W. Joines, R. Jirtle, M. Rafal, and D. Schaeffer, “Microwave power absorption differences between normal and malignant tissue,” Radiation Oncol. Biol. Phys ., vol. 6, pp. 681687, 1980.  5.
M. S. Hawley, A. Broquetas, L. Jofre, J. C. Bolomey, and G. Gaboriaud, “Microwave imaging of tissue blood content changes,” J. Biomed. Eng. , vol. 13, pp. 3744, 1991.  6.
K. R. Foster and E. D. Finch, “Microwave hearing: evidence for thermoacoustic auditory stimulation by pulsed microwaves,” Science , vol. 185, pp. 256258, 1974.  7.
R. Nikolaus, “Thermoacoustics,” Adv. Appl. Math. , vol. 20, pp. 135175, 1980.  8.
M. MartinezBurdalo, A. Martin, M. Anguiado and R. Villar, “Comparison of FDTDcalculated specific absorption rate in adults and children when using a mobile phone at 900 and 1800 MHz,” Phys. Med. Biol. , vol. 49, pp. 345354, 2004.  9.
Y. Xu and L. V. Wang, “Rhesus monkey brain imaging through intact skull with thermoacoustic tomography,” IEEE Trans. Ultrason. Ferroelectr. Freq. Control. , vol. 53, pp. 542548, Mar., 2006.  10.
D. Royer, “Mixed matrix formulation for the analysis of lasergenerated acoustic waves by a thermoelastic line source,” Ultras. Internat. , vol. 39, pp 345354, 2001.  11.
N.H. Noroy, D. Royer, M. Fink, “Transient elastic wave generation by an array of thermoelastic sources,” Appl. Phys. Lett , Vol.63, pp. 32763278, Dec 1993.  12.
J. Gamelin, A. Aguirre, A. Maurudis, F. Huang, D. Castillo and Q. Zhu, “Curved array photoacoustic tomography system for small animal imaging,” J. Biomed. Opt ., vol. 13, no.2, Mar.Apr., 2008.  13.
X. Jin, C. Li and L. V. Wang, “Effects of acoustic heterogeneities on transcranial brain imaging with microwaveinduced thermoacoustic tomography,” Med. Phys. , vol. 35, no. 7, 2008.  14.
SC. Lin, SJ. Shieh and M. J. Grimm, “Ultrasonic measurements of brain tissue properties,” Center for Disease Control conf. Preproceedings , pp. 2731,1997.  15.
W. A. Schneider, “Developments in seismic data processing analysis,” Geophysics , vol. 36, pp 10431073, 1971.  16.
W. A. Schneider, “Integral formulation for migration in two and three dimensions,” Geophysics , vol. 43, pp 4976, 1978.  17.
W. S. French, “Computer migration of oblique seismic reflection profiles,” Geophysics , vol. 51, pp. 961980, 1975.  18.
A. J. Berkout, Seismic migration imaging of acoustic energy by wave field extrapolation , Elseiver Science Publ. Co., Inc., 1982.  19.
J. Zhu and L. R. Lines, “Comparison of Kirchhoff and reverse time migration methods with applications to prestack depth imaging of complex structures,” Geophysics , vol. 63, pp. 11661176, 1998.  20.
F. J. Fry and J. E. Barger, “Acoustic Properties of Human Skull,” J. Acoust. Soc. Am. , vol. 63, pp. 15761590, 1978.  21.
C. Hemon, “Equations d'onde et modeles,” Geophys. Prosp , vol. 26, pp. 790821, 1978.  22.
G. A. McMechan, “Migration by extrapolation of timedependent boundary values,” Geophys. Prosp. , vol.31, pp. 412420, 1983.  23.
N. D. Whitmore, “Iterative depth migration by backward time propogation,” 53^{rd} Ann. Internat. Mtg. Of Soc. Expl. Geophys. , Expanded abstracts, pp. 827830, 1983.  24.
E. Baysal, D. D. Kosloff and J. W. C. Sherwood, “Reversetime migration,” Geophysics , vol. 48, pp. 15141524, 1983.  25.
K. Yee, “Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media,” IEEE Trans. Antennas Propag. , vol. 14, 1966.  26.
Liu, Q. H., “The PSTD algorithm: A timedomain method requiring only two cells per wavelength,” Microw. Opt. Technol. Lett. , Vol. 15, No. 3, Jun. 1997.  27.
P. Farmer, Z. Zheng and D. Jones, “The role of reverse time migration in imaging and model estimation,” The leading edge, Soc. Expl. Geophys , vol. 28, pp. 436441, 2009.  28.
Z. Liu and L. Liu, “Transcranial Thermoacoustic Tomography: A Comparison of Two Imaging Algorithms,” IEEE Trans. Med. Imag. , Vol 32, pp. 289294, 2013.  29.
X. Jin and L. V. Wang, “Thermoacousitc tomography with correction for acoustic speed variations,” Phys. Med. Biol. , vol. 51, pp. 64376448, 2006.  30.
Z. Liu and L. Liu, “A Novel Approach for Thermoacoustic Tomography by Kirchhoff Migration,” 9th International Conference on Theoretical and Computational Acoustics , Germany, 2010.  31.
L. Liu, K. He and L. V. Wang, “Transcranial ultrasonic wave propagation simulation: skull insertion loss and recovery,” in SPIE conf. , San Francisco, USA, 2007.  32.
R. E. Sheriff, Encyclopedic dictionary of applied geophysics , fourth edition, Society of Exploration Geophysicists, 2002.  33.
R. E. Sheriff and L. P. Geldart, Exploration seismology , second edition, Cambridge University Press, 1995.  34.
N. Lynnerup, J. G. Astrup and B. Astrup, “Thickness of the human cranial diploe in relation to age, sex and general body build,” Head and Face Medicine , 113, 2005.  35.
M. Hayner and K. Hynynen, “Numerical analysis of ultrasonic transmission and absorption of oblique plane waves through a human skull,” J. Acoust. Soc. Am. , vol. 110, pp. 33193330, Dec., 2001.  36.
J. C. Lin, “Microwaveinduced hearing: some preliminary theoretical observations,” J. Micro. Power , vol. 11, No. 3, 1976.  37.
X. Xu, P. Tikuisis and G. Giesbrecht, “A mathematical model for human brain cooling during coldwater neardrowning,” J Appl Physiol , vol. 86, pp. 265272, 1999.  38.
Y. Xu and L. V. Wang, “Signal processing in scanning thermoacoustic tomography in biological tissues,” Med. Phys. , vol. 28, pp. 15191524, Jul. 2001.