Photoacoustic tomography (PAT) is a hybrid ultrasound-mediated biophotonic imaging modality that combines the high contrast of optical imaging with the high-spatial resolution of ultrasound imaging.1–3 In PAT, biological tissues are illuminated with short laser pulses that results in the generation of internal acoustic wavefields via the thermoacoustic effect. The initial amplitudes of the induced acoustic wavefields are proportional to the spatially variant absorbed optical energy density in the tissue. The acoustic waves propagate through the tissue and are subsequently detected by use of a collection of wide-band ultrasonic transducers that are located outside the object. From knowledge of these data, a reconstruction algorithm is employed to estimate the initial induced pressure distribution or, equivalently if the Grueneisen parameter is known,2,4 the absorbed optical energy distribution within the tissue. Because optical contrast is dependent on hemoglobin concentration and related to the molecular constitution of tissue, PAT can reveal the pathological condition of the tissue,5,6 and therefore facilitate a wide-range of diagnostic tasks.7–11
Transcranial brain imaging represents an important application that may benefit significantly by the development of PAT methods. Existing high-resolution human brain imaging modalities such as x-ray computed tomography (CT) and magnetic resonance imaging (MRI) are expensive and employ bulky and generally nonportable imaging equipment. Moreover, x-ray CT employs ionizing radiation and is therefore undesirable for use with patients who require frequent monitoring of brain diseases or injuries. Ultrasonography is an established portable pediatric brain imaging modality, but its image quality degrades severely when employed after the closure of the fontanels. The photoacoustic (PA) signals recorded in a PAT experiment experience only a one-way transmission through the skull. Accordingly, they are generally less attenuated and aberrated than the echo data recorded in transcranial ultrasound imaging, which are contaminated by the effects of a two-way transmission through the skull. Moreover, much of the broadband PA signal energy resides at frequencies less than 1 MHz, and these relatively low-frequencies interact less strongly with skull bone12 than do higher frequency ultrasound beams that are typically employed in pure ultrasound imaging.
Transcranial PAT studies have revealed structure and hemodynamic responses in small animals13–15 and anatomical structure in human infant brains have been conducted.15–19 Because the skulls in those studies were relatively thin (), they did not significantly aberrate the PA signals and conventional back-projection methods were employed for image reconstruction. However, PA signals can be significantly aberrated by thicker skulls present in adolescent and adult primates. To render PAT an effective modality for use with transcranial imaging in large primates, including humans, it is necessary to develop image reconstruction methodologies that can accurately compensate for skull-induced aberrations of the recorded PA signals.
Toward this goal, Xin et al.20 proposed an image reconstruction method that sought to compensate for PA signal aberration associated with acoustic wave reflection and refraction within the skull. In that method, the skull was assumed to be acoustically homogeneous. Accordingly, the method did not explicitly account for scattering effects that arose from heterogeneities in the skull. As a result of the simplified skull model employed, only modest improvements in image quality were observed as compared to use of a standard back-projection-based reconstruction algorithm. Therefore, there remains an important need for the development of improved image reconstruction methodologies for transcranial PAT that are based upon more accurate models of the skull’s heterogeneous acoustic properties.
In this work, we propose and investigate a reconstruction methodology for transcranial PAT that employs detailed subject-specific descriptions of the acoustic properties of the skull to mitigate skull-induced blurring and distortions in the reconstructed image. The reconstruction methodology is comprised of two primary steps. In the first step, the spatially varying speed-of-sound (SOS) and mass density distributions of the to-be-imaged subject’s skull are determined by use of adjunct x-ray CT data. This is accomplished by use of a method that was developed previously to facilitate transcranial adaptive acoustic focusing for minimally invasive brain surgery.21 In the second step, the subject-specific SOS and density distributions are employed with a time-reversal image reconstruction method22 for estimation of the spatially variant initial amplitude of the thermoacoustically induced pressure signals within the brain.
The paper is organized as follows. In Sec. 2, the salient imaging physics and image reconstruction principles are briefly reviewed. The image reconstruction methodology is in Sec. 3, which includes a description of how the SOS and density maps of a skull are computed from adjunct x-ray CT data. Section 4 describes the results of the image reconstruction studies that employ a well-characterized phantom and a primate brain, both enclosed in a skull. The paper concludes with a summary and discussion of future work in Sec. 5.
Below we briefly review the imaging physics and image reconstruction principles relevant to transcranial PAT imaging. The reader is referred to Refs. 1–3 and 23 for comprehensive reviews of PAT image formation principles.
The method presented in this work accounts for refraction and diffraction of the photoacoustic wavefields caused by heterogeneities in the SOS and density distributions of the skull. The effects of shear waves propagating in the skull and attenuation of the photoacoustic wavefield in the skull on image formation have not been incorporated in this initial study and represent a topic for future study as discussed in Sec. 5. It should be noted that the density and SOS heterogeneities both have stronger effects and cause more acoustic energy to be reflected back into skull than do the shear wave mode-conversion or attenuation for acoustic plane-wave components incident on the skull surface at small ( or less) angles24 from the surface normal.
Consider tissue that is irradiated by the laser pulse at time and let denote the induced pressure wavefield at time at location . When acoustic attenuation and shear wave mode-conversion can be neglected, the following three coupled equations describe the forward propagation of in heterogenous acoustic media:25
Image Reconstruction Based on Time-Reversal
Let denote the measured pressure wavefield data, where denotes the location of the ’th ultrasonic transducer, indicates the number of measurement locations, and denotes the maximum time at which the pressure data were recorded. The PAT image reconstruction problem we address is to obtain an estimate of or, equivalently, , from knowledge of , , and . The development of image reconstruction methods for addressing this problem is an active area of research.26–28
In this work, we employed a k-space pseudospectral time-reversal image reconstruction algorithm that has been developed by Treeby and Cox.29 The reconstruction algorithm approximates Eqs. (1)–(3), neglecting terms proportional to , and solves discretized versions of the three coupled lossless acoustic backward in time with initial and boundary conditions specified as:
Unlike the approach proposed by Xin et al.,20 this algorithm is based on a full-wave solution to the acoustic wave equation for heterogeneous media, and can therefore compensate for scattering due to variations in and within the skull. Because the fast Fourier transform (FFT) algorithm is employed to compute the spatial derivatives in the k-space pseudospectral method,29 its computational efficiency is superior to competing methods such as finite-difference time-domain methods. The accuracy of the temporal derivatives can also be refined by use of k-space adjustments.22
Image Reconstruction Studies
Our methodology for aberration correction in transcranial PAT image reconstruction is described below along with details of our experimental studies that involved monkey skulls.
The reconstruction methodology is comprised of two primary steps. First, the spatially varying SOS and density distributions of the to-be-imaged subject’s skull are determined by use of adjunct x-ray CT data. These distributions are subsequently employed with the time-reversal image reconstruction method22 described above for estimation of within the brain tissue from knowledge of the measured data .
Description of Biological Phantoms
Two biological phantoms were employed in the experimental studies. The first phantom was the head of an eight-month old rhesus monkey that was obtained from the Wisconsin National Primate Research Center. The hair and scalp were removed from the skull. A second, more simple, phantom was constructed by removing the brain of the monkey and replacing it by a pair of iron needles of diameter 1 mm that were embedded in agar. This was accomplished by cutting off the calvaria to gain access to the brain cavity.
Estimation of the Skull’s SOS and Mass Density Distributions from CT Data
The wavefront aberration problem encountered in transcranial PAT is conjugate to one encountered in transcranial focusing of high-intensity ultrasound30–32 for therapy applications. Both problems involve a one-way propagation of ultrasound energy through the skull and both require that the wavefront aberrations induced by the skull be corrected. The problems differ in the direction of the propagating acoustic wavefields. The feasibility of utilizing skull information derived from adjunct x-ray CT image data to correct for wavefield aberrations in transcranial focusing applications has been demonstrated.21 As described below, we adopted this method for determining estimates of and , characterizing the acoustic properties of the subject’s skull, from adjunct x-ray CT data.
As described by Aubry et al.,21 the SOS and density maps of the skull can be estimated from a porosity map using mixture laws in a biphasic medium (bone/water). Let denote the value of the ’th voxel in the x-ray CT image, which is measured in Hounsfield units. A voxel-based representation of the porosity map, denoted as , can be established from knowledge of as21,33
Let and denote voxel-based representations of the skull’s mass density and SOS distributions. The density map can be estimated from the porosity map as12,21 According to Carter and Hayes,34 the elastic modulus of bone is proportional to the apparent density cubed as a first order approximation. This suggests a linear relationship between the speed-of-sound and the porosity: 12,21
The monkey skull that was part of both phantoms described in Sec. 3.1 was imaged using an x-ray CT scanner (Philips Healthcare, Eindhoven, The Netherlands) located at Washington University in St. Louis. Details regarding this system can be found in Ref. 35. Prior to imaging, three fiducial markers were attached to the skull to facilitate co-registration of the determined SOS and density maps with the reference frame of the PAT imaging system. The three fiducial markers [see Fig. 1(a)] were iron balls of diameter of 1.5 mm and were carefully attached to the outer surface of the skull. The fiducial markers were located in a transverse plane that corresponded to the to-be-imaged 2-D slice in the PAT imaging studies described below. In the x-ray CT studies, the tube voltage was set at 130 kV, and a tube current of 60 μA was employed. Images were reconstructed on a grid of 700 by 700 pixels of dimension . This pixel size is much less than the smallest wavelength (0.5 mm) detected by the ultrasound transducer used in the PAT imaging studies described below. This precision is sufficient to accurately model acoustic wave propagation in the skull by using the k-space pseudospectral methods.36 The reconstructed CT image is displayed in Fig. 1(a).
From knowledge of the CT image, the porosity map was computed according to Eq. (6). Subsequently, the density and SOS maps and were computed according to Eqs. (7) and (8). Images of the estimated and maps are displayed in Fig. 1(b) and 1(c). To corroborate the accuracy of the adopted method for estimating the skull’s SOS and density distributions from x-ray CT data, i.e., Eqs. (7) and (8), direct measurements, of the skull’s average SOS along ray-paths perpendicular to the skull surface at five locations were acquired. This was accomplished by use of a photoacoustic measurement technique depicted in Fig. 2(a). Additionally, the average density of the skull was computed and compared to the average computed from the values estimated from the x-ray CT data. The results show that the directly measured average SOS and density are very close (about 1 to 6%) to the estimated values from CT data. These results corrorborate the adopted method for estimating the skull’s SOS and density distributions from adjunct x-ray CT data. Details regarding these studies are contained in the Appendix.
PAT Imaging Studies: Data Acquisition
After the skull’s SOS and density distributions were estimated from the adjunct x-ray CT data, the two phantoms (that included the skulls) were imaged by use of the PAT imaging system shown in Fig. 2(b). Images of the two phantoms with the skull removed, i.e., images of the extracted monkey brain and crossed needles embedded in agar, were also acquired, which will serve as control images. The imaging system employed a 2-D scanning geometry and has been employed in previous studies of PAT imaging of monkey brains.37 The imaging plane and fiducial markers were chosen to be approximately 2 cm below the top of the skull, such that the imaging plane was approximately normal to the skull surface at that plane. The phantoms (crossed needles and the primate cortex) were moved to the imaging plane, so that the amount of acoustic energy refracted out of the imaging plane was minimized. Additionally, the system was aligned to ensure the scanning plane and the imaging plane coincided.
The phantoms were immersed in a water bath and irradiated by use of a tunable dye laser from the top (through the skull for the cases when it was present) to generate PA signals. The laser (NS, Sirah), which was pumped by a Q-switched Nd:YAG laser (PRO-35010, Newport) operating at a wavelength of 630 nm with a pulse repetition rate of 10 Hz, was employed as the energy source. The laser beam was expanded by use of a piece of concave lens and homogenized by a piece of ground glass before illuminating the target. The energy density of the laser beam on the skull was limited to (within the ANSI standard), which was further attenuated and homogenized by the skull before the laser beam reached the object.
As shown in Fig. 2, a circular scanning geometry with a radius of 9 cm was employed to record the PA signals. A custom-built, virtual point ultrasonic transducer was employed that had a central frequency of 2.25 MHz and a one-way bandwidth of 70% at . Additional details regarding this transducer have been published elsewhere.37 The position of the transducer was varied on the circular scan trajectory by use of a computer-controlled step motor. The angular step size was 0.9 degrees, resulting in measurement of at locations on the scanning circle. The total acquisition time was approximately 45 min.
The PA signals received by the transducer were amplified by a 50-dB amplifier (5072 PR, Panametrics, Waltham, MA), then directed to a data-acquisition (DAQ) card (Compuscope 14200; Gage Applied, Lockport, IL). The DAQ card was triggered by the Q-switch signal from the laser to acquire the photoacoustic signals simultaneously. The DAQ card features a high-speed 14-bit analog-to-digital converter with a sampling rate of . The raw data transferred by the DAQ card was then stored in the PC for imaging reconstruction.
PAT Imaging Studies: Image Reconstruction
Image reconstruction was accomplished in two steps: (1) registration of the SOS and density maps of the skull to the PAT coordinate system and (2) utilization of a time-reversal method for PAT image reconstruction in the corresponding acoustically heterogeneous medium.
The estimated SOS and density maps and were registered to the frame-of-reference of the PAT imaging as follows. From knowledge of the PAT measurement data, a scout image was reconstructed by use of a half-time reconstruction algorithm.38 This reconstruction algorithm can mitigate certain image artifacts due to acoustic aberrations, but the resulting images will, in general, still contain significant distortions. The PAT image of monkey head phantom (with brain present) reconstructed by use of the half-time algorithm is displayed in Fig. 1(d). Although the image contains distortions, the three fiducial markers are clearly visible. As shown in Fig. 1(a), the fiducial markers were also clearly visible in the x-ray CT image that was employed to estimate the SOS and density maps of the skull. The centers of the fiducial markers in the X-ray CT and PAT images were determined manually. From this information, the angular offset of the x-ray CT image relative to the PAT image was computed. The SOS and density maps were downsampled by a factor of two, to match the pixel size of the PAT images, and rotated by this angle to register them with the PAT images.
The re-orientated SOS and density maps were employed with the k-space time-reversal PAT image reconstruction algorithm developed by Cox et al.28 The numerical implementation of this algorithm provided in the Matlab k-Wave Toolbox29 was employed. The measured PA signals were pre-processed by a curvelet denoising technique prior to application of the image reconstruction algorithm.39 The initial pressure signal, , was reconstructed on a grid of of dimension 0.2 mm. For comparison, images were also reconstructed on the same grid by use of a back-projection reconstruction algorithm.40 This procedure was repeated to reconstruct images of both phantoms and the corresponding control phantoms (phantoms with skulls removed).
Images of Needle Phantom
The reconstructed images corresponding to the head phantom containing the needles are displayed in Fig. 3. Figure 3(a) displays the control image of the needles, without the skull present, reconstructed by use of the back-projection algorithm. Figure 3(b) and 3(c) displays reconstructed images of the phantom when the skull was present, corresponding to use of back-projection and time-reversal reconstruction algorithms, respectively. All images have been normalized to their maximum pixel value, and are displayed in the same gray-scale window. Due to the skull-induced attenuation of the high-frequency components of the PA signals, which was not compensated for in the reconstruction process, the spatial resolution of the control image in Fig. 3(a) appears higher than the images in Fig. 3(b) and 3(c). However, the image reconstructed by use of the time-reversal algorithm in Fig. 3(c) contains lower artifact levels and has an appearance closer to the control image than the image reconstructed by use of the back-projection algorithm in Fig. 3(b). This is expected, since the time-reversal algorithm compensates for variations in the SOS and density of the skull while the back-projection algorithm does not.
These observations are corrorborated by examination of profiles through the three images shown in Fig. 3(d), which correspond to the rows indicated by the superimposed dashed lines on the images. The solid black, dotted blue, and dashed red lines correspond to the reconstructed control image, and images reconstructed by use of the back-projection and time-reversal algorithms, respectively. The average full-width-at-half-maximum of the two needles in the images reconstructed by use of the time-reversal algorithm is reduced by 8% compared to the corresponding value computed from the images obtained via the back-projection algorithm.
Images of Monkey Brain Phantom
The reconstructed images corresponding to the head phantom containing the brain are displayed in Fig. 4. Figure 4(a) displays photographs of the cortex and outer surface of the skull. Figure 4(b) displays the control image (skull absent) reconstructed by use of the back-projection algorithm. The images of the complete phantom (skull present) reconstructed by use of the back-projection and time-reversal algorithms are shown in Fig. 4(c) and 4(d), respectively. All images have been normalized to their maximum pixel value, and are displayed in the same gray-scale window. As observed above for the needle phantom, the brain image reconstructed by use of the time-reversal algorithm in Fig. 4(d) contains lower artifact levels and has an appearance closer to the control image than the image reconstructed by use of the back-projection algorithm in Fig. 4(c).
This observation was quantified by computing error maps that represented the pixel-wise squared difference between the control and reconstructed images with the skull present. Figure 5(a) and 5(b) displays the error maps between the control image and the images reconstructed by use of the back-projection and time-reversal algorithms, respectively. The error maps were computed within the region interior to the skull, which is depicted by the red contours superimposed on Fig. 4(b)–4(d). Additionally, the root mean-squared difference (RMSD) was computed by computing the average values of the difference images. The RMSD corresponding to the back-projection and time-reversal results were 0.085 and 0.038. These results confirm that the image reconstructed by use of the time-reversal method, which compensated for the acoustic properties of the skull, was closer to the control image than the image produced by use of the back-projection algorithm.
Summary and Discussion
We have investigated a reconstruction methodology for transcranial PAT that employs detailed subject-specific descriptions of the acoustic properties of the skull to mitigate skull-induced distortions in the reconstructed image. Adjunct x-ray CT image data were employed to infer the spatially variant SOS and density distributions of the skull. Knowledge of these quantities was employed in a time-reversal image reconstruction algorithm to mitigate skull-induced aberrations of the measured PA signals. Our preliminary experimental results, which employed a primate skull, demonstrated that the reconstruction methodology can produce images with improved fidelity and reduced artifact levels as compared to a previously employed back-projection algorithm. This is an important step toward the application of PAT for brain imaging in human subjects.
The use of x-ray CT image data for estimating the skull’s SOS and density distributions was motivated by previous studies of transcranial ultrasound focusing.21 Assuming that the skull size and shape does not change, only a single CT scan is required to estimate the SOS and density maps, and does not need to be repeated for subsequent PAT imaging studies of that patient. Because of this, it may be possible to safely monitor brain injuries or conduct other longitudinal studies without repeated exposure to ionizing radiation. Moreover, it may be possible to use adjunct image data produced by alternative modalities such as magnetic resonance imaging41 or ultrasound tomography42 to estimate the required SOS and density maps.
There remain several important topics for future study that may further improve image quality in transcranial PAT. In this preliminary study, to accommodate the 2-D PAT imaging system and 2-D image reconstruction algorithm, the phantoms were moved to the image plane, approximately 2 cm below the top of the skull. Note that this is not the plane in which the cortical vessels are normally found; the primate brain was moved to align the cortical vessels with the imaging plane. For in-vivo transcranial PAT applications in which the cortical structure is of interest, the geometry of the skull necessitates a full 3-D treatment of the image reconstruction problem. The development of robust image reconstruction algorithms for this task is currently underway. Additionally, accurate measurement of the transducer’s electrical impulse response (EIR) and subsequent deconvolution of its effect on the measured data11 would improve image reconstruction accuracy. The possibility also exists that alternatives26,43 to the time-reversal image reconstruction algorithm employed in this study may yield improvements in image quality.
In terms of the imaging physics, it is expected that development and utilization of image reconstruction algorithms that can compensate for the effects of acoustic attenuation and shear wave mode-conversion will further improve image resolution, particularly for thicker skulls. To date, the manner in which shear waves propagating in the skull affect PAT has only been investigated quantitatively for stratified planar media and planar detection surfaces24,44 via computer-simulation studies. In that case, the effects were observed only in certain high-spatial resolution components of the imaged object. The effects of shear wave propagation and attenuation are both strongly dependent on the thickness of the skull. In this work, the average thickness of the skull was 3 mm and the distortions to the PA signal due to absorption and propagating shear waves were expected24 to be of second-order effect as compared to the distortions due to the variations in the SOS and density. For adult human skulls, where the skull can be thick, the relative importance of these effects in transcranial PAT remains to be investigated.
Validation of Speed-of-Sound and Density Maps
The density and speed-of-sound of the skull were also directly measured to corrorborate the accuracy of the adopted method for estimating the skull’s SOS and density distributions from x-ray CT data. By using the water displacement method, the measured average density of the monkey skull is . The average denisty of the skull can also be estimated from CT data:12,21 is the porosity of the ’th pixel of the skull, and is the total number of pixels of the skull in the CT image. The estimated average denisty of the skull, is very close (about 1%) to the measured value .
The SOS in the skull was measured using a photoacoustic approach, as shown in Fig. 2(b). A laser beam was split by use of a beam splitter and directed by mirrors to two convex lenses. The two convex lenses (, depth of focus ) focused the laser beam on the inner and outer surface of the skull, and the line connecting the two focused beam spots () was perpendicular to the skull surface. The ultrasonic transducer was placed coaxially with the laser spots; therefore, the average SOS between the two laser spots was calculated byFig. 1(a). The measured SOS values are displayed in the second column of Table 1.
The measured average SOS c¯pa via the PA experiment (column 2) and the estimated average SOS c¯ct from the CT measurements (column 3) for the five measurement locations [see Fig. 1(a)].
|Position||SOS (PA) in m/s||SOS (CT) in m/s|
The corresponding average SOS values were also computed by use of the x-ray CT image data and compared to the measured values. In order to determine the five locations on the CT image that correspond to the measurement locations described above, we measured the arc lengths between the fiducial markers and the measured locations. Then the average SOS at these locations can be estimated from CT data [derived from Eq. (8)]:12,21 and is the resolution of the CT image. The estimated SOS at these locations are shown in the last column of Table 1. The root mean square difference between the SOS inferred from the PA experiment and the SOS inferred from the CT data is .
Chao Huang, Robert W. Schoonover, and Mark A. Anastasio acknowledge support in part by NIH awards EB010049 and EB009715. Liming Nie, Zijian Guo, and Lihong V. Wang acknowledges support by the NIH awards EB010049, CA134539, EB000712, CA136398, and EB008085. L.V.W. has a financial interest in Microphotoacoustics, Inc. and Endra, Inc., which, however, did not support this work.