## Analysis of a photoacoustic imaging system by the crosstalk matrix and singular value decomposition

Optics Express, Vol. 18, Issue 11, pp. 11406-11417 (2010)

http://dx.doi.org/10.1364/OE.18.011406

Acrobat PDF (1065 KB)

### Abstract

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 OSA

## 1. Introduction

### 1.1 Background

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]

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]

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]

### 1.2 Singular value decomposition

**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):

**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

*et al*. initially as a way to recover and differentiate among Fourier coefficients used to describe an object in the frequency domain [16

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]

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]

*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

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]

*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 mm

^{3}). 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

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]

## Methods

### 2.1 Photoacoustic imaging system

*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 mm

^{3}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).

### 2.2 System calibration scan

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]

^{3}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×mm

^{3}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

**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

## 3. Results

### 3.1 Crosstalk sensitivity and aliasing

^{3}scan (10×10×10 voxels).

^{2}and is 16× 16 mm

^{2}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.

### 3.2 Singular value decomposition: Singular vectors

**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 mm

^{3}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.

## 4. Discussion

### 4.1 Crosstalk sensitivity and aliasing

**17**(17), 15228–15238 (2009). [CrossRef] [PubMed]

### 4.2 Singular value decomposition

**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.

### 4.3 Computation considerations

### 4.4 Imaging considerations

**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.

## 5. Conclusion

^{3}. The second was completed at a step size of 3 mm within an object space of 30×30×30 mm

^{3}. 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

## References and links

1. | M. Xu and L. V. Wang, “Photoacoustic imaging in biomedicine,” Rev. Sci. Instrum. |

2. | T. Lu, J. Jiang, Y. Su, R. K. Wang, F. Zhang, and J. Yao, Photoacoustic imaging: Its current status and future development” in |

3. | G. J. Diebold, T. Sun, and M. I. Khan, “Photoacoustic monopole radiation in one, two, and three dimensions,” Phys. Rev. Lett. |

4. | G. Paltauf, J. A. Viator, S. A. Prahl, and S. L. Jacques, “Iterative reconstruction algorithm for optoacoustic imaging,” J. Acoust. Soc. Am. |

5. | P. Liu, “The P-transform and photoacoustic image reconstruction,” Phys. Med. Biol. |

6. | C. G. A. Hoelen and F. F. de Mul, “Image reconstruction for photoacoustic scanning of tissue structures,” Appl. Opt. |

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 |

8. | M. Xu and L. V. Wang, “RF-induced thermoacoustic tomography” in |

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. |

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. |

11. | D. Modgil, M. A. Anastasio, and P. J. La Riviere, Photoacoustic image reconstruction in an attenuating medium using singular value decomposition” in |

12. | D. W. Wilson and H. H. Barrett, “Decomposition of images and objects into measurement and null components,” Opt. Express |

13. | P. Ephrat and J. J. L. Carson, Measurement of photoacoustic detector sensitivity distribution by robotic source placement” in |

14. | P. Ephrat, M. Roumeliotis, F. S. Prato, and J. J. Carson, “Four-dimensional photoacoustic imaging of moving targets,” Opt. Express |

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 |

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 |

17. | H. H. Barrett and H. Gifford, Cone-beam tomography with discrete data sets” in |

18. | J. Qi and R. H. Huesman, “Wavelet crosstalk matrix and its application to assessment of shift-variant imaging systems,” IEEE Trans. Nucl. Sci. |

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 |

**OCIS Codes**

(170.0110) Medical optics and biotechnology : Imaging systems

(170.5120) Medical optics and biotechnology : Photoacoustic imaging

**ToC Category:**

Medical Optics and Biotechnology

**History**

Original Manuscript: March 15, 2010

Revised Manuscript: May 7, 2010

Manuscript Accepted: May 7, 2010

Published: May 14, 2010

**Virtual Issues**

Vol. 5, Iss. 10 *Virtual Journal for Biomedical Optics*

**Citation**

Michael Roumeliotis, Robert Z. Stodilka, Mark A. Anastasio, Govind Chaudhary, Hazem Al-Aabed, Eldon Ng, Andrea Immucci, and Jeffrey J.L. Carson, "Analysis of a photoacoustic imaging system by the crosstalk matrix and singular value decomposition," Opt. Express **18**, 11406-11417 (2010)

http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-18-11-11406

Sort: Year | Journal | Reset

### References

- M. Xu and L. V. Wang, “Photoacoustic imaging in biomedicine,” Rev. Sci. Instrum. 77(4), 041101 (2006). [CrossRef]
- 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.
- 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]
- 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]
- P. Liu, “The P-transform and photoacoustic image reconstruction,” Phys. Med. Biol. 43(3), 667–674 (1998). [CrossRef] [PubMed]
- 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]
- 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, June18,2001- June 21 (SPIE), 74–80.
- 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), October23,2002- October 26 (Institute of Electrical and Electronics Engineers Inc), 1211–1212.
- 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]
- 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]
- 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 Sensing2009 (SPIE - The International Society for Optical Engineering), 71771B (7 pp.).
- 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]
- 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, January20,2008- January 23 (SPIE), Society of Photo-Optical Instrumentation Engineers (SPIE).
- 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]
- 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]
- 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]
- 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.
- 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]
- 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).

## Cited By |
Alert me when this paper is cited |

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

« Previous Article | Next Article »

OSA is a member of CrossRef.