Open Access
Add to BibSonomyAbstract
Photoacoustic imaging is a hybrid imaging modality capable of producing contrast similar to optical imaging techniques but with increased penetration depth and resolution in turbid media by encoding the information as acoustic waves. In general, it is important to characterize the performance of a photoacoustic imaging system by parameters such as sensitivity, resolution, and contrast. However, system characterization can extend beyond these metrics by implementing advanced analysis via the crosstalk matrix and singular value decomposition. A method was developed to experimentally measure a matrix that represented the imaging operator for a photoacoustic imaging system. Computations to produce the crosstalk matrix were completed to provide insight into the spatially dependent sensitivity and aliasing for the photoacoustic imaging system. Further analysis of the imaging operator was done via singular value decomposition to estimate the capability of the imaging system to reconstruct objects and the inherent sensitivity to those objects. The results provided by singular value decomposition were compared to SVD results from a de-noised imaging operator to estimate the number of measurable singular vectors for the system. These characterization techniques can be broadly applied to any photoacoustic system and, with regards to the studied system, could be used as a basis for improvements to future iterations.
©2010 Optical Society of America
1. Introduction
1.1 Background
Photoacoustic imaging (PAI) is a non-ionizing imaging modality that produces images based on the preferential absorption of optical energy in an absorber by means of the photoacoustic effect. The technique provides images of objects in turbid media with contrast similar to direct optical imaging techniques, but with increased resolution and penetration depth by encoding the optical information as acoustic waves [1,2]. PAI employs the use of a pulsed laser to diffusely irradiate a volume of interest. The optical energy is deposited rapidly allowing the thermal confinement condition to be met, which facilitates the thermo-elastic expansion of the absorbing structure leading to an outwardly propagating transient bipolar pressure wave [3]. Information is contained within the pressure wave regarding the location, size, shape, and optical properties of the absorbing objects [4]. Using the time-domain measurements acquired by acoustic transducers, an image of the distribution of optical absorbers inside the target volume can be inferred using an image reconstruction algorithm [5–10].
Typical metrics utilized to guide imaging system optimization tasks include sensitivity, resolution, and contrast. However, the characterization process can be extended beyond classic metrics by implementing techniques that generate higher level information. This includes the singular value decomposition (SVD) technique and the crosstalk matrix to extract additional system information. The SVD technique produces information concerning the geometry and sensitivity to objects that can be resolved by the system via decomposition of the imaging operator into a set of matrices that are representative of these systems qualities [11,12]. Understanding the geometry of objects that can be sensed by an imaging system is of particular importance when the imaging system acquires a limited number of data projections as system limitations will inevitably be significant when reconstructing images, as is the case with a staring, sparse-array PAI system recently described by our group [13–15]. The crosstalk matrix generates information that describes the spatially dependent sensitivity and aliasing (defined as the inability to distinguish expansion coefficients from each other) of an imaging system. These metrics are important to understanding system resolution in shift-variant imaging systems [16]. In broad terms, the application of the SVD technique and crosstalk matrix to any photoacoustic imaging system can provide a method to comprehensively understand properties of the imaging operator. Consequently, the results can be used to guide improvement and performance by optimization of transducer orientation and bandwidth as well as the number of data projections required to accurately reconstruct objects of relevant geometry.
1.2 Singular value decomposition
Image reconstruction is an inverse problem where the objective is to recover an image from a measured data set. In its discrete form, a noiseless imaging system can be expressed as Eq. (1):
where g is a vector that represents the measured data set, H is an imaging operator, and f is a vector that represents the unknown object(s) that produced the data in g. For our purposes, we assume the point source is band-limited to the dimensions of our expansion functions (the cubic voxel). During image reconstruction, Eq. (1) is solved for f given knowledge of g and H. Ideally, H would be invertible. However, it is generally found that for a real imaging system H is singular. For singular matrices, it can be shown that an M × N matrix, H, can be decomposed by means of Eq. (2):
where U is an M × M matrix, V is an N × N matrix, and both are nonsingular. The M × N matrix S is a diagonal matrix with non-zero diagonal entries representing the singular values of the imaging operator. The decomposition of H into these component matrices is known as the singular value decomposition. The rows of U and columns of V T are the orthonormal singular vectors that provide information regarding the imaging operator. Explicitly, the rows of U and columns of V T form an orthonormal basis for the measurement space and object space, respectively. For example, an object represented by a vector in object space can be projected onto the set of singular vectors in the matrix V T, where the results can be interpreted to indicate the capability with which the system can reproduce an object. It follows that each singular value of S relates the sensitivity of the imaging operator to the corresponding singular vector in V T.
1.3 Crosstalk matrix
The crosstalk concept was introduced by Barrett et al. initially as a way to recover and differentiate among Fourier coefficients used to describe an object in the frequency domain [16,17]. Crosstalk has since been expanded to the analysis of imaging systems in terms of expansion functions in the spatial domain as well as wavelet coefficients [18]. We restrict our analysis to the spatial domain since the expansion functions can be experimentally measured by our imaging system.
As was the case in the SVD analysis, an object can be approximately represented by a series of expansion functions in the spatial domain. In the case of our imaging system, the object is represented by cubic voxels. By constructing the measurement space matrix, the measured data and the voxel coefficients are related by Eq. (1). The crosstalk matrix then allows us to quantify the contribution of each voxel coefficient to the data and is expressed by:
with elements defined by:
where H T represents the transpose of H, j and j’ represent the index of the first and second voxel coefficient, k denotes the product of the time index for a given transducer and the index of the transducer, and K denotes the product of the total number time indices and the total number of transducers. There are two distinct challenges in attempting to recover a voxel coefficient from a discrete set of measurements. First, the voxel coefficient must make a significant contribution to the measurement data, g. Second, the contribution from each voxel coefficient must be distinguishable from contributions made by other voxel coefficients. The crosstalk matrix provides metrics for quantifying a system’s capacity to address these problems. If the constituents of each matrix entry are examined, it becomes apparent that the diagonal elements of the crosstalk matrix determine the sensitivity of the imaging system to the corresponding voxel location in object space. Additional information can be acquired by analyzing the off-diagonal entries in each row of the crosstalk matrix to distinguish among the spatially dependent contributions made by neighboring voxel coefficients, i.e. aliasing of information from other voxels into a voxel of interest. The crosstalk concept can be generalized to any expansion function provided the object can be adequately represented.
1.4 Objective
We have previously reported on a method to calibrate a 3D photoacoustic imaging system by way of the translation of a point source through the imaging volume [15]. This calibration process recorded the response of each transducer in terms of the amplitude, duration, and time-of-arrival of each received pressure wave. Due to technological limitations, the sparseness of the calibration scan was a significant shortcoming in obtaining a comprehensive system response. To more accurately implement the singular value decomposition technique and the crosstalk matrix, it was vital that our calibration process be improved in order to acquire the system response at a step-size approaching the resolution limit of the imaging system. According to Oraevsky et al. [19], the theoretical resolution limit due to transducer bandwidth is approximately 500 μm. However, as shown in previous works, the practical resolution of our system is closer to 1–2 mm [13]. In concurrence with this estimated resolution, it was our objective to acquire a calibration scan with a step-size on the order of the system resolution that spans a volume that approximately corresponds to the transducer sensitivity (25×25×25 mm3). The imaging operator derived from the calibration scan could then be used to comprehensively understand the performance and limitations of our system by SVD and crosstalk matrix analysis.
1.5 Approach
Our approach was to produce a dense calibration scan by means of improving the point source scanning technique published in previous work [15]. While the previous calibration scan was sparse (in relation to the system resolution), the new methodology was utilized to provide a dense system response on the order of the system resolution. Improvements were made to the data acquisition rate and point source translation interval to reduce the net time required to acquire the data from a particular location in object space. After a dense calibration scan was acquired, singular value decomposition and crosstalk analysis were performed on the imaging operator that provided information to characterize our photoacoustic imaging system with regards to the spatially dependent sensitivity, aliasing, as well as the geometry and orientation of objects that can be distinguished. As well, a de-noised imaging operator was constructed in order to compare the experimental imaging operator to a de-noised imaging operator containing the same shift-variant response as the experimental system but with greatly reduced noise.
Methods
2.1 Photoacoustic imaging system
The imaging system utilized 15 ultrasound transducers (model V304, 1” diameter, 2.25 MHz with fractional bandwidth of 65%, Panametrics-NDT, Waltham, Massachusetts) in a staring hemispherical arrangement. Transducers were mounted on 5 custom-built frames, each supporting 3 transducers at zenith angles of 22.5°, 45°, and 67.5°. The frames were designed such that the sensitivity of all 15 transducers were intended to overlap in a specified imaging volume of approximately 25×25×25 mm3 near the geometric center of the array. Laser illumination (“Surelite OPO Plus”, OPO-coupled Nd:YAG, Continuum, Santa Clara, California) was directed to a bifurcated fiber (400 μm diameter) such that half of each laser pulse was guided to a photodiode (to measure pulse-to-pulse variation) and the other half to an optical fiber immersed in the liquid (where the photoacoustic signal was generated) for a total of 16 channels collecting data (15 transducers, 1 photodiode). The pulse duration was 6 ns at a repetition rate of 10 Hz with a maximum laser output of approximately 100 mJ/pulse. Note that only a small fraction of the pulse was accepted by the fiber due to its small core size relative to the beam diameter (~1.5 cm). All calibration scans were done at 675 nm. Each transducer was electrically connected to a dedicated channel on a preamplifier card (custom built). The analog signals were acquired in parallel, converted to digital signals, and sent to a personal computer for analysis. The custom built data acquisition system sampled with 14-bit resolution at a frequency of 50 MHz. The PA system (with PA point source and optical fiber) is shown in Fig. 1(a) while a representative PA time series acquired during an experiment from a single transducer is shown in Fig. 1(b).
Fig. 1. (a) Isometric view of the hemispherical PA imaging array illustrating the transducer arrangement, placement of the liquid reservoir, and the optical fiber PA source. (b) Example of raw data acquired on a single acoustic transducer.
2.2 System calibration scan
In order to acquire calibration scans with step-size on the order of the system resolution, a number of improvements were made to the scan procedure in comparison to our previous work [15]. First, the linear slides responsible for translating the optical fiber through the absorbing liquid were replaced with a SCARA robot (Model E2C351S - UL, Epson), which could translate the source between data points faster and more reliably to reduce translation time. Second, the data acquisition cards were updated with USB 2.0 connections to relay data to the PC in a much shorter interval. With the upgrades, each scan point required approximately 15 seconds to complete. Subsequently, two calibration scans were performed with different imaging volumes and step-sizes. The first scan was completed with a 16–16–16 mm3 imaging volume and 2 mm step-size for a total of 512 data points. The second scan was performed on an imaging volume of 30×30×mm3 and 3 mm step-size for a total of 1000 data points. At each test position in the scan, the PA signal was averaged over 10 pulses and recorded simultaneously on all 15 transducers. After the calibration scan, the time series data for each transducer and grid location was analyzed off-line to obtain the imaging operator corresponding to each scan. Analysis included extracting the time series of a particular transducer and grid point, rectifying the time series, and then smoothing the time series using a moving average with a bin size of 40 points. Each time series was copied to the matrix representative of the imaging operator. Each row contained the concatenated time series for all 15 transducers corresponding to a position in object space. Therefore, the imaging operator had rows corresponding to the number of calibration grid points in object space and columns corresponding to the number of time points used to sample object space multiplied by the number of transducers used to collect the data. A de-noised imaging operator was constructed by removing the noise from the transducer responses prior to analysis. De-noising was performed by estimating the peak size, peak width and time of flight from the smoothed time series data. The peak size was found by locating the maximum in the time series data and the peak width estimated by differencing the location of the time points, which corresponded to half the value of the peak size. The time of flight was recorded as the temporal location of the peak. The parameters were then used to compute a synthetic time series consisting of zeros everywhere, except for the points representing the peak width centered upon the time of flight. These specific points were filled with a scaled and inverted parabola representative of the smoothed experimental time series. The inverted parabola is representative of the basis function used in the back projection model of the reconstruction algorithm from our previous work, which closely resembles the velocity potential of the bipolar pressure signal. The de-noised imaging operator was then constructed from these de-noised time series using the same method described above. The de-noised imaging operator then contained the same shift-variant response as the experimental imaging operator but with greatly reduced noise.
2.3 Singular value decomposition and singular vector correlation
Singular value decomposition of both imaging operators was performed in MATLAB via the built-in singular value decomposition function (svds, MATLAB version 7.8.0). The orientation (positive or negative) of the resulting singular vectors is not necessarily the same between imaging operators and is relatively unimportant when interpreting the physical meaning of the singular vectors. The inner product of the experimental singular vectors and de-noised singular vectors in matrix V T was computed after the SVD of both imaging operators. Because singular vectors parallel and orthogonal to each other were expected to result in a value of one and zero, respectively, we interpreted the inner product result as a correlation between the two singular vectors.
2.4 Crosstalk matrix
The experimental voxel crosstalk matrix was computed by multiplying the transpose of the experimental imaging operator by the experimental imaging operator.
3. Results
3.1 Crosstalk sensitivity and aliasing
After the crosstalk matrix was calculated for the large volume scan, the main-diagonal was reshaped to represent the location of each voxel in object space (to facilitate straightforward visualization of the data) and was plotted in Fig. 2. The map in Fig. 2 visually illustrates the sensitivity of the transducer arrangement to signals originating from each voxel location in object space for the 30×30×30 mm3 scan (10×10×10 voxels).
Fig. 2. Displays sensitivity of the PA system at each location in object space acquired from the main-diagonal of the crosstalk matrix corresponding to the 30×30×30 mm3. Both x and y axes represent voxel number in the y and z directions, respectively. Accordingly, each x-plane in object space is 10×10 voxels.
Figures 3(a) and 3(b) represent aliasing of signal originating at the center of the object space to all other voxel locations for the small and large scan, respectively. The voxel highlighted in Figs. 3(a) and 3(b) corresponds to the same location in object space (i.e. center point of object space). Figure 3(c) shows aliasing effects from a selected corner voxel to illustrate the shift-variant response of the PA system. The y-z plane in Figs. 3(b) and 3(c) is 30×30 mm2 and is 16× 16 mm2 in Fig. 3(a). The aliasing information was retrieved by reshaping the rows of the crosstalk matrix (in the same manner as described for Fig. 2). For example, aliasing in the center voxel is visualized by plotting row 455 of the crosstalk matrix for the large scan.
Fig. 3. (a) Illustrates aliasing from the center voxel for the 16×16×16 mm3 scan (each x-plane is 8×8 voxels) while (b) shows aliasing from the same position for the 30×30×30 mm3 scan (each x-plane is 10×10 voxels). (c) Shows representative aliasing plots from a voxel located at the corner of the imaging volume for the 30×30×30 mm3 scan (each x-plane is 10×10 voxels).
3.2 Singular value decomposition: Singular vectors
The decomposition of the imaging operator, according to Eq. (2), yields a set of orthonormal singular vectors that describes both the projection and object space of the PA system. An additional imaging operator was produced to better understand the results of the SVD. A de-noised imaging operator was produced by modifying the experimental imaging operator according to the strategy outlined in section 2.2. The signals acquired at each voxel in object space were replaced with a parabola of the same peak, width, and time of flight as the experimental signals as well as zeros in all other temporal recordings in order to examine the effects of noise on the SVD of the imaging operator. Column vectors of matrix V T were organized in the same manner as the data in Fig. 2 in order to aid in visualization of the results. The data corresponding to the experimental imaging operator for the 30×30×30 mm3 scan is displayed in Fig. 4(a). The de-noised imaging operator derived from the experimental imaging operator is shown in Fig. 4(b). The center plane of object space was then plotted in both the de-noised and experimental decompositions.
The normalized singular vectors of the de-noised and experimental imaging operators were multiplied and summed in order to correlate the similarity of each corresponding set of singular vectors. The absolute value of the product is presented as the correlation. The results of the correlation are shown in Fig. 5(a). Four additional imaging operators were constructed with systematic noise added to the de-noised imaging operator at values of ½, μ, 2, and 5 times the average system noise of the experimental imaging operator. The singular vectors of each imaging operator were again computed and projected onto the de-noised imaging operator values and resulted in data similar (but not shown) to that displayed in Fig. 5(a). In an ideal system, each set of singular vectors would correlate with a value of 1 while the multiplication of a singular vector with another singular vector in the basis would yield a value of 0 as the singular vectors would be orthogonal. The data in Fig. 5(a), as well as the correlation data not shown for the four additional imaging operators, were reordered and displayed in Fig. 5(b) to illustrate the correlation of singular vectors in descending order.
Fig. 4. (a) and (b) Displays the center y-z plane of the first 8 singular vectors acquired via experiment and de-noised, respectively. The field-of-view for each singular vector is 30×30 mm2. The singular vector number reads from left to right with the leftmost image representing singular vector 1.
Fig. 5. (a) Displays the correlation among the set of 1000 corresponding singular vectors in the de-noised and experimental matrices, (V)T. (b) Shows the same computation as in (a) but in descending order for an imaging operator with (i) ½ the intrinsic system noise, (ii) ¼ the intrinsic system noise, (iii) the experimental imaging operator, (iv) 2 times the intrinsic system noise, and (v) 5 times the intrinsic noise. The singular vector index changed with the order as the true singular vector number (corresponding to the matrix (V)T) was unknown after the initial projection operations were completed. The vertical axis in both (a) and (b) is shared.
4. Discussion
4.1 Crosstalk sensitivity and aliasing
As was demonstrated in an earlier publication [15], the sensitivity of the transducer array lacked uniformity through the imaging volume specified by the field-of-view of the transducer array. Results presented in the previous publication illustrated only individual transducer sensitivities, which provided only minor insight into the overall system sensitivity (cumulative sensitivity) to object space when interpreted together. Furthermore, in the previous work, the individual transducer sensitivities were incorporated into the image reconstruction process to correct for potential non-uniformity in image contrast that might arise due to non-uniform transducer sensitivity to object space. Here, the cumulative sensitivity was obtained directly by means of the main diagonal elements of the crosstalk matrix (Fig. 2) and confirmed the non-uniform sensitivity of the transducer array to object space. Although the cumulative sensitivity is not used for reconstruction purposes, it provides insight into the overall performance of the transducer array and can be used for optimization purposes (see below).
Figures 3(a) and 3(b) illustrate the aliasing from a similar location in object space. Because the step-size and volume of the defined object space is significantly different, the span of the aliasing effects appear more pronounced in Fig. 3(a). However, the span of aliasing measured as the FWHM of the highlighted voxel is similar in both maps. The span in x, y, and z is approximately 2/2 mm, 2/1 mm, and 2/2 mm, for Figs. 3(a), 3(b), respectively. In general, the aliasing contributions in object space are not isotropic because the PA system is shift-variant, which is evident when comparing aliasing from a centre voxel [Fig. 3(b)] to a corner voxel [Fig. 3(c)]. If we recall the process for calculating the crosstalk matrix, there are two independent contributions that lead to aliasing in object space. The signal from one voxel into another is influenced by both the extent of the overlapping temporal region among the two voxels as well as the shape of the expansion functions generated by the two voxels in the overlapping temporal region. For example, measured data generated from voxels in object space may possess a relatively small temporal overlap but could still contain large aliasing effects if the shape of the expansion functions in the overlapping region contains a relatively large area under the curve. This would result in a high crosstalk among voxels with relatively small temporal overlap. In our previous work, the backward model used during image reconstruction to compute the distribution of sources in object space accounted for the non-uniformity in transducer sensitivity, but assumed a symmetric expansion function. Accordingly, the aliasing effects expected in reconstructed images will not precisely correspond to that displayed in the crosstalk matrix.
Because of the shift-variant nature of the PA system, it is inaccurate to obtain global estimates of resolution and contrast via the crosstalk calculation. However, general behaviors of the system performance can be visualized and qualitative characterizations can be made. Qualitative metrics computed from the crosstalk matrix could potentially be used to change the placement of transducers in the array to reduce aliasing and enhance system performance. The aliasing effects computed via the crosstalk matrix are inherent to this PA system based on the frequency response of the transducers used as well as the relative position of the object space to the transducer arrangement. Opting for transducers of higher centre frequency would limit the aliasing effects seen but would fundamentally change the system performance in other, negative ways, such as reducing the penetration depth of the acoustic waves. To reduce aliasing effects to a practical limit, a greater number of transducers need be introduced to diminish the consequences of aliased signal.
4.2 Singular value decomposition
Recall the method for producing the de-noised imaging operator. Only the system noise was reduced after acquiring the imaging operator. It is clear that a strong correlation should exist between singular vectors of the same order when comparing the de-noised and experimental imaging operator provided the system noise has a minor effect on the imaging operator. In the case where the noise has a small effect on the geometry of the singular vector, it is expected that the correlation between the two singular vectors approaches one. However, when the system noise is significant, the correlation among singular vectors will be reduced considerably. In a mathematical context, the correlation between differing singular vectors should be zero because the basis of decomposed singular vectors is orthonormal. However, this relationship is not seen in practice due to system noise. The SVD of the imaging operator presents the singular vectors in ascending order of the corresponding singular value in matrix S. Therefore, there is potential for the system noise in the experimental imaging operator to impact the order of the decomposed singular vectors such that comparison of singular vectors of the same order will not necessarily represent a comparison of singular vectors of corresponding geometry. At some threshold, the impact of system noise is too significant for the system to accurately resolve the geometry of the singular vector and consequently the affected singular vector contributes no useful information when attempting to recover an object.
The correlation of the experimental singular vectors to the singular vectors obtained from the de-noised imaging operator is shown in Fig. 5(a). As expected, the lower order singular vector pairs had a high correlation (as high as 0.98 for the first pair of singular vectors), but the correlation dropped quickly as the order of the singular vector increased. Although the purpose of the plot was to gain insight into the number of measurable singular vectors, selection of a threshold correlation that delineated the number of measureable singular vectors was not obvious. For example, the distribution of correlation values appeared to have at least two components with points of inflection at indices 100 and 600. To better understand the sensitivity of the correlation distribution to measureable singular vectors, we compared the experimental findings to correlation values obtained from pairs of singular vectors, where defined amounts of noise were introduced into the imaging operator (i.e. de-noised imaging operator and the de-noised imaging operator with noise added back). The correlation distributions were reordered in descending order and displayed in Fig. 5(b) to facilitate interpretation. Curves (i) and (ii) in Fig. 5(b) were derived from imaging operators with less system noise (1/4 and 1/2, respectively) than the experimental imaging operator and consequently had a broader correlation distribution. Curves (iv) and (v) in Fig. 5(b), generated with 2 and 5 times the experimental noise, respectively, displayed a narrow correlation distribution. Taken together, the curves suggested that a correlation of 0.2 represented a reasonable threshold for delineating the number of measureable singular vectors due to the presence of a clearly defined inflection point below this correlation value for curves (i, ii, iv, and v). Therefore, we concluded that approximately 400 measurable singular vectors were present for the experimental imaging operator [using a 0.2 correlation threshold in Fig. 5(a)].
Although it is generally the goal of any imaging system to resolve as many singular vectors as possible, it should be emphasized that not all imaging tasks necessarily require a large number of measureable singular vectors to resolve objects in the field of view. This is illustrated by considering two canonical examples. (i) If the PA system is to macroscopically localize a spherical tumor mass in soft tissue, it may suffice to have a system that resolves a small number of singular vectors (perhaps several hundred) since the task is to reconstruct a single object of low morphological complexity within the field of view of the imaging system. (ii) A more complicated PA imaging task such as delineating microvasculature within a small animal may require many measurable singular vectors (perhaps several thousand) since many objects of complex morphology will be present in the field of view.
4.3 Computation considerations
It is important to note a practical shortcoming concerning the computation of both the crosstalk matrix and singular value decomposition. That is, the computation required to compute the associated matrices can be lengthy. The computation of the crosstalk matrix using the imaging operators was not intensive (because of the relatively small temporal domain, voxel number, and transducer count) while the computation of the SVD matrices with 1000 singular values required 53 minutes (Dell T7400 workstation: Dual Intel® Xeon® X5472 3.00 GHz, 8 GB Ram, Windows Vista-64). However, this could be a limiting factor for other PA systems that utilize a larger temporal domain, voxel number, and transducer count. It should be stressed that these computations need only be completed once to extract the required information.
4.4 Imaging considerations
In order to determine whether or not an imaging task fails or succeeds, an object described by the cubic voxel expansion function should be projected onto the set of computed singular vectors in the matrix V T. The fidelity with which the image is produced indicates the capability of the imaging system to fundamentally capture the information contained in the object. However, whether or not the task succeeds or fails is necessarily subjective based on the required task of the imaging system.
The link between crosstalk and SVD is important. While it is the SVD matrices that are directly responsible for producing an image from a measured set of data, it is the crosstalk matrix that provides insight in addressing why an imaging task failed or succeeded. It is the expansion coefficients describing an object that must be resolved in order to accurately capture the information when performing an imaging task. The capability of the system to do this is intrinsically contained within the SVD matrices; however, it is the crosstalk matrix that actually quantifies the spatially-dependent overlap of expansion coefficients (aliasing). If an imaging system subjectively fails an imaging task, the aliasing contribution extracted from the crosstalk matrix can be assessed and the imaging system can be improved where necessary.
5. Conclusion
A technique was developed to acquire a data set that described an imaging operator for a PAI system. Two scans were completed. The first contained a step size of 2 mm within an object space of 16×16×16 mm3. The second was completed at a step size of 3 mm within an object space of 30×30×30 mm3. Utilizing this data, computations to produce a voxel-based crosstalk matrix were made in order to characterize the spatially dependent sensitivity and aliasing. The lack of uniformity in the sensitivity confirmed the findings of our previous work. Singular value decomposition analysis was performed on the imaging operator to provide insight into the system’s sensitivity to objects of complex geometry. As well, insight was provided regarding the sensitivity of the imaging operator to noise. Ultimately, both techniques provided information that could be used to understand any PA system and could provide a means to improve future iterations of the imaging hardware.
Acknowledgements
The authors would like to acknowledge Lynn Keenliside and John Patrick for their valuable technical help. Michael Roumeliotis is supported by the London Health Sciences Centre through the Translational Breast Cancer Research Unit (TBCRU), the Canadian Institute for Health Research (CIHR) and by the University of Western Ontario (UWO). Research funding was provided by the Ontario Research Fund (ORF), the Canadian Foundation for Innovation (CFI), Natural Sciences and Engineering Research Council (NSERC), and MultiMagnetics Inc. EN was supported by a NSERC USRA. MAA and GC were supported in part by award NIH R01EB010049.
References and links
1. M. Xu and L. V. Wang, “Photoacoustic imaging in biomedicine,” Rev. Sci. Instrum. 77(4), 041101 (2006). [CrossRef]
2. T. Lu, J. Jiang, Y. Su, R. K. Wang, F. Zhang, and J. Yao, Photoacoustic imaging: Its current status and future development” in 4th International Conference on Photonics and Imaging in Biology and Medicine, September 03,2005 – September 06 (SPIE), National Natural Science Foundation of China; SPIE Russia Chapter; Int. Laser Center of M.V. Lomoson Moscow State Univ.; Bio-optics and Laser Medicine Comm. of Chinese Optics Soc.; Science and Techn. Garden of Tianjin University, China.
3. G. J. Diebold, T. Sun, and M. I. Khan, “Photoacoustic monopole radiation in one, two, and three dimensions,” Phys. Rev. Lett. 67(24), 3384–3387 (1991). [CrossRef] [PubMed]
4. 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). [CrossRef] [PubMed]
5. P. Liu, “The P-transform and photoacoustic image reconstruction,” Phys. Med. Biol. 43(3), 667–674 (1998). [CrossRef] [PubMed]
6. C. G. A. Hoelen and F. F. de Mul, “Image reconstruction for photoacoustic scanning of tissue structures,” Appl. Opt. 39(31), 5872–5883 (2000). [CrossRef]
7. D. Frauchiger, K. P. Kostli, G. Paltauf, M. Frenz, and H. P. Weber, Optoacoustic tomography using a two dimensional optical pressure transducer and two different reconstruction algorithms” in Hybrid and Novel Imaging and New Optical Instrumentation for Biomedical Applications, June 18,2001 – June 21 (SPIE), 74–80.
8. M. Xu and L. V. Wang, “RF-induced thermoacoustic tomography” in Proceedings of the 2002 IEEE Engineering in Medicine and Biology 24th Annual Conference and the 2002 Fall Meeting of the Biomedical Engineering Society (BMES / EMBS), October 23,2002 – October 26 (Institute of Electrical and Electronics Engineers Inc), 1211–1212.
9. K. P. Kostli, D. Frauchiger, J. J. Niederhauser, G. Paltauf, H. P. Weber, and M. Frenz, “Optoacoustic imaging using a three-dimensional reconstruction algorithm,” IEEE J. Sel. Top. Quantum Electron. 7(6), 918–923 (2001). [CrossRef]
10. P. Ephrat, L. Keenliside, A. Seabrook, F. S. Prato, and J. J. L. Carson, “Three-dimensional photoacoustic imaging by sparse-array detection and iterative image reconstruction,” J. Biomed. Opt. 13(5), 054052 (2008). [CrossRef] [PubMed]
11. D. Modgil, M. A. Anastasio, and P. J. La Riviere, Photoacoustic image reconstruction in an attenuating medium using singular value decomposition” in Photons Plus Ultrasound: Imaging and Sensing 2009 (SPIE - The International Society for Optical Engineering), 71771B (7 pp.).
12. D. W. Wilson and H. H. Barrett, “Decomposition of images and objects into measurement and null components,” Opt. Express 2(6), 254–260 (1998). [CrossRef] [PubMed]
13. P. Ephrat and J. J. L. Carson, Measurement of photoacoustic detector sensitivity distribution by robotic source placement” in 9th Conference on Photons Plus Ultrasound: Imaging and Sensing 2008, January 20,2008 – January 23 (SPIE), Society of Photo-Optical Instrumentation Engineers (SPIE).
14. P. Ephrat, M. Roumeliotis, F. S. Prato, and J. J. Carson, “Four-dimensional photoacoustic imaging of moving targets,” Opt. Express 16(26), 21570–21581 (2008). [CrossRef] [PubMed]
15. M. Roumeliotis, P. Ephrat, J. Patrick, and J. J. L. Carson, “Development and characterization of an omnidirectional photoacoustic point source for calibration of a staring 3D photoacoustic imaging system,” Opt. Express 17(17), 15228–15238 (2009). [CrossRef] [PubMed]
16. H. H. Barrett, J. L. Denny, R. F. Wagner, and K. J. Myers, “Objective assessment of image quality. II. Fisher information, Fourier crosstalk, and figures of merit for task performance,” J. Opt. Soc. Am. A 12(5), 834–852 (1995). [CrossRef]
17. H. H. Barrett and H. Gifford, Cone-beam tomography with discrete data sets” in Second International Meeting on Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine, 451–76.
18. J. Qi and R. H. Huesman, “Wavelet crosstalk matrix and its application to assessment of shift-variant imaging systems,” IEEE Trans. Nucl. Sci. 51(1), 123–129 (2004). [CrossRef]
19. A. A. Oraevsky, V. G. Andreev, A. A. Karabutov, and R. O. Esenaliev, Two-dimensional opto-acoustic tomography transducer array and image reconstruction algorithm,” Proc SPIE Int Soc Opt Eng 3601, 256–267 (1999).
