Open Access
Add to BibSonomyAbstract
Photoacoustic imaging is a non-ionizing imaging modality that provides contrast consistent with optical imaging techniques while the resolution and penetration depth is similar to ultrasound techniques. In a previous publication [Opt. Express 18, 11406 (2010)], a technique was introduced to experimentally acquire the imaging operator for a photoacoustic imaging system. While this was an important foundation for future work, we have recently improved the experimental procedure allowing for a more densely populated imaging operator to be acquired. Subsets of the imaging operator were produced by varying the transducer count as well as the measurement space temporal sampling rate. Examination of the matrix rank and the effect of contributing object space singular vectors to image reconstruction were performed. For a PAI system collecting only limited data projections, matrix rank increased linearly with transducer count and measurement space temporal sampling rate. Image reconstruction using a regularized pseudoinverse of the imaging operator was performed on photoacoustic signals from a point source, line source, and an array of point sources derived from the imaging operator. As expected, image quality increased for each object with increasing transducer count and measurement space temporal sampling rate. Using the same approach, but on experimentally sampled photoacoustic signals from a moving point-like source, acquisition, data transfer, reconstruction and image display took 1.4 s using one laser pulse per 3D frame. With relatively simple hardware improvements to data transfer and computation speed, our current imaging results imply that acquisition and display of 3D photoacoustic images at laser repetition rates of 10Hz is easily achieved.
©2011 Optical Society of America
1. Introduction
1.1 Background
Photoacoustic imaging (PAI) is a dual imaging modality that makes use of the contrast inherent to optical imaging techniques as well as the penetration depth and resolution of ultrasound techniques [1]. This is accomplished through illumination of optically absorbing objects with non-ionizing radiation generated from a pulsed laser. The optical energy is deposited rapidly, satisfying thermal and stress confinement conditions, which facilitates thermo-elastic expansion of the absorbing structure leading to an outwardly propagating transient bipolar pressure wave [2]. Characteristic information of the absorbing structures, including location, size, shape, and optical properties, is encoded in the propagating pressure waves [3]. Using time-domain measurements, photoacoustic images can be produced using a variety of reconstruction algorithms [4–9].
In an earlier publication, we introduced a technique that experimentally measured a matrix that represented the imaging operator for a photoacoustic imaging system [10]. Metrics intended to aid in the analysis of the imaging operator were then computed to provide information related to system performance. Among the techniques used was singular value decomposition (SVD). In a general context, the SVD technique provides information regarding the geometry of objects that can be accurately captured by an arrangement of transducers from a specified object space. This information is especially important when collecting a small number of data projections since the limited data will inevitably reduce the complexity of objects that can be accurately retrieved. This is the case with the staring, sparse-array PAI system described by our group [11–13].
While the previous work was important in introducing the process of experimentally acquiring the imaging operator, it did not address detail with respect to the number of object space singular vectors that can be reliably recovered and the impact these singular vectors directly have on image quality. This is vitally important as these singular vectors strictly determine the complexity of objects that can fundamentally be retrieved by a reconstruction algorithm, irrespective of the specific merits of the reconstruction technique. As well, the acquisition of an experimental imaging operator enables reconstruction of PA images to be computed directly by solving a linear system model, which translates time-dependent pressure measurements to an estimate of the object. In this paper, we examine the consequence of both variable transducer count and measurement space temporal sampling rate on the capability of the PA imaging system to reconstruct arbitrary objects in the object space. As well, a specific arrangement of transducers and measurement space temporal sampling rate is utilized to perform real-time imaging of a point source as it is translated through object space. In order to produce accurate and fast PA images, it is critical the effect of both transducer count and temporal sampling rate are understood in order to capture an imaging operator with only the precisely required information. It is these parameters that directly determine the number of elements in the imaging operator. As a consequence, improper sampling of the object space can lead not only to infidelity in image reconstruction but also increased computation speed.
1.2 Singular value decomposition
Image reconstruction is an inverse problem, where the objective is to recover an estimate of the object that produced a measured data set. A linear model is typically used to describe the imaging system to first order. This is expressed mathematically as
where g is a vector that represents the measured data set, H is the imaging operator, and f is a vector that represents a finite-dimensional approximation of the unknown object(s) that produced the data in g. It is generally found that H is singular and cannot be inverted directly. For singular matrices, it can be shown that an M x N matrix, H, can be decomposed by means of Eq. (2):where U is an M x M matrix, V is an N x N matrix, and both are nonsingular. The M x N matrix S is a diagonal matrix with non-zero diagonal entries representing the singular values of the imaging operator. The rows of U and columns of V T are the orthonormal singular vectors that form a complete linear vector space (in their respective dimensions) describing the measurement space and discretized object space, respectively. It follows that each singular value of S relates the sensitivity of the imaging operator to the corresponding singular vectors in U and V T.1.3 Estimate of effective singular values
Upon decomposing the imaging operator, the vectors provided in VT are linearly independent. However, by examining the associated magnitude of the singular values in matrix S, it is clear not all vectors contribute equally to the overall system response. In fact, some do not effectively contribute at all to the reconstruction of an object [14]. It is the matrix rank (number of linearly independent rows) of the imaging operator that indicates the singular vectors that contribute usefully to image reconstruction. A number of techniques have been proposed to determine the rank of a matrix in the context of a real imaging operator [15–17]. Konstantinides et al. used a statistical model based on the signal-to-noise ratio of a measured imaging operator. The technique proved to be very accurate in determining the rank of the matrix.
As referenced in [15], a single criterion can be used to evaluate the effective rank of the imaging operator. Mathematically,
where α represents the magnitude of the singular value and t is the effective rank of the matrix. While this threshold is seemingly subjective, it can be evaluated more precisely by stating the rank of the matrix is estimated to be equal to t if,Via the results of [15–17], it is clear that in an experimental imaging operator the singular values in S become exceedingly small when the remaining columns of VT do not contribute useful information about object space.An estimate of the effective number of singular vectors that contribute usefully to image reconstruction is critical when the image reconstruction is performed by solving directly for the vector f in Eq. (1). In this paper, the SVD of the imaging operator, H, was computed and utilized to compute a regularized pseudoinverse solution to Eq. (1). The regularized solution utilized only the singular values and associated singular vectors that were identified as corresponding to useful information as described above.
1.4 Experimental objective
We previously reported on a method to acquire a matrix that represented the imaging operator of a 3D photoacoustic imaging system [10,13]. At each point in object space, the amplitude, duration, and time of flight of the pressure wave was recorded for each transducer. This provided the foundation of the imaging operator, which described the system response to a point source for each transducer at each location in object space. Due to technological limitations, the previous calibration scan did not optimally cover the entire object space at a step-size on the order of the system resolution. It was necessary to improve the data acquisition rate in order to facilitate a calibration scan with many more points in object space. In our present study, we have augmented the PA system to now include a total of 30 transducers.
Once the dense calibration scan was acquired, subsets of the transducer arrangement could be generated. The imaging operators contained the full set of 30 transducers or subsets with 25, 20, 15, 10 or 5 transducers. Similar subsets of the experimental imaging operator were created with different measurement space temporal sampling rates with a constant count of 15 transducers. The SVD of each imaging operator was subsequently computed. Different imaging operators were compiled from the same calibration scan in order to eliminate variability between scans. This approach had the added benefit that any variation in matrix rank was due solely to the arrangement, number of transducers, or sensitivity of the transducers contributing to each imaging operator.
2. Methods
2.1 Photoacoustic imaging system
The imaging system was consistent with the experimental setup referred to in our earlier publication [10]. Data acquisition was accomplished with 2 custom-built 16-channel preamplifier cards electronically connected to 4 custom-built 8-channel data acquisition cards held within the same chassis. A total of 32 channels collected data (30 transducers, 1 photodiode, 1 channel was unused). Data acquisition, transfer to the PC, image reconstruction and display was controlled with a custom application developed within National Instruments LabviewTM v8.5. The imaging system utilized a total of 30 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 custom frames, each supporting 3 transducers at evenly spaced azimuthal angles (36° between columns). Fifteen (15) of the transducers were mounted on frames with zenith angles of 22.5°, 45°, and 67.5°. The remaining 15 transducers were mounted on frames placed azimuthally between the original frames, at zenith angles of 33.75°, 56.25°, and 78.75°. A schematic of the transducer array is shown in Fig. 1(a) .
Fig. 1 (a) Isometric view of the hemispherical PA imaging array illustrating the transducer arrangement. Columns with transducers lightly shaded in green correspond to zenith angles of 22.5°, 45°, and 67.5° while columns with transducers lightly shaded blue correspond to zenith angles of 33.75°, 56.25°, and 78.75°. (b) Represents an unfolded schematic of (a) whereby each plane of 5 transducers is referenced P1 through P6, with P1 representing the bottom-most row and P6 the top-most row. (c) Summarizes the results of the estimate of effective matrix rank on the 6 different transducer arrangements. (d) Semi-logarithmic plot of the magnitude of the singular values versus singular value index for each of the 6 transducer arrangements. Left-most curve corresponds to the 5 transducer imaging operator while right-most to the 30 transducer imaging operator with intermediate curves represent increasing transducer count from left to right.
2.2 System calibration scan
As in section 2.1, the calibration arrangement was developed in an earlier publication and is cited here for reference [10]. The calibration scan was completed with a 30x30x30 mm3 object space and 1.5 mm step-size for a total of 8,000 voxel locations. It was originally intended to complete scans of varying voxel count and step-size. However, calibration scans exceeding 8,000 voxels represented the practical limit of our experimental methodology. For the experimental image reconstruction, a second scan of dimension 16x16x7 mm3 was produced with 30 transducers at a measurement space temporal sampling rate of 10 MHz.
2.3 Singular value decomposition
The processed photoacoustic time series were compiled to produce a matrix representative of the imaging operator. For a given voxel, the time series recorded by each transducer were concatenated to produce one row in the matrix. Therefore, a single imaging operator contained rows equal to the number of voxels in the scan (8,000) and columns equal to the number of transducers multiplied by the number of time points used to record the photoacoustic data.
Singular value decomposition of both imaging operators was performed in MATLAB® (The Mathworks, Inc., version 7.8.0, Natick, Massachusetts) via the built-in singular value decomposition function (svds). A technical limitation was reached with large imaging operators. Specifically, as the size of the matrix approached 10GB, decomposition by SVD in MATLAB failed due to the available memory in the workstation. As a work-around, we employed a 10-point average to reduce the dimensions of the measurement space by 10-fold. Imaging operators acquired at this temporal sampling rate (5 MHz) were used to compare the change in matrix rank as the transducer count was varied. Afterward, one particular arrangement of 15 transducers was selected and an n-point average was applied to examine the effects of measurement space temporal sampling rate on the matrix rank.
2.4 Regularized pseudoinverse image reconstruction
Image reconstruction was performed by solving directly for the vector f in Eq. (1). The SVD component matrices of Eq. (2) were used to represent the imaging operator in solving for the vector f. The built-in MATLAB function (pinv) was used to compute the pseudoinverse of the matrix S, referred to as S +. Due to noise present in the imaging operator, singular values with indices greater than 90 percent of the computed matrix rank were set to zero when computing S+. The pseudoinverse was then multiplied by the data set measured experimentally in g, yielding an estimate of the object f.
A second reconstruction technique was used for basis of comparison on a single system setup (30 transducers, 5 MHz). This technique was an iterative algebraic reconstruction technique, which estimated the image by computing the difference between the data set, g, and the right side of Eq. (1), which was then added to a master of the image. While this model is computationally expensive, it was deemed instructive to show the results of the simulated data sets being reconstructed via another approach.
2.5 Real-time photoacoustic imaging
The robot used to perform the calibration scan was also used to translate the optical fiber (photoacoustic point-like source) through the object space of the PA imaging system at a speed of 0.40 mm/s during an imaging experiment. One image of the point-like source was captured every 1.4 s. The interval time represented a practical limit related to the laser pulse duration (9 ns), the data acquisition time (100 µs), the time taken to trigger and transfer 32 channels of data from the 4 DAQ cards to the PC (~1.2 s), and the time taken to perform the matrix multiplication, archive the data to disk, and display the data to PC screen (~0.2 s). Three slices of the 3D image data (xy, xz and yz) were displayed on the workstation after every image reconstruction and before the next laser pulse. Repeated laser pulses resulted in a 3D photoacoustic movie (i.e. real-time 4D photoacoustic imaging).
3. Results
3.1 Estimate of matrix rank
The PA system is shown in Fig. 1(a) while an unfolded representation of the transducer arrangement shown in Fig. 1(b). The schematic shown in Fig. 1(b) was used to describe the varying transducer arrays for which matrix rank was computed. Each plane (P1 through P6) of Fig. 1(b) represents the transducers at the same zenith elevation. The results of the matrix rank estimate (based on Eq. (4)) are summarized in Fig. 1(c). The raw data - a semi-logarithmic plot of the magnitude of the singular values versus singular value index for each of the transducer arrangements - is shown in Fig. 1(d). It is evident in Fig. 1(d) that an abrupt decline in singular value magnitude is present in each of the six curves.
The results shown in Figs. 1(c) and 1(d) are summarized and plotted in Fig. 2(a) . An analysis similar to that shown in Fig. 1(d) was performed to estimate the matrix rank on imaging operators with variable measurement space temporal sampling rate (raw data not shown). Imaging operators were produced with measurement space temporal sampling rates of 25, 12.5, 10, 5, 2.5, 1.25, 1 and 0.5 MHz, all with the 15 transducer arrangement according to the setup described in Fig. 1(c). The matrix rank was then estimated for each of the 8 temporal sampling rates and plotted in Fig. 2(b). Three (3) distinct regions are outlined in Fig. 2(b) by the dotted lines. Area (i) is the region where the system is undersampled (according to the Nyquist criterion and an estimate of the system resolution). Area (ii) corresponds to the sampling regime that meets the Nyquist criterion and where the object space dimension is greater than the measurement space dimension. Finally, area (iii) describes a region where the Nyquist criterion is met and the measurement space dimension exceeds the object space dimension. To further clarify the effective matrix rank estimate in Fig. 2(b), linear regression analysis was performed on the data contained within region (ii). A clear linear trend is shown in both plots with both coefficients of determination greater than or equal to 0.999. While this trend was expected, it is confirmed here experimentally for a photoacoustic imaging system with low transducer count. To emphasize the non-linear behavior of the matrix rank in region (i) and (iii), Fig. 2(c) shows a plot of the residual errors for each temporal sampling rate. Residual error was computed by comparing the expected matrix rank (based on the slope in Fig. 2(b)) to the measured matrix rank (estimated using Eq. (4)). Finally, sample singular vectors are plotted in Fig. 2(d) for the imaging operator of 30 transducers and 5 MHz temporal sampling rate. Each image depicts the centre plane of singular vector 1, 10, 3632, and 4036 for image (i) through (iv), respectively. Images (i) and (ii) were selected to show discernable geometry in the singular vectors at low orders, while (iii) corresponds to the singular vector at 90% of the matrix rank, and (iv) shows the vector at one index beyond the matrix rank.
Fig. 2 (a) Displays the estimated matrix rank for variable transducer count and arrangements. (b) Displays the estimated matrix rank for variable measurement space temporal sampling rates. Linear regression for (b) was performed only on the 4 data points contained within region (ii). The line is shown throughout the entire figure to show the expected value of the matrix rank. (c) Compares the expected matrix rank to the measured rank and is plotted as a percent error to highlight the deviation from linearity in regions (i) and (iii). (d) Provides a visual interpretation of the geometry associated with selected singular vectors for the imaging operator corresponding to the 30 transducer, 5 MHz temporal sampling rate. Images (i) through (iv) correspond to singular vectors of index 1, 10, 3632, and 4036.
3.2 Image reconstruction of objects with different transducer count
Three (3) separate objects were created and used as phantoms to be reconstructed using the technique outlined in Section 2.4. A point, line, and multi-point source were created using experimental data acquired during the calibration scan and the linear superposition principle. However, in order to avoid reconstructing data directly from the imaging operator, realistic system noise (normally distributed) was added to the data set, g, to produce unique data. These time-dependent pressure measurements were then multiplied by the pseudoinverse obtained from component matrices (U T, S +, V) for each transducer count. Figure 3 shows the centre z-plane of the image reconstruction of each of the sources as a function of transducer count. The second column from the right shows the results of the iterative reconstruction technique applied to each of the 3 object types. A phantom of the objects is displayed in the right-most column.
Fig. 3 Reconstruction of a simulated point, line, and multi-point source for each of the 6 transducer arrangements. The iterative technique was implemented on the system with 30 transducers and 5 MHz temporal sampling rate (shown in the column second from right). Ideal image based on phantom is shown in the right-most column.
3.3 Image reconstruction of objects with different measurement space temporal sampling rate
The same objects used in Section 3.2 were again used as the experimental objects to reconstruct by the pseudoinverse of imaging operators (15 transducers) with variable measurement space temporal sampling rate. The result of the image reconstruction for each of the 8 temporal sampling rates (and phantom) is shown in Fig. 4 .
Fig. 4 Reconstruction of a simulated point, line, and multi-point source for each of the 8 measurement space temporal sampling rates (MHz).
3.4 Image reconstruction of a point source in real-time (1.4 seconds per frame)
The photoacoustic point-like source was imaged in 3D at 1.4 seconds intervals. Visual display of the point source is shown in Fig. 5 with each of the three orthogonal planes represented to illustrate the 3D nature of the image reconstruction process. The source had a starting position (10 mm, 11 mm, 4 mm) and was translated to position (6 mm, 11 mm, 4 mm).
Fig. 5 4D real-time photoacoustic imaging experiment where data acquisition and image reconstruction were performed in real time. Images were captured from a photoacoustic point-like source translated in the negative x-direction at a velocity of 0.40 mm/s. The interval between 3D photoacoustic images was 1.4 s and rate-limited by the data acquisition transfer speed, computational, and data storage overhead. The first row shows an xy-plane (z = 4 mm). The second row shows the same reconstruction in a zy-plane (x = 7 mm) and the third row displays an xz-plane (y = 11 mm). Each column represents the 3D image data acquired at a specific point in time (indicated along bottom).
4. Discussion
4.1 Estimate of matrix rank
The criterion in Eq. (4) worked well at producing a reliable estimate of the matrix rank. The ratio between consecutive singular values at the index of the matrix rank threshold was much greater than any other consecutive singular values. For example, the 30 transducer imaging operator had a maximum ratio of approximately 12 orders of magnitude while the ratio of any other 2 sequential singular values was never greater than 2 orders of magnitude. This was also seen in each of the 8 imaging operators with variable measurement space temporal sampling rate examined. Using Eq. (4), the effective matrix rank of the imaging operators obtained from a collection of transducer arrays of differing transducer count very closely followed a linear trend, as expected. In both Figs. 2(a) and 2(b), the slope estimated by the linear regression indicates the increase in matrix rank per change in transducer count and temporal sampling rate, respectively. While this is an important parameter, it does not translate directly to an improvement in image quality that would result from such an increase in singular vector count, as this would be dependent on the object being captured by the imaging system.
The results acquired by varying the measurement space temporal sampling rate were unexpected at the onset. The same linear increase in matrix rank was found as the temporal sampling rate increased. However, the linear trend was only observed in region (ii) of Fig. 2(b) for temporal sampling rates from 2.5 MHz to12.5 MHz. As estimated in a previous publication, the approximate resolution of the system is on the order of 1-2 mm [18]. According to Oraevsky et al. [19], the minimum temporal sampling rate required to capture information at a spatial scale of 1 mm is approximately 2.25 MHz. The two temporal sampling rates directly above (2.5 MHz) and below (1.25 MHz) this threshold identified the position at which the matrix rank deviated from the linear trend. Figure 2(c) was shows this abrupt deviation from linearity, which is not explicitly clear by simply interpreting Fig. 2(b).
When the measurement space dimension exceeded the object space dimension, as was the case for region (iii), the limit on the maximum number of resolvable singular vectors became 8000 rather than the measurement space dimension. Accordingly, this places a limit on the effective matrix rank. Again, this deviation from linearity is depicted clearly in Fig. 2(c).
Below the sampling rate required to maintain system resolution, it is perhaps intuitive that the matrix rank would decrease. However, in region (ii) the presence of a linear trend is less intuitive. Although the Nyquist sampling criterion is rigorously met in region (ii), there is an averaging of the information contained in the signal recorded in the imaging operator. While fundamentally the point source is still adequately resolved at lower temporal sampling rates (but still above the Nyquist criterion), dynamic range is decreased since fewer measurement points are used to describe object space at each grid point. In this context, a greater dynamic range signifies an increase in the number of bits used to describe the intensity of a voxel. When the temporal sampling rate is increased, the effective number of bits of used to represent each voxel is increased. Similar behavior is observed when applying a backprojection algorithm to the raw measurement data. As the degree of oversampling increases, a greater number of measurement points contribute to each voxel resulting in greater dynamic range.
The singular vectors plotted in Fig. 2(d) are shown in order to aid in the visualization of the vector content in a representative imaging operator. Coherent geometry is clear in images (i) and (ii), corresponding to singular vector indices 1 and 10. While a pattern is not discernable in images (iii) and (iv), it does not indicate that the vectors are void of information. First, only the centre plane is shown making it difficult to visualize geometry that could be coherent out of the plane. But more importantly, the practical limit of the voxel step-size makes it difficult to capture and observe the high frequency content that may actually reside in the higher order vectors. While image (iv) is beyond the matrix rank and consequently shows a singular vector in the null space, it is not clear how this impacts the limits of image reconstruction for this PA system.
The matrix rank is essential to applying generalized matrix inverse (GMI) reconstruction algorithms as it precisely indicates the object space singular vectors that contribute to resolving the unknown object. This has been investigated in other imaging modalities such as Single Photon Emission Computed Tomography (SPECT), indicating predictably that image quality is not improved when incorporating the use of singular vectors beyond the threshold outlined in Eq. (4) [17]. While the remaining object space singular vectors that exist in the null space of the imaging operator are linearly independent of the vectors above the threshold, the system cannot produce any object that is a combination of these vectors unless the reconstruction algorithm employs a priori information about the object (as can be achieved with iterative reconstruction methods [20]). However, even singular vectors within the limits of the matrix rank can be attributed to noise when the imaging operator is acquired experimentally. In these experimental imaging operators, if all vectors were included in the reconstruction up to the matrix rank, the noise in the experimental data set became correlated to the vectors within the threshold of the matrix rank.
4.2 Image reconstruction of objects with different transducer count
For the reconstruction of a point, line, and multi-point objects shown in Fig. 3, the image fidelity increased as the transducer count increased. The reconstruction of the point source with 5 transducers was not a reliable reproduction of the phantom object. However, the task was still adequately completed if using the FWHM criterion. That is, no other voxel in the image exceeded 50% intensity of the voxel containing the point source. Yet, the same was not true for objects of increasing complexity. Both the line and multi-point source were not adequately reconstructed by the lower transducer count systems due to view aliasing. It should be noted that, for example, the strong artifacts seen in the 5 transducer count reconstructions were actually still present in the 30 transducer count reconstructions but were much less significant due to the acquisition of a greater number of projections. As illustrated by reconstruction of the different object types, the number of singular vectors required to successfully construct an image was task specific. As object complexity increased, the matrix rank required to produce an accurate image of the object had to be increased as well.
Reconstruction utilizing the iterative reconstruction technique resulted in images of relatively poor quality. This was performed only for the imaging system with 30 transducers and 5 MHz temporal sampling rate. The computation time required to complete the task was approximately 74 seconds. It was thought this technique may produce higher quality images in comparison to the pseudoinverse technique, however, this technique utilizes only the raw imaging operator and it is expected that the effects of noise are significant in comparison to a reconstruction technique that regularizes the pseudoinverse. As well, the reconstruction time is prohibitive with respect to a system oriented towards real-time imaging.
4.3 Image reconstruction of objects with different measurement space temporal sampling rate
Below the Nyquist threshold, the point source did not reconstruct well. Although the location of the point source was correctly reconstructed, the images were corrupted with streaking and ray aliasing artifacts. Likewise, objects of greater complexity (i.e. line and grid of point sources) did not reconstruct with good fidelity. However, the reconstruction quality improved for all three object types as the temporal sampling rate increased. The best reconstruction performance was achieved for the highest temporal sampling rate of 25 MHz for an array of 15 transducers. For all three object types, a distinct transition in reconstruction quality was made at the threshold between regions (i) and (ii) of Fig. 2(b). Below the Nyquist threshold, the images did not accurately depict the object. This was observed most clearly in the case of the point source that exhibited streaking artifacts in each of the three images belonging to region (i). Similar results were observed in the reconstruction of the array of point sources. However, the artifacts were less obvious as the streaking originated from many points, obscuring the overall effect from any particular point source in the grid. Certainly, in all three object types, reliable estimates of the object were not made until the temporal sampling rate exceeded the Nyquist threshold. This corresponded to a matrix rank of approximately 2000 or more.
4.4 Real-time imaging of a photoacoustic point-like source
The reconstruction of the point-like source was performed successfully at a period of 1.4 seconds per frame. While the object was not complex, Fig. 5 demonstrates that photoacoustic signal acquisition, transfer, image reconstruction and display can be performed successfully on real objects. It should be noted that the images shown in Fig. 5 have not been modified through any means such as thresholding. Each 3D image was reconstructed from data captured after a single laser pulse. Therefore, each 3D image effectively represented a snap shot of the point source during a 9 ns time interval (the pulse length of the laser). This unprecedented capability should provide a path toward new applications such as 3D photoacoustic lifetime imaging and 3D tracking of moving targets (e.g. extensions of the 2D techniques described by Ashkenazi [21] and Su [22], respectively).
4.5 Implications to previous work
For larger matrices, the image reconstruction with the pseudoinverse significantly reduces the image reconstruction time when off-line image reconstruction is tolerable. In previous 4D photoacoustic imaging experiments, where iterative image reconstruction was utilized the image reconstruction time per frame was several minutes [12]. With knowledge of the imaging operator and a one-time computation of the pseudoinverse, the image reconstruction is now reduced to, at most, a few seconds per frame.
4.6 Technical considerations for faster 3D frame rates
Real-time imaging and display encompasses both the data acquisition time as well as the reconstruction of an image from the data acquired. Until now, the time limiting process was image reconstruction because iterative reconstruction techniques were utilized that spanned minutes to produce an image. Solving directly for the image vector has improved the reconstruction time to approximately 0.05 s and the data transfer speed has become the rate limiting process. For our current system, the data transfer rate permits 5120 points from 32 channels to be transferred in approximately 1.2 s. In order to reduce this time, modifications to the acquisition process can be made. First, since only 2000 of the total 5120 points per channel are used in the imaging operator, implementing a delay between the laser pulse and the start of data capture would provide a reduction in data transfer time from 1.2 s to 0.47 s. Second, improvements could be made to transfer the data from each card in parallel. With the current system, this would provide an additional 4-fold decrease in transfer time. Separately, the transfer process could be entirely changed to utilize a proprietary parallel data bus at 50 MHz that could reduce the transfer times to ~1.5 ms. With the improvements suggested, both acquisition and image reconstruction could be completed in under 0.1 s.
The utilization of a solution to Eq. (1) to perform image reconstruction is a significant transition towards 3D photoacoustic imaging where the data is captured and images are reconstructed in real time (e.g. Fig. 5). Because the pseudoinverse of the imaging operator can be computed before data acquisition, only matrix multiplication of the pseudoinverse to the measured data set is required to produce the image. Even with no software optimization, for a matrix size of 1792 x 6000, the matrix multiplication typically required approximately 0.05 seconds in MATLAB. This computation time could be reduced with GPU hardware acceleration. Therefore, real-time 3D photoacoustic imaging at the laser repetition rate (10 Hz) should be achievable using pseudoinverse image reconstruction with a known imaging operator and improved computational capabilities.
5. Conclusion
Our PAI system was modified to contain 30 transducers arranged in a hemispherical geometry. An imaging operator was experimentally acquired with a step-size of 1.5 mm within an object space of 30x30x30 mm3. By decomposing the imaging operator and analyzing the magnitude of the singular values, the effective matrix rank of the imaging operator was estimated. The results clearly show a linear increase in the effective matrix rank as both the transducer count and measurement space temporal sampling rate were increased. Image reconstruction of a variety of objects was performed by computing the pseudoinverse of each decomposed imaging operator and multiplying it with the measured data set of each phantom object. As expected, image fidelity was greater in cases with more transducers and increased temporal sampling rate. Image reconstruction was also performed in a period of 1.4 s per frame on a moving point-like source to demonstrate real-time acquisition, reconstruction and display of 3D photoacoustic images. Extension of the real-time 3D photoacoustic imaging to a frame rate of 10Hz should be easily achieved with straightforward improvements to data transfer time and computation speed.
Acknowledgments
The authors would like to acknowledge Lynn Keenliside and John Patrick for their technical support. 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 Ministry of Research and Innovation (ORF-RE) through the Ontario Preclinical Imaging Consortium (OPIC), the Canadian Foundation for Innovation (CFI), Natural Sciences and Engineering Research Council (NSERC), the Canadian Institutes for Health Research (CIHR), the Lawson Health Research Institute, and MultiMagnetics Inc. Mark Anastasio was supported in part by National Institutes of Health (NIH) awards EB010049 and EB009715.
References and links
1. L. V. Wang, “Ultrasound-mediated biophotonic imaging: a review of acousto-optical tomography and photo-acoustic tomography,” Dis. Markers 19(2-3), 123–138 (2003-2004). [PubMed]
2. 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]
3. R. A. Kruger, D. R. Reinecke, and G. A. Kruger, “Thermoacoustic computed tomography—technical considerations,” Med. Phys. 26(9), 1832–1837 (1999). [CrossRef] [PubMed]
4. 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,” Proc. SPIE 4434, 74–80 (2001). [CrossRef]
5. M. Frenz, K. P. Kostli, G. Paltauf, H. Schmidt-Kloiber, and H. P. Weber, “Reconstruction technique for optoacoustic imaging,” in Biomedical Optoacoustics II, January 23–24, 2001 (SPIE), 130–137.
6. C. G. A. Hoelen and F. F. M. de Mul, “Image reconstruction for photoacoustic scanning of tissue structures,” Appl. Opt. 39(31), 5872–5883 (2000). [CrossRef] [PubMed]
7. 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]
8. 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]
9. M. H. Xu and L. V. Wang, “Universal back-projection algorithm for photoacoustic computed tomography,” Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 71(1), 016706 (2005). [CrossRef] [PubMed]
10. M. Roumeliotis, R. Z. Stodilka, M. A. Anastasio, G. Chaudhary, H. Al-Aabed, E. Ng, A. Immucci, and J. J. L. Carson, “Analysis of a photoacoustic imaging system by the crosstalk matrix and singular value decomposition,” Opt. Express 18(11), 11406–11417 (2010). [CrossRef] [PubMed]
11. 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]
12. 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]
13. 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]
14. 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]
15. K. Konstantinides and K. Yao, “Statistical analysis of effective singular values in matrix rank determination,” IEEE Trans. Acoust. Speech Signal Process. 36(5), 757–763 (1988). [CrossRef]
16. K. Konstantinides, “Threshold bounds in SVD and a new iterative algorithm for order selection in Ar models,” IEEE Trans. Signal Process. 39(5), 1218–1221 (1991). [CrossRef]
17. D. J. Kadrmas, E. C. Frey, and B. M. W. Tsui, “An SVD investigation of modeling scatter in multiple energy windows for improved SPECT images,” IEEE Trans. Nucl. Sci. 43(4), 2275–2284 (1996). [CrossRef] [PubMed]
18. 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–23, 2008 (SPIE).
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 3601, 256–267 (1999). [CrossRef]
20. J. Provost and F. Lesage, “The application of compressed sensing for photo-acoustic tomography,” IEEE Trans. Med. Imaging 28(4), 585–594 (2009). [CrossRef] [PubMed]
21. S. Ashkenazi, “Photoacoustic lifetime imaging of dissolved oxygen using methylene blue,” J. Biomed. Opt. 15(4), 040501 (2010). [CrossRef] [PubMed]
22. J. Su, A. Karpiouk, B. Wang, and S. Emelianov, “Photoacoustic imaging of clinical metal needles in tissue,” J. Biomed. Opt. 15(2), 021309 (2010). [CrossRef] [PubMed]
