Spatiotemporal antialiasing in photoacoustic computed tomography

Abstract

Among the various aspects of the present disclosure is the provision of systems and methods of imaging using photoacoustic computed tomography.

Description

FIELD

Certain aspects generally pertain to photoacoustic imaging and, more specifically, to photoacoustic computed tomography systems and methods.

BACKGROUND

Photoacoustic tomography provides high-resolution images beyond the optical diffusion limit by combining optical absorption contrast and ultrasonic spatial resolution. By converting highly scattered photons into ultrasonic waves, which are much less scattered than light in biological tissues, photoacoustic tomography can be a useful technique for forming high-resolution images of the optical properties of biological tissues.

SUMMARY

Certain aspects pertain to photoacoustic computed tomography with spatiotemporal antialiasing techniques as can be used, for example, to image biological tissues.

Certain aspects pertain to photoacoustic computed tomography methods that include acquiring photoacoustic data recorded by one or more data acquisition devices of a photoacoustic computed tomography system. The methods also include applying location-based temporal filtering to the photoacoustic data acquired. The location-based temporal filtering is based at least in part on geometry of the ultrasonic transducer array of the photoacoustic computed tomography system. The methods also include reconstructing one or more photoacoustic images from the filtered photoacoustic data. In some cases, an upper cutoff frequency is determined based at least in part on the geometry of the ultrasonic transducer array and one or more lowpass filters with the upper cutoff frequency determined is applied. According to one aspect, spatial interpolation may also be applied after location-based temporal filtering.

Certain aspects pertain to photoacoustic computed tomography methods that include determining whether a reconstruction location is at or inside, or outside an alias-free region. If it is determined that the reconstruction location is outside the alias-free region, location-based temporal filtering is applied to the photoacoustic data acquired. The application of location-based temporal filtering includes applying one or more lowpass filters with an upper cutoff frequency based at least in part on the geometry of the ultrasonic transducer array. The methods also include reconstructing one or more photoacoustic images from the filtered photoacoustic data. According to one aspect, spatial interpolation may also be applied after location-based temporal filtering.

Certain aspects pertain to non-transitory computer readable media for generating one or more photoacoustic images of a specimen being imaged, the non-transitory computer readable media, when read by one or more processors, is configured to perform operations. The operations performed include: acquiring photoacoustic data recorded by one or more data acquisition devices of a photoacoustic computed tomography system, applying location-based temporal filtering to the photoacoustic data acquired, wherein the location-based temporal filtering is based at least in part on geometry of the ultrasonic transducer array of the photoacoustic computed tomography system, and reconstructing one or more photoacoustic images from the filtered photoacoustic data.

Certain aspects pertain to non-transitory computer readable media for generating one or more photoacoustic images of a specimen being imaged, the non-transitory computer readable media, when read by one or more processors, is configured to perform operations. The operations performed include: (i) determining whether a reconstruction location is at or inside, or outside an alias-free region, (ii) if it is determined that the reconstruction location is outside the alias-free region, applying location-based temporal filtering to the photoacoustic data acquired, wherein the applying location-based temporal filtering includes applying one or more lowpass filters with an upper cutoff frequency based at least in part on the geometry of the ultrasonic transducer array, and (iii) reconstructing one or more photoacoustic images from the filtered photoacoustic data.

These and other features are described in more detail below with reference to the associated drawings.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1A is a schematic diagram depicting an analysis of spatial aliasing in spatial sampling by the full-ring ultrasonic transducer array, according to implementations.

FIG. 1B is a schematic diagram depicting an analysis of spatial aliasing in image reconstruction, according to one implementation.

FIG. 1C is a schematic diagram illustrating three regions in the field-of-view of the full-ring transducer array that represent regions of different types of spatial aliasing, according to one implementation.

FIG. 2A is a schematic diagram depicting spatial aliasing appearing during image reconstruction in a PACT system with spatiotemporal antialiasing with a full-ring ultrasonic transducer array, according to one implementation.

FIG. 2B is a schematic diagram representing mitigating spatial aliasing in image reconstruction using spatial interpolation, according to an implementation.

FIG. 2C is a schematic diagram representing mitigating spatial aliasing in spatial sampling using location-dependent temporal filtering, according to an implementation.

FIG. 3 is a schematic illustration of a portion of a surface of a linear ultrasonic transducer array, according to an implementation.

FIG. 4 is a schematic illustration of the portion of the linear ultrasonic transducer array in FIG. 3 depicting the imaging domain S₀, the one-way Nyquist zone S₁, and the two-way Nyquist zone S₂, according to an aspect.

FIG. 5 is a schematic diagram of components of a PACT system with spatiotemporal antialiasing, according to certain implementations.

FIG. 6 is signal flow diagram of signals communicated between components of a PACT system with spatiotemporal antialiasing, according to implementations.

FIG. 7 is a flowchart depicting operations of a PACT method with spatiotemporal antialiasing that implements one or more antialiasing techniques, according to an implementation.

FIG. 8 is a flowchart depicting operations of a PACT method with spatiotemporal antialiasing that applies location-based temporal filtering, according to an implementation.

FIG. 9A is a flowchart of operations of a universal back-projection process that can be used to reconstruct either a 2D image or a 3D image, according to an implementation.

FIG. 9B is a flowchart of additional operations of the universal back-projection process in FIG. 9A as used for the 3D image, according to an implementation.

FIG. 10A is an illustration of a numerical simulation of a transducer array and a simple phantom, according to an implementation.

FIG. 10B is an illustration of a numerical simulation of the transducer array in FIG. 10A with a complex numerical phantom, according to an implementation.

FIG. 11 depicts ground truth slices, according to an implementation.

FIG. 12 depicts reconstructed images, according to an implementation.

FIG. 13 is a graph of standard deviations of pixel values in regions of interest in reconstructed images in FIG. 11, according to an implementation.

FIG. 14 is an illustration of ground truth of a simple initial pressure distribution for a simple phantom, according to an implementation.

FIG. 15A is a reconstructed image of the simple phantom object in FIG. 14 using a first method including UBP reconstruction, according to an implementation.

FIG. 15B is a reconstructed image of the simple phantom object in FIG. 14 using a second method including UBP reconstruction and spatial interpolation, according to an implementation.

FIG. 15C is a reconstructed image of the simple phantom object in FIG. 14 using a third method including UBP reconstruction, spatial interpolation, and location-dependent temporal filtering, according to an implementation.

FIG. 16 is a graph with a plot of standard deviations of the regions of interest in FIG. 15A using (1) UBP reconstruction (2) UBP reconstruction and spatial interpolation, and (3) UBP reconstruction, spatial interpolation, and location-dependent temporal filtering (TF).

FIG. 17 is a graph with a plot of profiles of lines using the three methods: (1) UBP reconstruction (2) UBP reconstruction and spatial interpolation, and (3) UBP reconstruction, spatial interpolation (SI), and location-dependent temporal filtering (TF).

FIG. 18 is a graph with a plot of profiles of lines Q using image reconstruction (1) with UBP reconstruction (2) UBP reconstruction and spatial interpolation, and (3) UBP reconstruction, spatial interpolation, and location-dependent temporal filtering.

FIG. 19 is an illustration of ground truth of a simple initial pressure p₀ distribution for a first complex phantom, according to an implementation.

FIG. 20A is a reconstructed image of the first complex phantom object shown in FIG. 19 using a first antialiasing method including universal back projection reconstruction, according to an implementation.

FIG. 20B is a reconstructed image of the first complex phantom object in FIG. 19 using the second method including UBP reconstruction and spatial interpolation, according to an implementation.

FIG. 20C is a reconstructed image of the first complex phantom object in FIG. 19 using a third method including UBP reconstruction, spatial interpolation, and location-dependent temporal filtering, according to an implementation.

FIG. 21 is a graph with a plot of standard deviations of three regions of interest obtained using: (1) a first method with UBP reconstruction alone, (2) a second method with UBP reconstruction and spatial interpolation, and (3) a third method with UBP reconstruction, spatial interpolation, and location-dependent temporal filtering, according to implementations.

FIG. 22 is a graph with a plot of profiles of line P obtained using three methods, according to implementations.

FIG. 23 is a graph with a plot of profiles of line Q obtained using three methods, according to implementations.

FIG. 24 is an illustration of ground truth of a simple initial pressure p₀ distribution for a second complex phantom, according to an implementation.

FIG. 25A is a reconstructed image of the second complex phantom object shown in FIG. 24 using a first antialiasing method including universal back projection reconstruction, according to an implementation.

FIG. 25B is a reconstructed image of the second complex phantom object in FIG. 24 using the second method including UBP reconstruction and spatial interpolation, according to an implementation.

FIG. 25C is a reconstructed image of the second complex phantom object in FIG. 24 using a third method including UBP reconstruction, spatial interpolation, and location-dependent temporal filtering, according to an implementation.

FIG. 26 is a graph with a plot of standard deviations of three regions of interest: (1) a first method with UBP reconstruction alone, (2) a second method with UBP reconstruction and spatial interpolation, and (3) a third method with UBP reconstruction, spatial interpolation (SI), and location-dependent temporal filtering (TF), according to implementations.

FIG. 27 is a graph with a plot of profiles of line P using the three methods, according to implementations.

FIG. 28 is a graph with a plot of profiles of line Q using the three methods, according to implementations.

FIG. 29A is a reconstructed image using universal back propagation without either spatial interpolation or location-based temporal filtering, according to an implementation.

FIG. 29B is a reconstructed image of the same region in FIG. 29A using universal back propagation with spatial interpolation, according to an implementation.

FIG. 29C is a reconstructed image of the same region in FIG. 29A using universal back propagation with spatial interpolation and location-based temporal filtering, according to an implementation.

FIG. 30 is a graph with a plot of profiles of line P using the three methods, according to implementations.

FIG. 31 is a graph with a plot of profiles of line Q using the three methods, according to implementations.

FIG. 32 is a graph with a plot of the difference between the normalized spectra

$\frac{\left|F\left(f_1{f}_2\right)\left(f\right)\right|}{\underset{f}{\max }\left|F\left(f_1{f}_2\right)\left(f\right)\right|}\mathrm{and}\frac{\left|F\left(f_2\right)\left(f\right)\right|}{\underset{f}{\max }\left|F\left(f_2\right)\left(f\right)\right|},$ according to one implementation.

FIG. 33A is a graph with a plot of a normalized estimated point source response, according to one implementation.

FIG. 33B is a graph with a plot of a normalized frequency spectrum of the estimated point source response, according to an implementation.

FIG. 33C is a graph with a plot of a normalized frequency spectrum of the noise, according to an implementation.

FIG. 33D is a graph with a plot of a frequency-dependent signal-to-noise ratio, according to an implementation.

FIG. 33E is a graph with a plot of frequency-dependent SNR of the derivative of the point source response, according to an implementation.

FIG. 33F is a graph with a plot of normalized frequency spectra, according to an implementation.

FIG. 34 is a schematic diagram of a full-ring transducer array, according to an implementation.

FIG. 35 is a graph with a plot of examples of calculated errors in sample size approximation, according to an implementation.

FIG. 36 is a graph with a plot of examples of calculated errors in sample size approximation, according to an implementation.

These and other features are described in more detail below with reference to the associated drawings.

DETAILED DESCRIPTION

Different aspects are described below with reference to the accompanying drawings. The features illustrated in the drawings may not be to scale. In the following description, numerous specific details are set forth in order to provide a thorough understanding of the presented embodiments. The disclosed embodiments may be practiced without one or more of these specific details. In other instances, well-known operations have not been described in detail to avoid unnecessarily obscuring the disclosed embodiments. While the disclosed embodiments will be described in conjunction with the specific embodiments, it will be understood that it is not intended to limit the disclosed embodiments.

I. Introduction to Photoacoustic Computed Tomography (PACT)

Photoacoustic computed tomography (PACT) is an imaging modality that ultrasonically images optical contrast via the photoacoustic effect. PACT can provide tomographic images (e.g., two-dimensional (2D) section(s) and/or 3D volume(s) constructed from 2D sections) of biological tissues and other samples. By converting highly scattered photons into ultrasonic waves, which are much less scattered than light in biological tissues, PACT can form high-resolution images of the optical properties of these tissues at different depths. Some examples of photoacoustic imaging of biological tissues can be found in Kruger, R. A., Kuzmiak, C. M., Lam, R. B., Reinecke, D. R., Del Rio, S. P. and Steed, D., “Dedicated 3D photoacoustic breast imaging,” Med. Phys., vol. 40, no. 11 (2013) (hereinafter “Kruger”), Tzoumas, S., et al., “Eigenspectra optoacoustic tomography achieves quantitative blood oxygenation imaging deep in tissues,” Nat. Commun., vol. 7 (2016), Dean-Ben, X. L. et al., “Functional optoacoustic neuro-tomography for scalable whole-brain monitoring of calcium indicators,” Light Sci. Appl., vol. 5, no. 12, p. e16201, (2016), Li, L. et al., “Single-impulse panoramic photoacoustic computed tomography of small-animal whole-body dynamics at high spatiotemporal resolution,” Nat. Biomed. Eng., vol. 1, no. 5, p. 0071 (2017), Matsumoto, Y., et al., “Label-free photoacoustic imaging of human palmar vessels: a structural morphological analysis,” Sci. Rep., vol. 8, no. 1, p. 786 (2018), Lin, L. et al., “Single-breath-hold photoacoustic computed tomography of the breast,” Nat. Commun., vol. 9, no. 1, p. 2352 (2018), Li, L., “Small near-infrared photochromic protein for photoacoustic multi-contrast imaging and detection of protein interactions in vivo,” Nat. Commun., vol. 9, no. 1, p. 2734 (2018), and Wu, Z. et al., “A microrobotic system guided by photoacoustic computed tomography for targeted navigation in intestines in vivo,” Sci. Robot., vol. 4, no. 32, p. 0613, (2019), which are hereby incorporated by reference in their entireties. Some examples of PACT methods and systems can be found in U.S. patent application Ser. No. 16/798,204, titled “PHOTOACOUSTIC COMPUTED TOMOGRAPHY (PACT) SYSTEMS AND METHODS,” filed on Feb. 21, 2020; PCT publication WO2018102446, titled “SINGLE-IMPULSE PANORAMIC PHOTOACOUSTIC COMPUTED TOMOGRAPHY (SIP-PACT), filed on Nov. 29, 2017; U.S. Patent Publication US 2016/0262628, titled “PHOTOACUSTIC COMPUTED TOMOGRAPHY WITH AN ACOUSTIC REFLECTOR.” filed on Nov. 12, 2014; which are hereby incorporated by reference in their entireties.

PACT methods generally include a data acquisition phase and an image reconstruction phase. During the data acquisition phase, photon-induced acoustic waves (sometimes referred to herein as “photoacoustic waves”) are detected by an ultrasonic transducer array or its scanning equivalent, such as a full-ring ultrasonic transducer array, a linear array, or two-dimensional array. During the image reconstruction phase, the detected acoustic signals are used to reconstruct the sample's optical absorption via inverse reconstruction algorithms. Some examples of inverse reconstruction algorithms that can be used include: (i) forward-model-based iterative methods, (ii) time-reversal methods, and (iii) the universal back-projection (UBP) method. Examples of the universal back-projection (UBP) method can be found in Kruger, Li, L. et al., “Single-impulse panoramic photoacoustic computed tomography of small-animal whole-body dynamics at high spatiotemporal resolution,” Nat. Biomed. Eng., vol. 1, no. 5, p. 0071 (2017), Lin, L. et al., “Single-breath-hold photoacoustic computed tomography of the breast,” Nat. Commun., vol. 9, no. 1, p. 2352 (2018), Xu M. and Wang, L. V., “Universal back-projection algorithm for photoacoustic computed tomography,” Phys. Rev. E, vol. 71, no. 1, p. 016706 (2005), Song, L., Maslov, K. I., Bitton, R., Shung, K. K., and Wang, L. V., “Fast 3-D dark-field reflection-mode photoacoustic microscopy in vivo with a 30-MHz ultrasound linear array,” J. Biomed. Opt., vol. 13, no. 5, p. 054028 (2008), Dean-Ben, X. L. and Razansky, D., “Portable spherical array probe for volumetric real-time optoacoustic imaging at centimeter-scale depths,” Opt. Express, vol. 21, no. 23, pp. 28062-28071 (2013), and Pramanik, M., “Improving tangential resolution with a modified delay-and-sum reconstruction algorithm in photoacoustic and thermoacoustic tomography,” JOSA A, vol. 31, no. 3, pp. 621-627 (2014), which are hereby incorporated by reference in their entireties. Examples of forward-model-based iterative methods can be found in Wang, K., et al. “Investigation of iterative image reconstruction in three-dimensional optoacoustic tomography,” Phys. Med. Biol., vol. 57, no. 17, p. 5399 (2012), Huang, Chao et al., “Full-wave iterative image reconstruction in photoacoustic tomography with acoustically inhomogeneous media,” IEEE Trans. Med. Imaging, vol. 32, no. 6, pp. 1097-1110 (June 2013), Mitsuhashi, K., et al., “A forward-adjoint operator pair based on the elastic wave equation for use in transcranial photoacoustic computed tomography,” SIAM J. Imaging Sci., vol. 10, no. 4, pp. 2022-2048 (2017), Treeby, B. E. and Cox, B. T., “k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields,” J. Biomed. Opt., vol. 15, no. 2, p. 021314, (2010), Mitsuhashi, K., Wang, K. and Anastasio, M. A., “Investigation of the far-field approximation for modeling a transducer's spatial impulse response in photoacoustic computed tomography,” Photoacoustics, vol. 2, no. 1, pp. 21-32 (2014), Han, Y., Ntziachristos, V. and Rosenthal, A., “Optoacoustic image reconstruction and system analysis for finite-aperture detectors under the wavelet-packet framework,” J. Biomed. Opt., vol. 21, no. 1, p. 016002 (2016), Arridge, S., et al., “Accelerated high-resolution photoacoustic tomography via compressed sensing,” ArXiv Prepr., ArXiv160500133 (2016), Han, Y. et al., “Three-dimensional optoacoustic reconstruction using fast sparse representation,” Opt. Lett., vol. 42, no. 5, pp. 979-982 (2017), Schoeder, S. et al. “Optoacoustic image reconstruction: the full inverse problem with variable bases,” Proc. R. Soc. A., vol. 474, no. 2219, p. 20180369 (2018), and Matthews, T. P., Poudel, J., Li, L., Wang, L. V. and Anastasio, M. A., “Parameterized Joint Reconstruction of the Initial Pressure and Sound Speed Distributions for Photoacoustic Computed Tomography,” SIAM J. Imaging Sci., vol. 11, no. 2, pp. 1560-1588 (2018), which are hereby incorporated by reference in their entireties. Examples of time reversal methods can be found in Xu, Y. and Wang, L. V., “Time reversal and its application to tomography with diffracting sources,” Phys. Rev. Lett., vol. 92, no. 3, p. 033902 (2004), Treeby, B. E., Zhang, E. Z. and Cox, B., “Photoacoustic tomography in absorbing acoustic media using time reversal,” Inverse Probl., vol. 26, no. 11, p. 115003 (2010), Cox, B. T. and Treeby, B. E., “Artifact trapping during time reversal photoacoustic imaging for acoustically heterogeneous media,” IEEE Trans. Med. Imaging, vol. 29, no. 2, pp. 387-396 (2010), Treeby, B. E., Jaros, J. and Cox, B. T., “Advanced photoacoustic image reconstruction using the k-Wave toolbox,” in Photons Plus Ultrasound: Imaging and Sensing 2016, vol. 9708, p. 97082P (2016), and Ogunlade, O. et al., “In vivo three-dimensional photoacoustic imaging of the renal vasculature in preclinical rodent models,” Am. J. Physiol.-Ren. Physiol., vol. 314, no. 6, pp. F1145-F1153 (2017), which are hereby incorporated by reference in their entireties. Examples of universal back-projection (UBP) method can be found in Kruger, Li, L. et al., “Single-impulse panoramic photoacoustic computed tomography of small-animal whole-body dynamics at high spatiotemporal resolution,” Nat. Biomed. Eng., vol. 1, no. 5, p. 0071 (2017), Lin, L. et al., “Single-breath-hold photoacoustic computed tomography of the breast,” Nat. Commun., vol. 9, no. 1, p. 2352 (2018), Xu M. and Wang, L. V., “Universal back-projection algorithm for photoacoustic computed tomography,” Phys. Rev. E, vol. 71, no. 1, p. 016706 (2005), Song, L., Maslov, K. I., Bitton, R., Shung, K. K., and Wang, L. V., “Fast 3-D dark-field reflection-mode photoacoustic microscopy in vivo with a 30-MHz ultrasound linear array,” J. Biomed. Opt., vol. 13, no. 5, p. 054028 (2008), Dean-Ben, X. L. and Razansky, D., “Portable spherical array probe for volumetric real-time optoacoustic imaging at centimeter-scale depths,” Opt. Express, vol. 21, no. 23, pp. 28062-28071 (2013), and Pramanik, M., “Improving tangential resolution with a modified delay-and-sum reconstruction algorithm in photoacoustic and thermoacoustic tomography,” JOSA A, vol. 31, no. 3, pp. 621-627 (2014), which are hereby incorporated by reference in their entireties.

PACT based on a full-ring ultrasonic transducer array is widely used for small animal wholebody and human organ imaging, thanks to its high in-plane resolution and full-view fidelity. Some discussion of PACT used to for small animal wholebody and/or human organ imaging can be found in Li, L. et al., “Single-impulse panoramic photoacoustic computed tomography of small-animal whole-body dynamics at high spatiotemporal resolution,” Nat. Biomed. Eng., vol. 1, no. 5, p. 0071 (2017), Xu, Y., Xu, M. and Wang, L. V., “Exact frequency-domain reconstruction for thermoacoustic tomography. II. Cylindrical geometry,” IEEE Trans. Med. Imaging, vol. 21, no. 7, pp. 829-833, 2002, which are hereby incorporated by reference in their entireties. To help avoid any artifacts occurring in image reconstruction in PACT, the ultrasonic transducer array should provide dense spatial sampling around the object to satisfy the Nyquist sampling theorem. That is, the spatial sampling interval on the specimen surface should be less than half of the lowest detectable acoustic wavelength. Otherwise, artifacts may appear in image reconstruction, a problem called spatial aliasing. In practice, due to the high cost of a transducer array with a large number of elements or limited scanning time, spatially sparse sampling is common.

PACT systems and methods with spatiotemporal antialiasing described herein may mitigate artifacts caused by spatial undersampling without having to physically increase the number of elements in the transducer array to satisfy Nyquist sampling criterion. Moreover, in some implementations, the PACT systems and/or methods with spatiotemporal antialiasing may compensate for undersampling without substantially compromising image resolution. PACT systems/methods with spatiotemporal antialiasing may be advantageous in reducing equipment costs by being able to utilize less expensive ultrasonic transducer arrays that sparsely sample while still avoiding artifacts and maintaining high resolution in photoacoustic images.

In the following sections, the sources of spatial aliasing in a full-ring geometry PACT and linear array PACT and applications to other array geometries are discussed. The PACT systems and methods with spatiotemporal antialiasing described here implement one or more techniques that may mitigate artifacts caused by spatial aliasing without physically increasing the number of elements in an ultrasonic transducer array.

Based on spatiotemporal analysis, two sources of spatial aliasing are: spatial sampling and image reconstruction. To classify different cases of spatial aliasing based on locations (of source points and reconstruction locations), three zones are defined in the imaging domain of the ultrasonic transducer array: (1) a detection zone S₀, (2) a one-way Nyquist zone S₁, and (3) a two-way Nyquist zone S₂. Zone S₁ is contained in S₀, while S₂ is contained in S₁. Spatial aliasing in spatial sampling does not appear for objects inside S₁, but appears for objects outside S₁. Spatial aliasing in image reconstruction does not appear for objects and reconstruction locations inside S₂, but appears for other combinations of objects and reconstruction locations. Zone S₂ is the original aliasing-free region, which can be extended to S₁ by one or more antialiasing techniques described herein. The aliasing free region can also be further extended from S₁ to S₁′ by one or more antialiasing techniques described herein.

Two categories of antialiasing techniques include: spatial interpolation and temporal filtering. Spatial interpolation can be used to eliminate spatial aliasing in image reconstruction without compromising the spatial resolution, effectively extending the alias-free zone from S₂ to S₁. Temporal filtering that uses an upper cutoff frequency based on Nyquist sampling criterion can remove aliasing in spatial sampling, while compromising the spatial resolution in the affected regions. Thus, this combination of spatial interpolation with temporal filtering using the upper cutoff frequency based on Nyquist sampling criterion removes spatial aliasing in both spatial sampling and image reconstruction, while compromising the spatial resolution.

In certain aspects, as a balance between spatial antialiasing and high resolution, an antialiasing filter is designed based on the reconstruction location's distance to the center of the ultrasonic transducer array. In a full-ring transducer array implementation, the reconstruction location's distance is to the full-ring array center. Temporal filtering that implements an antialiasing filter designed based on the reconstruction location's distance to the center of the ultrasonic transducer array is referred to herein as location-dependent temporal filtering (also sometimes referred to as “radius-dependent” temporal filtering).

In these aspects, spatial interpolation can extend the alias-free zone from S₂ to S₁ and applying location-dependent temporal filtering along, the alias-free zone is further extended from S₁ to S₁′. The lower the upper cutoff frequency used in the location-based temporal filtering, the larger the extended region S₁′ and the poor the imaging resolution. As provided in Sections IV and V, the results from numerical simulations and in vivo experiments show that location-dependent temporal filtering along with spatial interpolation may be effective in mitigating artifacts caused by spatial aliasing in either image reconstruction or spatial sampling.

Antialiasing techniques that include spatial interpolation and/or location-based temporal filtering are discussed below with reference to a transducer array having circular geometry (e.g., full-ring array) and a linear geometry, for which expressions of for zones S₂ and S₁ are provided. It would be understood that the disclosure is not so limiting and that these antialiasing techniques may be used for other array geometries.

II. PACT Systems and Methods with Spatiotemporal Antialiasing

Spatial aliasing may be generally classified into two categories: (1) aliasing in spatial sampling and (2) aliasing in image reconstruction. In certain implementations, PACT systems and methods with spatiotemporal antialiasing described herein implement antialiasing techniques including: (i) spatial interpolation and/or (ii) location-dependent temporal filtering that is dependent on the geometry of the ultrasonic transducer array. Spatial interpolation can mitigate (reduce or eliminate) spatial aliasing in image reconstruction. Location-dependent temporal filtering can mitigate spatial aliasing in spatial sampling. In one aspect, location-dependent temporal filtering implements an antialiasing filter with one or more lowpass filters having an upper cutoff frequency that is based on the distance between the reconstruction location and the ultrasonic transducer array. In one implementation, an antialiasing technique determines whether the reconstruction location (and associated source point) is within the present alias-free zone and if it is determined that the reconstruction location is outside the present alias-free zone, an antialiasing filter is applied to the raw data to extend the present alias-free zone to an expanded alias-free zone to include the reconstruction location. In certain implementations, the PACT methods and systems with spatiotemporal antialiasing mitigate aliasing artifacts in image reconstruction.

Although the spatiotemporal antialiasing techniques are described in this section with reference to a transducer array having a circular geometry and a transducer array having a linear geometry, it would be understood that the disclosure is not so limiting, and that these techniques would apply to other array geometries. For example, based on the basic 1D circular and linear geometries, one can apply 2D/3D geometries, such as, a planar geometry and the spherical geometry, through decomposition. Moreover, although the antialiasing techniques are sometimes described being implementing with the universal back-projection (UBP) described in Xu, M. and Wang, L. V., “Universal back-projection algorithm for photoacoustic computed tomography,” Phys. Rev. E, vol. 71, no. 1, p. 016706 (2005), it would be understood that other inverse reconstruction techniques can be used, such as those referenced in Section I. In addition, the antialiasing techniques described herein may apply to other image reconstruction methods. For example, spatial interpolation may be used in time reversal methods to generate a dense enough grid for numerical computation. Furthermore, location-dependent temporal filtering may be incorporated into a wave propagation model and be used in time reversal methods and/or iterative methods to mitigate aliasing in spatial sampling according to other aspect.

A. Photoacoustic Image Reconstruction

In a homogeneous medium, a photoacoustic wave can be expressed as:

$\begin{array}{cc}p\left(r,t\right)=\frac{1}{4{\pi c}^2}\frac{\partial }{\partial t}\left(\frac{1}{\mathrm{ct}}{\int }_V{\mathrm{dr}}^{\prime }{p}_0\left(r^{\prime }\right)\delta \left(t-\frac{||r^{\prime }-r||}{c}\right)\right).& \left(\mathrm{Eqn}.1\right)\end{array}$ A detailed discussion of the expression in Eqn. 1 can be found Wang, L. V. and Wu, H., Biomedical optics: principles and imaging. John Wiley & Sons (2012), Zhou, Y., Yao, J. and Wang, L. V., “Tutorial on photoacoustic tomography,” J. Biomed. Opt., vol. 21, no. 6, p. 061007, 2016, which are hereby incorporated by reference in their entireties. In Eqn. 1, p(r, t) is the pressure at location r and time t, c is the speed of sound, V is the volumetric space occupied by the tissue, and p₀(r′) is the initial pressure at r′. Eqn. 1 can be rewritten as:

$\begin{array}{cc}p\left(r,t\right)=\frac{1}{4{\pi c}^2}{\int }_V\frac{p_0\left(r^{\prime }\right)}{||r^{\prime }-r||}\frac{\partial }{\partial t}\delta \left(t-\frac{||r^{\prime }-r||}{c}\right){\mathrm{dr}}^{\prime }.& \left(\mathrm{Eqn}.2\right)\end{array}$ Discreizing Eqn 2 in space, provides:

$\begin{array}{cc}p\left(r_n,t\right)=\frac{1}{4{\pi c}^2}\sum _{m=1}^M{v}_m\frac{p_0\left(r^{\prime }\right)}{||r_m^{\prime }-r_n||}\frac{\partial }{\partial t}\delta \left(t-\frac{||r^{\prime }-r||}{c}\right),n=1,2,\dots ,N.& \left(\mathrm{Eqn}.3\right)\end{array}$ Eqn. 3 refers to m source points distributed at r_m′, to =1, 2, . . . M, and N point detection elements distributed at r_n, n=1, 2, . . . , N. The term v_m is the volume of the m-th source point. The response of an ultrasonic transducer array can be described as:{circumflex over (p)}(r_n,t)=p(r_n,t)*_th_e(t),n=1,2, . . . ,N.  (Eqn. 4)

Here, {circumflex over (p)}(r_n, t) is the pressure impinging on the n-th point detection element (transducer element) at time t, and h_e (t) is the ultrasonic transducer's electric impulse response. Substituting Eqn. 3 into Eqn. 4, the following is obtained:

$\begin{array}{cc}\hat{p}\left(r_n,t\right)=\frac{1}{4{\pi c}^2}\sum _{m=1}^M{v}_m\frac{p_0\left(r^{\prime }\right)}{||r_m^{\prime }-r_n||}{h}_e^{\prime }\left(t-\frac{||r^{\prime }-r||}{c}\right).& \left(\mathrm{Eqn}.5\right)\end{array}$ The term

$h_e^{\prime }\left(t-\frac{||r_m^{\prime }-r_n||}{c}\right)$ is a function of both time and space, where the first prime denotes a temporal derivative.

When photoacoustic signals received by the ultrasonic transducer array are digitized by one or more data acquisition systems (DAQs), an antialiasing filter (e.g. one or more low-pass filters) with an upper cutoff frequency selected for a sufficiently high temporal sampling rate may be used to mitigate temporal aliasing in certain implementations. For simplicity, the time variable is assumed to be continuous and spatial variables are discretized, which allows for a discussion of spatial sampling (sometimes referred to herein as “SS”). For different transducer detection geometries (e.g., planar, spherical, and cylindrical surfaces of transducer array), image flapping the initial pressure p₀ (r″) using the Universal back-projection (UBP) algorithm found in Xu, M. and Wang, L. V., “Universal back-projection algorithm for photoacoustic computed tomography,” Phys. Rev. E, vol. 71, no. 1, p. 016706 (2005), provides:

$\begin{array}{cc}{p}_0\left(r^{″}\right)={\int }_{{\Omega }_0}b\left(r,t=\frac{||r^{″}-r||}{c}\right)\frac{d\Omega }{{\Omega }_0},& \left(\mathrm{Eqn}.6\right)\end{array}$ where the back-projection term,

$b\left(r,t\right)=2p\left(r,t\right)-2t\frac{\partial p\left(r,t\right)}{\partial t},d\Omega =\frac{\mathrm{dS}}{||r^{″}-r{||}^2}\frac{n_S\left(r\right)·\left(r^{″}-r\right)}{||r^{″}-r||},$ is the solid angle for the detection element at r with respect to reconstruction location r″, dS is the detection element surface area, and n_S(r) is the ingoing normal vector. The total solid angle is denoted as Ω₀. The actual or true pressure p(r, t) may be approximated by the detected pressure {circumflex over (p)}(r_n, t), leading to the following discretized form of Eqn. 6:

$\begin{array}{cc}{\hat{p}}_0\left(r^{″}\right)\approx \sum _{n=1}^N{w}_n\hat{b}\left(r_n,t=\frac{||r^{″}-r_n||}{c}\right).& \left(\mathrm{Eqn}.7\right)\end{array}$ Here, {circumflex over (p)}₀(r″) is the reconstructed initial pressure an

$\hat{b}\left(r_n,t\right)=2\hat{p}\left(r_n,t\right)-2t\frac{\partial \hat{p}\left(r_n,t\right)}{\partial t}$ is the back-projection term computed from the detected pressure by the ultrasonic transducer array. The weights w_n, n=1, 2, . . . , N come from

$\frac{d\Omega }{{\Omega }_0}$ in Eqn. 6.B. Spatiotemporal Analysis of Spatial Aliasing in Spatial Sampling (SS) and Image Reconstruction (IR)

To determine the source of spatial aliasing in a PACT system, spatiotemporal analysis is used to analyze the detection geometry of a transducer array. For a particular detection geometry, the spatial sampling frequency may be calculated using reciprocals of the distances between adjacent sampling positions of the transducer elements in the array. In this section, a spatiotemporal analysis of a full-ring transducer array with point detection elements is used as an example. It would be understood that this spatiotemporal analysis can be applied to other array geometries such as a linear array as discussed in Section II(C)(2). The full-ring ultrasonic transducer array has a radius of R and has N evenly-distributed detection elements. The upper cutoff frequency of the full-ring ultrasonic transducer array is ƒ_c, and the corresponding lower cutoff wavelength is

${\lambda }_c=\frac{c}{f_c}.$

An example of a method for estimating an upper cutoff frequency ƒ_c that would eliminate spatial aliasing due to spatial sampling in the entire detection zone (imaging domain) of the transducer array is provided in Section VII. In certain implementations, a PACT system/method implements an antialiasing filter with one or more low pass filters (e.g., a third-order lowpass Butterworth filter and a sinc filter, both with same cutoff frequency of ƒ_c) with an upper cutoff frequency that is lower than this estimated upper cutoff frequency to substantially mitigate aliasing due to spatial sampling while maintaining the best level of image resolution as discussed in Section II(C). The higher the cutoff frequency, the better the image resolution. In certain implementations, the highest possible cutoff frequency is selected, and thus have the best image resolution.

For a full-ring transducer array, the detection zone, S₀, of the full-ring transducer array is defined as:S₀={r|∥r∥≤R}.  (Eqn. 8)where the full-ring transducer array has a radius, R and r′ is a source point or a reconstruction location.

FIG. 1A, FIG. 1B, and FIG. 1C are schematic diagrams depicting an analysis of spatial aliasing in a PACT system with a full-ring ultrasonic transducer array, according to certain implementations. The diagrams include a detection zone, S₀, 104, with an outer circle representing the detection surface of a full-ring transducer array having a radius R and N evenly-distributed transducer elements. The center O of detection zone, S₀, 104, represents the origin of a coordinate system used in image reconstruction.

FIG. 1A is a schematic diagram depicting an analysis of spatial aliasing in spatial sampling (sometimes referred to herein as “SS”) by the full-ring ultrasonic transducer array, according to implementations. The diagram includes a representative transducer detection element with location r and a representative source point with location r′.

FIG. 1B is a schematic diagram depicting an analysis of spatial aliasing in image reconstruction (sometimes referred to herein as “IR”), according to one implementation. The diagram includes a detection (transducer) element with a detection element location r. The diagram also includes two reconstruction locations r″ and r₁″, and a source point location r′.

FIG. 1C is a schematic diagram illustrating three regions in the field-of-view of the full-ring transducer array that represent where different types of aliasing might arise, according to one implementation. The diagram includes an alias-free two-way Nyquist zone, S₂, 101, which contains all source points and reconstruction locations where image reconstruction, e.g., reconstruction including universal back projection (UBP), yields no spatial aliasing. The diagram also includes a one-way Nyquist zone, S₁, 102, within which aliasing does not exist in spatial sampling, but may exist in image reconstruction. The diagram also includes a third zone, S₀, which is the imaging domain. Within the third zone, S₀, 103 aliasing may exist due to spatial sampling.

First, aliasing from spatial sampling of the full-ring transducer array is analyzed. When the detection element at location r varies discretely, the step size along the perimeter of the full-ring transducer array is

$\frac{2\pi R}{N}.$ The tangential direction to the perimeter is shown as a dotted line in FIG. 1A, which is perpendicular to the vector r from the center O of the circle. When considering a source point at r′, and extending the line segment r−r′ as a dashed line extension in FIG. 1A. Vectors −r′ and r−r′ form an angle while vector r−r′ forms an angle α′ with the tangential dotted line that is perpendicular to the vector r. The angle formed by vectors −r and r′−r can be expressed as

${\alpha }^{\prime }-\frac{\pi }{2}.$ The local sampling step size of ∥r−r′∥ may be approximated as:

$\begin{array}{cc}\frac{2\pi R}{N}{\mathrm{cos\alpha }}^{\prime }.& \left(\mathrm{Eqn}.9\right)\end{array}$

The absolute value means the length of the local sampling step size, while the sign (positive or negative of the value in Eqn. 9 means the sampling direction. An example of a method that can be used to approximate the sampling step size is discussed in Section VIII. From Eqn. 5, at a given time t, and with the lower cutoff wavelength λ_c, the Nyquist criterion can be expressed as:

$\begin{array}{cc}\frac{2\pi R\left|{\mathrm{cos\alpha }}^{\prime }\right|}{N}<\frac{{\lambda }_c}{2}.& \left(\mathrm{Eqn}.10\right)\end{array}$ The Law of Sines can be used to transform this inequality to a constraint for the source point location r′, to:

$\begin{array}{cc}\frac{R}{\sin \beta }=\frac{r^{\prime }}{\sin \left({\alpha }^{\prime }-\frac{\pi }{2}\right)}=\frac{r^{\prime }}{-{\mathrm{cos\alpha }}^{\prime }}.& \left(\mathrm{Eqn}.11\right)\end{array}$ Here r′=∥r′∥. Using Eqn. 11, Eqn. 9 can be transformed to:

$\begin{array}{cc}\frac{2{\pi R\mathrm{cos\alpha }}^{\prime }}{N}=\frac{-2{\pi r}^{\prime }\mathrm{sin\beta }}{N}.& \left(\mathrm{Eqn}.12\right)\end{array}$ Combining Inequality in Eqn. 10 and Eqn. 12, the following is obtained:

$\begin{array}{cc}{r}^{\prime }<\frac{{\mathrm{N\lambda }}_c}{4\pi \left|\mathrm{sin\beta }\right|},& \left(\mathrm{Eqn}.13\right)\end{array}$ which must be satisfied for any β∈[0,2π). When

$\beta =\frac{\pi }{2}\mathrm{or}\frac{3\pi }{2},$ Eqn. 13 leads to the smallest upper limit of r′:

$\begin{array}{cc}{r}^{\prime }<\frac{{N\lambda }_c}{4\pi },& \left(\mathrm{Eqn}.14\right)\end{array}$ and the length of the local step size maximizes to

$\frac{2{\pi r}^{\prime }}{N}.$ If r′ satisfies Eqn. 14, the length of the local step size of ∥r−r′∥ given by Eqn. 12 also satisfies the Nyquist criterion:

$\begin{array}{cc}\frac{2\pi R\left|{\cos \alpha }^{\prime }\right|}{N}=\frac{2{\pi r}^{\prime }\left|\mathrm{sin\beta }\right|}{N}\le \frac{2{\pi r}^{\prime }}{N}<\frac{{\lambda }_c}{2}.& \left(\mathrm{Eqn}.15\right)\end{array}$

The one-way Nyquist zone S₁ may be defined as:

$\begin{array}{cc}{S}_1=\left\{r^{\prime }||r^{\prime }||<\frac{{N\lambda }_c}{4\pi }\right\}.& \left(\mathrm{Eqn}.16\right)\end{array}$ The boundary of one-way Nyquist zone S₁ can be considered a virtual detection surface, where the sampling spacing is scaled down from the actual detection spacing by

$\frac{R_1}{R}$ with R₁ being the radius of one-way Nyquist zone S₁. For any source point inside zone S₁, there is no spatial aliasing due to spatial sampling because the sampling spacing of the transducer array is less than half of the lower cutoff wavelength.

Next, spatial aliasing due to image reconstruction is analyzed. Eqn. 5 is substituted into Eqn. 6, obtaining:

$\begin{array}{cc}{\hat{p}}_0\left(r^{″}\right)\approx \frac{1}{2{\pi c}^2}\sum _{n=1}^N{w}_n\sum _{m=1}^M{v}_m\frac{p_0\left(r_m^{\prime }\right)}{||r_m^{\prime }-r_n||}\left(1-\left(t+\frac{||r_m^{\prime }-r_n||}{c}\right)\frac{\partial }{\partial t}\right)h_e^{\prime }\left(\frac{||r^{″}-r_n||}{c}-\frac{||r_m^{\prime }-r_n||}{c}\right).& \left(\mathrm{Eqn}.17\right)\end{array}$ Using a differential operation, the following is obtained:

$\begin{array}{cc}\left(1-\left(t+\frac{||r_m^{\prime }-r_n||}{c}\right)\frac{\partial }{\partial t}\right)=h_e^{\prime }\left(t\right)-\left(t+\frac{||r_m^{\prime }-r_n||}{c}\right)h_e^{″}\left(t\right).& \left(\mathrm{Eqn}.18\right)\end{array}$

The expression

$\left(1-\left(t+\frac{||r_m^{\prime }-r_n||}{c}\right)\frac{\partial }{\partial t}\right)h_e^{\prime }\left(\frac{||r^{″}-r_n||}{c}-\frac{||r_m^{\prime }-r_n||}{c}\right)$ in Eqn. 17 includes two terms:

$h_e^{\prime }\left(\frac{||r^{″}-r_n||}{c}-\frac{||r_m^{\prime }-r_n||}{c}\right)\mathrm{and}\frac{||r^{″}-r_n||}{c}{h}_e^{″}\left(\frac{||r^{″}-r_n||}{c}-\frac{||r_m^{\prime }-r_n||}{c}\right).$ As provided in Section VI, the difference between the spectra of

$\frac{||r^{″}-r_n||}{c}{h}_e^{″}\left(\frac{||r^{″}-r_n||}{c}-\frac{||r_m^{\prime }-r_n||}{c}\right)$ and

$h_e^{″}\left(\frac{||r^{″}-r_n||}{c}-\frac{||r_m^{\prime }-r_n||}{c}\right)$ is negligible. In addition, as shown in Section VII, h_e″(t) has the same upper cutoff frequency as h_e′(t). Thus, if h_e′(t) has an upper cutoff frequency, ƒ_c,

$\left(1-\left(t+\frac{||r_m^{\prime }-r_n||}{c}\right)\frac{\partial }{\partial t}\right)h_e^{\prime }\left(t\right)$ has approximately the same upper cutoff frequency.

In FIG. 1B, extensions of the line segments r−r″, r−r₁′, and r−r′ form angles α″, α₁″ and α′, respectively, with the tangential dotted line that is perpendicular to r. Vectors r″ and r′ form an angle

$\gamma =\arccos \left(\frac{r^{″}}{R}\right)+\arccos \left(\frac{r^{\prime }}{R}\right),$ where r″=∥r″∥ and r′=∥r′∥. Points r″ and r₁″ are on the same circle centered at r. Given a reconstruction location r″ and a source point location r′ as shown in FIG. 1B, the sampling step size of ∥r−r∥−∥r′−r∥ needs only be analyzed while the detection element location r varies. The lengths of the step sizes of both ∥r″−r∥ and ∥r′r∥ reach maxima when r″−r and r′−r are perpendicular to r″ and r′, respectively. If the angle γ between vectors r″ and r′ satisfies

$\gamma =\arccos \left(\frac{r^{″}}{R}\right)+\arccos \left(\frac{r^{\prime }}{R}\right),$ where r″=∥r″∥ and r′=∥r′∥, then the lengths of the step sizes of ∥r′−r∥ and ∥r′-∥ achieve maxima of

$\frac{2{\pi r}^{″}}{N}\mathrm{and}\frac{2{\pi r}^{\prime }}{N},$ respectively, with r at the same location, as shown in FIG. 1B. In addition, as r passes this location clockwise, ∥r″−r∥ increases while ∥r′−r∥ decreases. Thus, the length of the step size of ∥r″−r∥-∥r′−r∥ achieves its maximum of:

$\begin{array}{cc}\frac{2{\pi r}^{″}}{N}+\frac{2{\pi r}^{\prime }}{N}=\frac{2\pi \left(r^{″}+r^{\prime }\right)}{N}.& \left(\mathrm{Eqn}.19\right)\end{array}$ The Nyquist sampling criterion provides that:

$\begin{array}{cc}\frac{2\pi \left(r^{″}+r^{\prime }\right)}{N}<\frac{{\lambda }_c}{2},& \left(\mathrm{Eqn}.20\right)\end{array}$ which is equivalent to:

$\begin{array}{cc}{r}^{″}+r^{\prime }<\frac{N{\lambda }_c}{4\pi }.& \left(\mathrm{Eqn}.21\right)\end{array}$

The physical propagation of the photoacoustic wave from the object being imaged to the transducer detector elements is succeeded by a time-reversal propagation for the image reconstruction. The combined region encompasses a disc with a radius of r″+r′. On the perimeter of this disc, the Nyquist sampling criterion provides that the sampling spacing be less than half of the lower cutoff wavelength, i.e., Eqn. 20.

The two-way Nyquist zone S₂ may be defined as:

$\begin{array}{cc}{S}_2=\left\{r^{\prime }|r^{\prime }<\frac{N{\lambda }_c}{8\pi }\right\}.& \left(\mathrm{Eqn}.22\right)\end{array}$ The boundary of two-way Nyquist zone S₂ may be considered a virtual detection surface where the sampling spacing is scaled down from the actual detection spacing by

$\frac{R_2}{R}$ with R₂ being the radius of two-way Nyquist zone S₂. If the source points are inside zone S₂ and there is image reconstruction at points within zone S₂, then

$r^{\prime }<\frac{N{\lambda }_c}{8\pi }\mathrm{and}{r}^{″}<\frac{N{\lambda }_c}{8\pi },$ respectively, and the inequality of Eqn. 21 is satisfied. Thus, there is no spatial aliasing in image reconstruction, and zone S₂ can be considered an aliasing-free region.

Due to the finite duration of the ultrasonic transducer's temporal response, the function

$\left(1-\left(t+\frac{r_m^{\prime }-r_n}{c}\right)\frac{\partial }{\partial t}\right)h_e^{\prime }\left(t\right)$ has a nonzero value only when t is within a finite interval, denoted as T_e. The broader the bandwidth of h_e(t), the shorter T_e. When

$\frac{r^{″}-r}{c}-\frac{r^{\prime }-r}{c}$ is out of T_e, signals from source point r′ that are detected by the element at r have no contribution to the reconstruction at r″.

Assuming that

$\frac{r^{″}-r}{c}-\frac{r^{\prime }-r}{c}$ belongs to T_e, and with source points inside S₂, there may still be aliasing for reconstruction locations in a region outside zone S₂, but inside zone S₁. For example, where both r″ and r′ are on the boundary of S₂, and the length of the step size of ∥r−r∥-∥r′−r∥ achieves the maximum value of

$\frac{2\pi R\left(\cos {\alpha }^{″}-\cos {\alpha }^{\prime }\right)}{N}=\frac{{\lambda }_c}{2},$ there may be aliasing outside zone S₂, but inside zone S₁. In this example, α″ and α′ denote the angles formed by the line segments r−r″ and r−r′, respectively, with the tangential dotted line that is perpendicular to r.

As shown in FIG. 1B, the reconstruction location r″ is moved to a new position r₁″ outside S₂, but inside S₁. The distance ∥r₁″−r∥=∥r″−r∥ is kept constant and thus,

$\frac{r_1^{″}-r}{c}-\frac{r^{\prime }-r}{c}=\frac{r^{″}-r}{c}-\frac{r^{\prime }-r}{c}$ still belongs to T_e. As r″ moves to r₁″, the angle α″ decreases to α₁″. Both α″ and α₁″ belong to

$\left(0,\frac{\pi }{2}\right),$ then cos α₁″: >cos α″. Thus, for the local step size of ∥r₁″−r∥-∥t″−r∥, there is the estimation

$\frac{2\pi R\left(\cos {\alpha }_1^{″}-\cos {\alpha }^{\prime }\right)}{N}>\frac{2\pi R\left(\cos {\alpha }^{″}-\cos {\alpha }^{\prime }\right)}{N}=\frac{{\lambda }_c}{2},$ which means that spatial aliasing appears in image reconstruction. Switching the source point and reconstruction locations, the analysis is repeated and a similar conclusion is drawn: with source points inside zone S₁, but outside zone S₂, spatial aliasing may occur when doing reconstruction.

FIG. 1C illustrates an example of the relative sizes of the three zones: (1) the two-way Nyquist zone S₂, (2) one-way Nyquist zone, S₁, with no spatial aliasing due to spatial sampling, and (3) detection zone, S₀, with spatial aliasing due to spatial sampling and image reconstruction. The zones S₀, S₁, and S₂ can be defined by Eqn. 9, Eqn. 16, and Eqn. 22 respectively. Spatial aliasing due to spatial sampling do not appear for source points and reconstruction locations inside one-way Nyquist zone S₁ (e.g., zone S₁, 102, FIG. 1C), but do appear for objects outside one-way Nyquist zone, S₁. Spatial aliasing in image reconstruction does not appear for source points of objects and reconstruction locations inside two-way Nyquist zone S₂, (e.g., zone S₂, 102, in FIG. 1C), but appears for other combinations of objects and image reconstruction locations.

FIGS. 2A, 2B, and 2C are schematic diagrams of different combinations of source points and image reconstruction locations in the imaging domain of a full-ring ultrasonic array that are subject to no spatial aliasing or spatial aliasing due to spatial sampling and/or image reconstruction, according to certain aspects. The diagrams include a two-way Nyquist zone S₂ 201, an initial one-way Nyquist zone S₁ 202, and an imaging domain S₀ 204. The outer boundary of the imaging domain S₀ 204 represents the detection surface of the full-ring transducer array having a radius R and having N evenly-distributed transducer elements. The center O of the imaging domain S₀ 204 represents the origin of a coordinate system used in image reconstruction. The first line radiating from the origin O represents the range of source locations for spatial sampling, while the second line radiating from the tip of the first line represents the range of reconstruction locations for image reconstruction. A solid line means no aliasing, while a dotted line means aliasing. The boundaries for the zones S₀, S₁, and S₂ may be defined by Eqn. 9, Eqn. 16, and Eqn. 22 respectively.

FIG. 2A is a schematic diagram depicting spatial aliasing appearing during image reconstruction (e.g., image reconstruction that uses universal back projection) in a PACT system with a full-ring ultrasonic transducer array, according to one implementation. FIG. 2A illustrates three zones (regions): (1) an innermost two-way Nyquist zone, S₂, which is an aliasing-free region (2) a one-way Nyquist zone, S₁, does not have spatial aliasing due to spatial sampling, and (3) a detection zone, S₀, that has spatial aliasing in spatial sampling and in image reconstruction.

FIG. 2B is a schematic diagram representing mitigating spatial aliasing in image reconstruction using spatial interpolation, according to an implementation. The large arrow enlarging two-way Nyquist zone S₂ 201 to zone S₁ 202 depicts the mitigation resulting from spatial interpolation. In this example, spatial interpolation is shown to eliminate spatial aliasing in the three instances of image reconstruction, extending the two-way Nyquist zone S₂ to the one-way Nyquist zone S₁ to provide a larger aliasing-free region. This mitigation of spatial aliasing in image reconstruction is represented by the change from the dotted lines in the three instances in FIG. 2A to the solid lines in FIG. 2B.

FIG. 2C is a schematic diagram representing mitigating spatial aliasing in spatial sampling using location-dependent temporal filtering, according to an implementation. The large arrow depicts the extension of zone S₁ 202 to zone S₁′ 203 due to mitigation based on location-dependent temporal filtering. Location-dependent temporal filtering extends (adjusts) the initial one-way Nyquist zone S₁ to a larger adjusted one-way Nyquist zone S₁′. By further applying spatial interpolation to this example, the zone S₁′ becomes an aliasing-free region. As can be seen, zone S₁′ 203 is not extended all the way to edge zone S₀ 204 in order to maintain high resolution imaging.

In certain implementations, location-dependent temporal filtering extends (adjusts) an initial one-way Nyquist zone S₁ to a larger adjusted one-way Nyquist zone S₁′ that does not cover the area of the entire two-way Nyquist zone S₀. In one example, the one-way Nyquist zone S₁′ covers an area that is more than 95% of the area of the two-way Nyquist zone S₀. In one example, the one-way Nyquist zone S₁′ covers an area that is more than 90% of the area of the two-way Nyquist zone S₀. In one example, the one-way Nyquist zone S₁′ covers an area that is more than 80% of the area of the two-way Nyquist zone S₀. In one example, the one-way Nyquist zone S₁′ covers an area that is in a range from about 80% to 100% of the area of the two-way Nyquist zone S₀.

C. PACT Systems and Methods with Spatiotemporal Antialiasing

In certain aspects, PACT systems and methods with spatiotemporal antialiasing mitigate (substantially reduce or eliminate) spatial aliasing in image reconstruction using spatial interpolation and mitigate spatial aliasing in spatial sampling using a location-dependent temporal filtering technique. The location-dependent temporal filtering is designed based on the geometry of the ultrasonic transducer array being implemented. For example, in one implementation, it is determined whether a distance between the reconstruction location and the transducer array is within or outside a first distance or radius from the transducer array in a one-way Nyquist zone where spatial aliasing does not occur and then only apply location-based temporal filtering if the reconstruction location is outside the one-way Nyquist zone. According to one aspect, the location-dependent temporal filtering technique is designed to both maintain high resolution of the images being reconstructed while substantially mitigating spatial aliasing due to undersampling.

Spatial aliasing solely in image reconstruction, but not in spatial sampling, can be mitigated using spatial interpolation. For example, without spatial aliasing, spatially continuous photoacoustic signals can be accurately recovered from spatially discrete photoacoustic signals (e.g., digitized signals from one or more DAQs) through, e.g., Whittaker-Shannon spatial interpolation. In one aspect, if there is no spatial aliasing in spatial sampling, a PACT method implements spatial interpolation and no spatial aliasing occurs in reconstructing the image using spatially continuous photoacoustic signals.

1. Spatiotemporal Antialiasing for Circular Transducer Arrays (e.g., Full-Ring Array)

Although this section uses a circular transducer array as an example, one or more of the operations discussed in this section may also apply to a linear array or other array geometry. In one aspect, PACT methods with spatiotemporal antialiasing include image reconstruction operations that utilize spatial interpolation with an operation that numerically increases the number of detection elements in the transducer array. At any given time, t, a function sampled at

${\theta }_n=\frac{2\pi n}{N},$ n=0, 1, 2, . . . , N−1 may be defined as:ƒ_R(θ)={circumflex over (p)}(r,t),  (Eqn. 23)

where r=(R cos θ, R sin θ), θ∈[0,2π). Referring to FIG. 1C, for objects inside zone S₁, there is no aliasing in spatial sampling. The function ƒ_R(θ) can be recovered from ƒ_R (θ_n), n=0, 1, 2, . . . , N−1 through spatial interpolation such as, e.g., the Whittaker-Shannon spatial interpolation. To extend the aliasing-free region S₂ to a larger S₂′, the number of detection (transducer) elements may be numerically doubled from N′=2N during a spatial interpolation operation. Substituting N′ for N in Eqn. 22, a larger aliasing-free region S₂′ is obtained:

$\begin{array}{cc}{S}_2^{\prime }=\left\{r^{\prime }|r^{\prime }<\frac{N^{\prime }{\lambda }_c}{8\pi }\right\}=\left\{r^{\prime }|r^{\prime }<\frac{N{\lambda }_c}{4\pi }\right\}=S_1.& \left(\mathrm{Eqn}.24\right)\end{array}$

With source and reconstruction locations inside S₂′, image reconstruction has no aliasing. Since S₁=S₂′, the spatial interpolation operation mitigated spatial aliasing in image reconstruction. For source points outside one-way Nyquist zone, S₁, spatial sampling has aliasing, and spatial interpolation cannot recover the lost information. The spatial aliasing in spatial sampling may be mitigated by temporal lowpass filtering, which expands zone S₁ to a larger region of zone S₁′:

$\begin{array}{cc}{S}_1^{\prime }=\left\{r^{″}\text{|}r^{″}<r^{\prime }\right\}\mathrm{with}{r}^{\prime }>\frac{N{\lambda }_c}{4\pi }.& \left(\mathrm{Eqn}.25\right)\end{array}$

In certain aspects, PACT systems/methods with spatiotemporal antialiasing perform one or more location-dependent temporal filleting operations before performing spatial interpolation and image reconstruction operations. The one or more location-dependent temporal filtering operations process the photoacoustic signals using one or more lowpass filters with an upper cutoff frequency. In one aspect, one or more location-dependent temporal filtering operations apply an antialiasing filter with one or more lowpass filters having an upper cutoff frequency of

$f_c^{\prime }=\frac{\mathrm{Nc}}{4\pi r^{\prime }}$ where N is the number of elements in the transducer array, c is the speed of sound, r′ is the distance between the center of the transducer array and the reconstruction location. By Replacing λ_c in Eqn. 16 with

${\lambda }_c^{\prime }=\frac{c}{f_c^{\prime }},$ the zone S₁ is extended to:

$\begin{array}{cc}{S}_1^{\prime }=\left\{r^{″}|r^{″}<\frac{N{\lambda }_c^{\prime }}{4\pi }\right\}=\left\{r^{″}\text{|}r^{″}<r^{\prime }\right\}.& \left(\mathrm{Eqn}.26\right)\end{array}$

By applying spatial interpolation during image reconstruction, for source points inside S₁′, any points inside S₁′, of an image can be reconstructed without aliasing artifacts. This one-way Nyquist zone S₁′ may be further extended through temporal lowpass filtering. In one example, spatial antialiasing may extend the zone S₁′ depicted in FIG. 2C to the whole region S₀ by applying an antialiasing filter with one or more lowpass filters having an upper cutoff frequency of

$f_c^{\prime }=\frac{\mathrm{Nc}}{4\pi r^{\prime }}.$ Directly extending S₁′ to S₀ may, however, compromise image resolution. Since lowpass filtering removes high-frequency signals (i.e. with frequencies above the upper cutoff frequency), reconstructed images may be blurred if an upper cutoff frequency of

$f_c^{\prime }=\frac{\mathrm{Nc}}{4\pi r^{\prime }}$ is used.

In certain implementations, PACT systems/methods with spatial antialiasing implement location-dependent temporal filtering that balances image resolution with reducing aliasing from spatial sampling to extend the zone S₁. To reconstruct the image at r′∈S₀, these PACT systems/methods including an antialiasing filter with one or more lowpass filters that are designed based at least in part on the distance r′=∥r′∥ of an object of interest being imaged at an imaging plane) to the center of the transducer array.

For example, the PACT methods may include one or more location-dependent temporal filtering operations. During these operations, it is determined whether

$r^{\prime }<\frac{N{\lambda }_c}{4\pi },$ where N is the number of elements in the transducer array, λ_c is the lower cutoff wavelength, and r′ is the distance between the center of the transducer array and the imaging location. In the region where

$r^{\prime }<\frac{N{\lambda }_c}{4\pi },$ there is no spatial aliasing in spatial sampling. If it is determined that

$r^{\prime }<\frac{N{\lambda }_c}{4\pi },$ spatial interpolation is applied and then image reconstruction is performed. If it is determined that

$r^{\prime }\ge \frac{N{\lambda }_c}{4\pi },$ the photoacoustic signals are filtered first with an antialiasing filter having one or more low pass filters with an upper cutoff frequency of

$f_c^{\prime }=\frac{\mathrm{Nc}}{4\pi r^{\prime }},$ and then spatial interpolation is applied, and then image reconstruction. N is the number of elements in the transducer array, c is the speed of sound, r′ is the distance between the center of the transducer array and the reconstruction location.

2. Spatiotemporal antialiasing for Linear Transducer Arrays

In another aspect, a PACT system with spatiotemporal antialiasing includes a linear array with N (e.g., N≥4) transducer elements, a pitch d, and a lower cutoff wavelength λ_c<2d. (e. g., N≥4). FIG. 3 is a schematic illustration of a portion of a surface of a linear ultrasonic transducer array 300 with N elements and pitch d, according to an implementation. The illustration shows a center O of the linear ultrasonic transducer array 300. A Cartesian coordinate is formed with the array's center O as the origin, the array direction as a y-axis, and its normal vector as a z-axis. Two adjacent ransducer element locations r_n, r_n+1, and a source point location r′ are also shown. Two symmetric hyperbolas (with sampling step size

$\frac{{\lambda }_c}{2})$ with z≥0 are shown as solid curves and their asymptotes as dashed lines. The one-way Nyquist zone S₁ is formed by the points above the two asymptotes crossing z-axis.

FIG. 4 is a schematic illustration of the portion of the surface of the linear ultrasonic transducer array 300 in FIG. 3 depicting the imaging domain S₀, the one-way Nyquist zone S₁, and the two-way Nyquist zone S₂ outlined with three lines, respectively. The illustration includes outlined boundaries 302, 301 of the one-way Nyquist zone S₁, and the two-way Nyquist zone S₂, respectively.

Using the linear array center O as the origin, the linear transducer array 300 as one axis and its normal as another axis, a cartesian coordinate can be constructed as shown in FIG. 3. For a source point location r′=(y, z) and two adjacent element locations

$r_n=\left(\left(n-\frac{N-1}{2}\right)d,0\right)\mathrm{and}{r}_{n+1}=\left(\left(n-\frac{N-3}{2}\right)d,0\right),$ n=0, . . . , N−2, the sampling step size can be expressed as:

$\begin{array}{cc}r^{\prime }-r_n-r^{\prime }-r_{n+1}=\sqrt{z^2+{\left(y-\left(n-\frac{N-1}{2}\right)d\right)}^2}-\sqrt{z^2+{\left(y-\left(n-\frac{N-3}{2}\right)d\right)}^2}.& \left(\mathrm{Eqn}.27\right)\end{array}$

The sampling step size is

$\frac{{\lambda }_c}{2}$

$\begin{array}{cc}r^{\prime }-r_n-r^{\prime }-r_{n+1}=\frac{{\lambda }_c}{2},& \left(\mathrm{Eqn}.28\right)\end{array}$ which is a hyperbola with a standard form:

$\begin{array}{cc}\frac{{\left(y-\left(n-\frac{N-2}{2}\right)d\right)}^2}{{\left(\frac{{\lambda }_c}{4}\right)}^2}-\frac{z^2}{{\left(\frac{d}{2}\right)}^2-{\left(\frac{{\lambda }_c}{4}\right)}^2}=1.& \left(\mathrm{Eqn}.29\right)\end{array}$

For any source points (e.g., r_n, r_n+1) that lie in between the two curves of this hyperbola, the spatial Nyquist criterion given by:

$\begin{array}{cc}||r^{\prime }-r_n||-||r^{\prime }-r_{n+1}||<\frac{{\lambda }_c}{2},& \left(\mathrm{Eqn}.30\right)\end{array}$ is satisfied. Based on this geometry, the spatial Nyquist criterion for all transducer element pairs can be simplified into two cases: the leftmost one

$||r^{\prime }-r_0||-||r^{\prime }-r_1||<\frac{{\lambda }_c}{2},$ and the rightmost one

$||r^{\prime }-r_{N-2}||-||r^{\prime }-r_{N-1}||<\frac{{\lambda }_c}{2},$ Illustrations of these two hyperbolas with z≥0 are shown in FIG. 3 as solid line curves and their asymptotes are shown as dashed lines. Due to the symmetry between these the two hyperbolas, they intersect with z axis at the same point

$\left(0,z_1^{\prime }\right)=\left(0,\sqrt{\left({\left(\frac{2d}{{\lambda }_c}\right)}^2-1\right)\left({\left(\frac{N-2}{2}d\right)}^2-{\left(\frac{{\lambda }_c}{4}\right)}^2\right)}\right),$ while their asymptotes intersect with the z axis at:

$\begin{array}{cc}\left(0,z_1\right)=\left(0,\frac{N-2}{2}d\sqrt{{\left(\frac{2d}{{\lambda }_c}\right)}^2-1}\right).& \left(\mathrm{Eqn}.31\right)\end{array}$ The one-way Nyquist zone can be approximated by:

$\begin{array}{cc}{S}_1=\left\{\left(y,z\right)|z>\max \left\{\begin{array}{c}\left|y-\frac{N-2}{2}d\right|\\ \left|y+\frac{N-2}{2}d\right|\end{array},\right\}\sqrt{{\left(\frac{2d}{{\lambda }_c}\right)}^2-1}\right\},& \left(\mathrm{Eqn}.32\right)\end{array}$ which is the region above the two intersecting asymptotes, as shown in FIG. 3. At y=0, the approximation error achieves the maximum value:

$\begin{array}{cc}\left|z_1-z_1^{\prime }\right|\le \frac{{\left(\frac{{\lambda }_c}{4}\right)}^2\sqrt{{\left(\frac{2d}{{\lambda }_c}\right)}^2-1}}{\sqrt{{\left(N-2\right)}^2{d}^2-{\left(\frac{{\lambda }_c}{2}\right)}^2}}<\frac{{\left(\frac{{\lambda }_c}{4}\right)}^2\sqrt{{\left(\frac{2d}{{\lambda }_c}\right)}^2-1}}{d\sqrt{\left(N-1\right)\left(N-3\right)}}.& \left(\mathrm{Eqn}.33\right)\end{array}$ For example, if c=1.5 mm·μs⁻¹, ƒ_c=4.5 MHz, N=256, and d=0.25 mm, z₁≈35.5 mm and

$\frac{2\left|z_1-z_1^{\prime }\right|}{{\lambda }_c}<7.4\times {10}^{-4},$ which proves the accuracy of using S₁ as the one-way Nyquist zone. Further, the two-way Nyquist zone may be approximated using:

$\begin{array}{cc}{S}_2=\left\{\left(y,z\right)|z>\max \left\{\begin{array}{c}\left|y-\frac{N-2}{2}d\right|\\ \left|y+\frac{N-2}{2}d\right|\end{array},\right\}\sqrt{{\left(\frac{4d}{{\lambda }_c}\right)}^2-1}\right\},& \left(\mathrm{Eqn}.34\right)\end{array}$

The PACT system with spatiotemporal antialiasing can implement spatial interpolation to extend S₂ to S₁, and location-dependent temporal filtering can be implemented to further extend S₁.

III. Examples of PACT Systems and Methods with Spatiotemporal Antialiasing

A. Examples of PACT Systems with Spatiotemporal Antialiasing

FIG. 5 is a schematic diagram of components of a PACT system with spatiotemporal antialiasing 500, according to certain implementations. The PACT system with spatiotemporal antialiasing 500 includes one or more light sources 510 (e.g., a pulsed laser) that provide illumination (e.g., generate pulsed or modulated light) and an optical system 520 having one or more optical components configured to propagate the illumination from the one or more light sources to a specimen 530 during operation. At the instant in time shown in FIG. 5, the specimen 530 is at a location at which it can receive the illumination during operation. It would be understood that at other times the specimen 530 may not be at this location as denoted by the dotted line.

The optical system 520 includes one or more optical components (e.g., lens(es), optical filter(s), mirror(s), beam steering device(s), beam-splitter(s), optical fiber(s), relay(s), and/or beam combiner(s)) configured to propagate and/or alter light from the one or more light sources 510 to provide the illumination to the specimen 530. The optical system 520 is in optical communication with the one or more light sources 510 to receive illumination during acquisition. In one aspect, the optical system 520 is configured to convert a light beam from the one or more light sources 510 into shaped illumination such as donut-shaped illumination (sometimes referred to herein as “donut illumination” or a “donut beam”). Donut-shaped illumination and other circular illumination may be generated by employing an engineered diffuser such as, e.g., an EDC15 diffuser made by RPC Photonics®.

In some implementations, a PACT system with spatiotemporal antialiasing includes one or more light sources (e.g., a pulsed laser) that can generate pulsed or modulated illumination such as, e.g., pulsed or modulated light at a near-infrared wavelength or a narrow band of near-infrared wavelengths. For example, a light source may be a pulsed laser that can generate near infrared pulses having a wavelength or narrow band of wavelengths in a range from about 700 nm to about 1000 nm. As another example, a light source may be a pulsed laser that can generate near infrared pulses having a wavelength or narrow band of wavelengths in a range from about 600 nm to about 1100 nm. In yet another example, a light source may be a pulsed laser that can generate near infrared pulses with a wavelength or narrow band of wavelengths greater than 760 nm. In yet another example, a light source may be a pulsed laser that can generate near infrared pulses with a wavelength or narrow band of wavelengths greater than 1000 nm. In another example, a light source may be a pulsed laser that can generate a 1064-nm laser beam. A commercially-available example of a suitable pulsed laser is the PRO-350-10, Quanta-Ray® laser with a 10-Hz pulse repetition rate and 8 ns-12 ns pulse width sold by Spectra-Physics®. The low optical attenuation of 1064 nm light or other near infrared light can be used to deeply penetrate to, e.g., a depth of 4 cm, into biological tissues. Imaging of biological tissues using near infrared light is discussed in Smith, A. M., Mancini, M. C. & Nie, S., “Bioimaging: second window for in vivo imaging,” Nat. Nanotechnol. 4, 710-711 (2009), which is hereby incorporated by reference in its entirety. Alternatively, a light source may be a continuous wave laser source that is chopped, modulated and/or gated. In implementations that include a light source in the form of a pulsed laser, the pulse repetition rate may be about 10-Hz in some cases, about 20-Hz in other cases, about 50-Hz in other cases, and about 100-Hz in other cases. In another case, the pulse repetition rate is in a range from about 10-Hz to about 100-Hz.

In one aspect, the one or more light sources of a PACT system with spatiotemporal antialiasing include a tunable narrow-band pulsed laser such as, e.g., one of a quantum cascade laser, an interband cascade laser, an optical parametric oscillator, or other pulsed laser that can be tuned to different narrow bands (e.g., a near-infrared band). In another aspect, the one or more light sources may include a pulsed laser of a single wavelength or approximately a single wavelength. In another aspect, the one or more light sources may include multiple lasers of the same type. For example, multiple lasers of the same type may be used where each of the lasers has a lower power. In another aspect, the one or more light sources may include a combination of different types of lasers. For example, an optical parametric oscillator combined with an Nd:YAG laser may be used in one implementation.

In one aspect, the optical system is configured to a light beam from the one or more light sources into shaped illumination, e.g., a donut beam that may be used to circumferentially illuminate the specimen. For example, the optical system may include an axicon lens (e.g., an axicon lens having 25 mm diameter and a 160° apex angle) followed by an engineered diffuser (e.g., EDC-10-A-2s made by RPC Photonics) to convert a light beam from the one or more light sources into a donut beam. The axicon lens may be positioned to receive a single laser beam propagated from a pulsed laser source and may be configured to convert the beam into a ring having a thickness and diameter and the engineered diffuser may be configured to expand the ring into a donut beam. The donut beam may be used to provide mass energy in homogenized, uniform illumination deep into tissue. An example of a technique for generating donut-shaped illumination can be found in U.S. patent application Ser. No. 16/464,958, titled “SINGLE-IMPULSE PANORAMIC PHOTOACOUSTIC COMPUTED TOMOGRAPHY (SIP-PACT),” and filed on Nov. 29, 2017, which is hereby incorporated by reference in its entirety.

Returning to FIG. 5, the PACT system with spatiotemporal antialiasing 500 also includes an ultrasonic transducer array 540 for sampling photoacoustic waves, an optional mechanism 570 for moving (translating and/or rotating) the ultrasonic transducer array 540, one or more pre-amplifiers 550 for boosting (increasing amplitude) photoacoustic signals communicated from the ultrasonic transducer array 540, one or more data acquisition systems (DAQs) 560 for digitizing the one or more boosted photoacoustic signal, and a computing device 580 for receiving the digitized photoacoustic data from the DAQ(s) 560. The computing device 580 includes an antialiasing filter 581 including one or more lowpass filters that can perform location-dependent temporal filtering of the digitized photoacoustic data communicated from the DAQ(s) 560. The computing device 580 also includes one or more processors and/or other circuitry 582 in electrical communication with the anti-aliasing filter 581 to receive photoacoustic data and configured to perform operations such as spatial interpolation and image reconstruction, an optional (denoted by dotted line) display 586 in electrical communication with the one or more processors and/or other circuitry 582, and a computer readable media 584. The computing device 580 may also communicate one or more control signals to the antialiasing filter 581 to control the upper cutoff frequency or to trigger temporally filtering. For example, the computing device 580 may execute instructions with logic that determine whether the distance between the center of the ultrasonic transducer array and the reconstruction location, r′ is less than

$\frac{N{\lambda }_c}{4\pi }$ where N is the number of elements in the transducer array

$\left(r^{\prime }<\frac{N{\lambda }_c}{4\pi }\right),$ λ_c is the lower cutoff wavelength, and r′ is the distance between the center of the transducer array and the reconstruction location. If it is determined that r′ is less than

$\frac{N{\lambda }_c}{4\pi },$ then spatial interpolation operation(s) are applied and image reconstruction performed without temporal filtering. If it is determined that. If r′ is greater than or equal to

$\frac{N{\lambda }_c}{4\pi },$ then an upper cutoff frequency is determined based on the geometry of the ultrasonic transducer array. For example, the upper cutoff frequency may be calculated as

$f_c^{\prime }=\frac{\mathrm{Nc}}{4{\pi r}^{\prime }},$ where N is the number of elements in the transducer array, c is the speed of sound, r is the distance between the center of the transducer array and the imaging plane. The computing device may then send control signals to trigger the anti-aliasing filter to temporally filter the photoacoustic data. In one example, cutoff frequencies of 3.8 Hz or smaller are used. In other examples, a cutoff frequency in a range between about 1.5 Hz (low frequency transducer used for transcranial human brain imaging) to about 40 MHz (high frequency transducer used for superficial tissue imaging) may be used. In one implementation, a GPU will be used for acceleration of both location-dependent temporal filtering and image reconstruction.

In certain implementations, the PACT system with spatiotemporal antialiasing includes an antialiasing filter (e.g., anti-aliasing filter 581 in FIG. 5) including one or more lowpass filters for location-dependent temporal filtering one or more photoacoustic signals from the DAQ(s). The antialiasing filter may include a single lowpass filter or a combination of the same type or different types of lowpass filters. Some examples of lowpass filters that be implemented include a Butterworth filter, a Cheyshev filter and a Bessel filter. The antialiasing filter is in electrical communication with the one or more DAQs of the PACT system to receive digitized photoacoustic data and remove frequency components with frequencies higher than an upper cutoff frequency ƒ_c from digitized data communicated from the DAQ(s). The antialiasing filter may be a component of a computing device, a component of the controller, a component of one of the one or more DAQs, part of another system component, or a standalone device. Examples of one or more operations that can be used to perform location-dependent temporal filtering are described in Section III(B).

Returning to FIG. 5, the ultrasonic transducer array 540 is coupled or otherwise in acoustic communication with the specimen 530 to be able to receive one or more photoacoustic waves induced by the illumination from the one or more light sources 510 and output and/or record one or more photoacoustic signals. The one or more pre-amplifiers 550 are in electrical communication (directly or via other circuitry) with the ultrasonic transducer array 540 to receive one or more photoacoustic signals. The pre-amplifier(s) 550 can boost one or more photoacoustic signals received and output boosted photoacoustic signals. The DAQ(s) 560 is configured to process the photoacoustic signals, for example, digitize and/or record the digitized photoacoustic data. The DAQ(s) 560 may include at least one digitizer to digitize the signals.

During operation, the ultrasonic transducer array 540 is coupled to or otherwise in acoustic communication with the specimen 530 during signal acquisition. In some cases, an acoustic medium such as an acoustic gel, water, or other medium capable of conveying ultrasound pulses, is provided at least partially between the specimen 530 and the ultrasonic transducer array 540. In other cases, the acoustic medium may be omitted. The ultrasonic transducer array 540 is acoustically coupled to the specimen 530 to be able to detect photoacoustic waves induced by illumination and sample photoacoustic signals. These photoacoustic signals are indicative of the optical absorption of the specimen by the illumination. In one implementation, the PACT system with spatiotemporal antialiasing may also include a tank that is at least partially filed with acoustic medium such as a water tank (e.g., an acrylic water tank). The specimen being imaged may be located directly in the acoustic medium or in a portion of the tank that is submerged or otherwise located in the acoustic medium.

The ultrasonic transducer array 540 includes a plurality of N transducers (sometimes referred to herein as “transducer elements”) operable to collect photoacoustic signals, e.g., in parallel. Each transducer element has an aperture (e.g., a flat-rectangular aperture). The transducer elements have a height, a width, and a pitch. In one case, the pitch is about 1.35 mm. In one case, the width is 0.65 mm. In another case, the pitch is in a range of 1.20 mm to 1.50 mm. In another case, the height is about 5 mm. In another case, the height is in a range of 2 mm to 10 mm. The N transducer elements in the ultrasonic transducer array 540 are arranged in, e.g., a circular array such as a full-ring array, a linear array, or a two-dimensional array. In one implementation, a full-ring ultrasonic transducer array is employed, e.g., to be able to provide panoramic detection. In this case, the full-ring ultrasonic transducer array (e.g., a 512-element full-ring ultrasonic transducer) includes transducer elements distributed along the circumference of a ring having a diameter and an inter-element spacing. The ring diameter may be at least 220 mm in one aspect, may be at least 200 mm in one aspect, or may be at least 250 mm in one aspect. In one aspect, the ring diameter is in a range of about 150 mm to about 400 mm. The inter-element spacing may be less than or equal to about 1.0 mm in one aspect, less than or equal to 0.7 mm in one aspect, less than or equal to 1.5 mm in one aspect, or less than or equal to 2.0 mm in one aspect. In one aspect, the inter-element spacing is in a range of 0 mm to about 5 mm.

In one aspect, the ultrasonic transducer array includes one or more unfocused transducer elements. One or more of the unfocused transducer elements may have a diffraction angle of about 10 degrees in one case, a diffraction angle in a range of about 5 degrees to about 30 degrees in another case, a diffraction angle of about 20 degrees in another case, a diffraction angle in a range of about 5 degrees to about 30 degrees in another case. In one case, each of the unfocused transducer elements has a central frequency in a range of 0.50 MHz to 2.25 MHz and a one-way bandwidth of more than 50%. In another aspect, each of the unfocused transducer elements has a central frequency in a range of 2.25 MHz to 10 MHz and a one-way bandwidth of more than 50%.

Returning to FIG. 5, the PACT system 500 includes an optional (denoted by dashed line) mechanism 570 coupled to or otherwise operably connected to the ultrasonic transducer array 540 to be able to move (translate and/or rotate) the ultrasonic transducer array 540 to one or more positions during operation, e.g., along an axis in one or both directions. For example, the mechanism 570 may be able to move ultrasonic transducer array 540 to two or more different positions along an axis (e.g., between z₁ and z₂ along a z-axis). In addition, the mechanism 570 may be able to hold the ultrasonic transducer array 540 at each position for a period of time. Some examples of time periods that may be used include about 10 seconds, about 15 seconds, and about 20 seconds. In one case, the time period may be in a range from about 10 seconds to about 20 seconds. The step size between positions may be defined by the elevational resolution of the 2d reconstructed image for the PACT implementation being used. For example, if the elevational resolution is 3-4 cm thick, the step size may be 3-4 cm. An example of a commercially-available mechanism for linear movement is a linear stage KR4610D made by THK America, Inc.

The one or more DAQ(s) 560 record photoacoustic signals at time intervals defined by a sampling frequency. In one example, the sampling frequency is in a range from about 4 MHz to about 100-Hz. In another example, the sampling frequency is 40 MHz.

According to one aspect, the one or more DAQs and one or more pre-amplifiers of a PACT system provide one-to-one mapped associations with the transducers in the ultrasonic transducer array. These one-to-one mapped associations allow for fully parallelized data acquisition of all ultrasonic transducer channels and avoids the need for multiplexing after each laser pulse excitation or other modulated or pulsed excitation illumination. With one-to-one mapped associations between pre-amplifiers and transducer elements, each ultrasound transducer element in the array is in electrical communication with one dedicated pre-amplifier channel (also referred to as “preamp channel”). The one dedicated pre-amplifier channel is configured to amplify only photoacoustic signals detected by the one associated/mapped ultrasound transducer. These one-to-one mapped associations between the transducers and the pre-amplifier channels allow for parallelized pre-amplification of the photoacoustic signals detected by the plurality of transducers in the ultrasound transducer array. With one-to-one mapped analog-to-digital sampling, each pre-amplifier is operatively coupled to a corresponding dedicated data channel of an analog-to-digital sampling device in a DAQ to enable parallelized analog-to-digital sampling of the plurality of pre-amplified PA signals. The pre-amplified PA signals produced by each individual preamp channel are received by a single dedicated data channel of the at least one analog-to-digital sampling devices. Any suitable number of pre-amplifier devices and/or DAQ devices may be used to provide the one-to-one mapping. For example, a PACT system may include four 128-channel DAQs (e.g., SonixDAQ made by Ultrasonix Medical ULC with 40 MHz sampling rate, 12-bit dynamic range, and programmable amplification up to 51 dB) in communication with four 128-channel pre-amplifiers to provide simultaneous one-to-one mapped associations with a 512-element transducer array. This PACT system can acquire photoacoustic signals from a cross section within 100 μs without multiplexing after each laser pulse excitation. The plurality of pre-amplifier channels may be directly coupled to the corresponding plurality of ultrasound transducers or may be coupled with electrical connecting cables. In one aspect, wireless communication may be employed.

Each of the one or more pre-amplifiers of a PACT system with spatiotemporal antialiasing may be set to a pre-amplifier gain that may be determined by one or more factors. For example, the pre-amplifier gain may be determined based on one or more of a minimum signal-to-noise ratio and one or more operating parameters of the data acquisition and processing system components such as analog-to-digital sampling devices (digitizers) of the DAQs, signal amplifiers, buffers, and the computing device. In one aspect, the pre-amplifier gain is in a range that is high enough to enable transmission of the photoacoustic signals with minimal signal contamination, but below a gain that may saturate the dynamic ranges of the DAQ system used to digitize the photoacoustic signals amplified by the pre-amplifier(s). In certain aspects, the gain of the plurality of pre-amplifier channels may be at least about 5 dB, at least about 7 dB, at least about 9 dB, at least about 11 dB, at least about 13 dB, at least about 15 dB, at least about 17 dB, at least about 19 dB, at least about 21 dB, at least about 23 dB, at least about 25 dB, or at least about 30 dB.

According to one aspect, a PACT system with spatiotemporal antialiasing includes a plurality of multi-channel preamplifiers (e.g., 128-channel preamplifiers) in electrical communication with an N-element ultrasonic transducer array and a plurality of multi-channel data acquisition systems (DAQs) in electrical communication with the plurality of pre-amplifiers respectively. Each of the DAQs is in communication with one of the preamplifiers. The plurality of multi-channel preamplifiers and the plurality of DAQs may also be in one-to-one mapping association with the N element ultrasonic transducer array. For example, a set of four 128-channel preamplifiers may be in electrical communication with a 512-element ultrasonic transducer array and a set of four 128-channel data acquisition systems (DAQs) may be in electrical communication with four 128-channel pre-amplifier(s) respectively. The four sets of 128-channel DAQs provide simultaneous one-to-one mapped associations with the 512-element ultrasonic transducer array to enable acquiring photoacoustic signals from a cross section without multiplexing after each laser pulse excitation.

Returning to FIG. 5, the PACT system with spatiotemporal antialiasing 500 also includes a computing device 580 having one or more processors or other circuitry 582, an optional (denoted by dashed line) display 586 in electrical communication with the processor(s) 582, and a computer readable media (CRM) 584 in electronic communication with the processor(s) or other circuitry 582. The computer readable media (CRM) 584 may be, e.g., a non-transitory computer readable media. Optionally (denoted by dashed line), the computing device 580 is in electronic communication with the light source(s) 510 to send control signals to trigger the illumination (e.g., triggering laser pulses). The computing device 580 is in electrical communication with the one or more DAQs 560 to receive data transmissions and/or to send control signal(s). Optionally (denoted by dashed line), the computing device 580 is in electronic communication with the one or more pre-amplifiers 550 to send control signal(s), e.g., to adjust amplification. The computing device 580 is in electrical communication with the antialiasing filter 581 to send control signals to adjust the cutoff frequency. The electrical communication between system components of the PACT system 500 may be in wired and/or wireless form. One or more of the electrical communications between components of the PACT system 500 may be able to provide power in addition to communicate signals. The computing device 580 may be, for example, a personal computer, an embedded computer, a single board computer (e.g. Raspberry Pi or similar), a portable computation device (e.g. tablet), a controller, or any other computation device or system of devices capable of performing the functions described herein. The computing device 580 is in electronic communication with the mechanism 570 to send control signals to control the movement and/or hold positions of the ultrasonic transducer array 540. The processor(s) 582 are in electrical communication with the CRM 584 to store and/or retrieve data such as the photoacoustic signal data. The one or more processor(s) and/or other circuitry 582 are in electrical communication with the optional display 586 to display data. Although not shown, the computing device 580 may also include a user input component for receiving data from a user.

The one or more processors and/or other circuitry 582 execute instructions stored on the CRM 584 to perform one or more operations of PACT methods. In certain implementations, the processor(s) 582 and/or other circuitry 582 execute instructions to perform one or more of: 1) determining an upper cutoff frequency for the antialiasing filter to trigger the antialiasing filter to comment temporal filtering; 2) communicating control signals to one or more components of the PACT system with spatiotemporal antialiasing 500, and 3) performing one or more reconstruction operations to reconstruct a two-dimensional or three-dimensional photoacoustic image of the specimen 530 using photoacoustic data. For example, the processor(s) and/or other circuitry 582 and/or one or more external processors may execute instructions that communicate control signals to the mechanism 570 to scan the ultrasonic transducer array 540 along a z-axis between to two elevations (3D mode) or move the ultrasonic transducer array 540 to one or more different elevations (2D mode) and send control signals to the digitizer in the DAQ(s) 560 to simultaneously record photoacoustic signals received by ultrasonic transducer array 540 from the illuminated regions of the specimen.

In some implementations, the PACT system includes one or more communication interfaces (e.g., a universal serial bus (USB) interface). Communication interfaces can be used, for example, to connect various peripherals and input/output (I/O) devices such as a wired keyboard or mouse or to connect a dongle for use in wirelessly connecting various wireless-enabled peripherals. Such additional interfaces also can include serial interfaces such as, for example, an interface to connect to a ribbon cable. It should also be appreciated that the various system components can be electrically coupled to communicate with various components over one or more of a variety of suitable interfaces and cables such as, for example, USB interfaces and cables, ribbon cables, Ethernet cables, among other suitable interfaces and cables.

FIG. 6 is signal flow diagram of signals communicated between components of a PACT system with spatiotemporal antialiasing 600, according to implementations. The PACT system 600 includes one or more light sources 610 (e.g., a pulsed laser) and an optical system (not shown) that is configured to propagate/alter illumination from the one or more light sources 610 to a specimen being imaged. The PACT system 600 also includes an ultrasonic transducer array 640 that is coupled to or otherwise in acoustic communication with a specimen to receive photoacoustic signals induced by the illumination. The PACT system 600 also includes a optional (denoted by dotted line) mechanism 670 (e.g., a linear scanner) coupled to or otherwise operably connected to the ultrasonic transducer array 640 that can move the ultrasonic transducer array 640 to one or more elevational positions and/or scan the ultrasonic transducer array 640 between two elevational positions. The PACT system 600 also includes one or more pre-amplifiers 650 and one or more data acquisition systems (DAQs) 660 in one-to-one mapped association with the transducers in the ultrasonic transducer array 640. The one or more pre-amplifiers 650 are in electrical communication with the ultrasonic transducer array 640 to receive one or more photoacoustic signals and the DAQ(s) 660 are in electrical communication with the pre-amplifier(s) 650 to receive one or more boosted photoacoustic signals.

The PACT system with spatiotemporal antialiasing 600 also includes a computing device 680 having one or more processors or other circuitry and a computer readable media (CRM) in electronic communication with the processor(s). The PACT system 600 includes an anti-aliasing filter that may reside on the computer readable media. The antialiasing filter 681 is in electrical communication with the DAQ(s) 650 to receive photoacoustic data. The computing device 680 includes instructions residing on the computer readable media that can be executed to perform functions of the system 600 such as image reconstruction, spatial interpolation and determining an upper cutoff frequency of the anti-aliasing filter 681 and determining when to trigger temporal filtering.

The PACT system 600 also includes a controller 685 in electronic communication with the DAQ(s) 660 and the mechanism 670 to send control signals. The controller 685 may include one or more processors. To synchronize components of the PACT system 600, the one or more light sources 610 are configured to transmit a trigger signal to the controller 684 that triggers transmission of control signals to the DAQ(s) 660 to record signals and the mechanism 670 to move the ultrasonic transducer array 640. The computing device 680 is in communication with the controller 685 to be able to transmit control signals. The electrical communication between system components of the PACT system 600 may be in wired and/or wireless form. The electrical communications may be able to provide power in addition to communicate signals in some cases. During operation, digitized radio frequency data from the antialiasing filter 681 may be first stored in an onboard buffer, and then transferred to the computing device 680, e.g., through a universal serial bus 2.0. The DAQ(s) 660 may be configured to record photoacoustic signals within a time period, e.g., 100 μs, after each laser pulse excitation.

B. Examples of PACT Methods with Spatiotemporal Antialiasing

The PACT methods described in this section can be used to obtain one or more two-dimensional images and/or one or more three-dimensional volumetric images. The operations of these PACT methods can be performed by one or more components of a the PACT system with spatiotemporal antialiasing such as, e.g., the PACT system with spatiotemporal antialiasing 500 shown in FIG. 5 or the PACT system with spatiotemporal antialiasing 600 shown in FIG. 6. One or more of the operations may be performed executing instructions retrieved from computer readable media.

FIG. 7 is a flowchart depicting operations of a PACT method with spatiotemporal antialiasing that applies one or more antialiasing techniques (e.g., spatial interpolation and location-based temporal filtering), according to various implementations. Some examples of antialiasing techniques include one or both of spatial interpolation and location-based temporal filtering. For example, one implementation is directed to a method that applies spatial interpolation, but not location-based temporal filtering, and implements universal back-propagation for image reconstruction. As another example, one implementation is directed to a method that applies both spatial interpolation and location-based temporal filtering, and then applies universal back-propagation during image reconstruction. Spatial interpolation may be applied after location-based on temporal filtering in this example.

At operation 710, photoacoustic data is acquired. For example, one or more control signals may be sent to component(s) of the PACT system to perform data acquisition operations or recorded photoacoustic data may be retrieved from a computer readable media such as an onboard buffer.

In one aspect, the PACT system may control operations of system components to perform data acquisition. During the data acquisition phase, the recording of photoacoustic data by the one or more data acquisition systems (DAQs) may be synchronized with light pulses from the one or more light sources. To synchronize data acquisition with light pulses from the light source(s), the external trigger from the light source(s) may be used to trigger recording by the data acquisition systems and/or trigger movement of any mechanism that may be moving the transducer array according to one aspect. For example, during data acquisition in an exemplary two-dimensional (2D) mode, the ultrasonic transducer array may be moved to one or more elevational positions (e.g., different locations z₁, z₂, . . . along a z-axis of a linear stage) and held in each elevational position for a period of time during which the one or more DAQs record data. Some examples of time periods that may be used include: about 10 seconds, about 15 seconds, and about 20 seconds. In another example, the time period is in a range of about 10 seconds to about 20 seconds. At each cross section, photoacoustic data is continuously recorded at a certain sampling rate to monitor the cross section. As another example, during data acquisition in an exemplary three-dimensional (3D) mode, the ultrasonic transducer array is scanned through multiple scanning steps (elevational positions) through a depth.

During an exemplary data acquisition phase, digitized radio frequency data from the DAQs may be stored in an onboard buffer, and then transmitted to the computing device through e.g., a universal serial bus 2.0 or a universal serial bus 3.0. The photoacoustic data is recorded by the one or more data acquisition systems at a sampling frequency. In one aspect, the sampling frequency is 40 MHz. In another aspect, the sampling frequency may be in a range from 4 MHz to 100 MHz. The one or more data acquisition systems may be set to record photoacoustic data within a particular time period (e.g., 100 μs, 200 μs, or 300 μs) after each illumination e.g., laser pulse excitation. In certain implementations, a PACT system with spatiotemporal antialiasing is equipped with a one-to-one mapped signal amplification and data acquisition (DAQ) systems or DAQ circuits to the transducer elements. If there is one-to-one mapping associations, the one or more data acquisition systems can acquire photoacoustic signals from a cross-section of a specimen without multiplexing after illumination e.g., after each laser pulse excitation.

Returning to FIG. 7, at operation 720, one or more antialiasing techniques are implemented. For example, one or both of (i) location-based temporal filtering, and (ii) spatial interpolation may be implemented. If both spatial interpolation and location-based temporal filtering are implemented, location-dependent temporal filtering is applied before spatial interpolation according to one aspect. In some implementations, one or both of (i) location-based temporal filtering, and (ii) spatial interpolation may be used with UBP reconstruction.

An example of a spatial interpolation technique that may be used is the frequency domain implementation of Whittaker-Shannon interpolation. Given a time t, signals from all the detection elements formed a vector with length N. A fast Fourier transformation (FFT) is applied to the vector and zeros were padded behind the highest frequency components to double the vector length to 2N in a new vector. The inverse FFT is then applied to the new vector to spatially interpolate the data. Varying time t, one can generate a new set of signals from the detected signals. In this example, spatial interpolation numerically doubled the number of detection elements.

In certain implementations, the PACT methods implements location-based temporal filtering by applying an antialiasing filter with one or more low pass filters to remove components with frequencies higher than an upper cutoff frequency ƒ_c′. The antialiasing filter may include a single lowpass filter or a combination of the same type or different types of lowpass filters. Some examples of lowpass filters that be implemented include a Butterworth filter, a Cheyshev filter, and a Bessel filter. In one example, the antialiasing filter includes at least in part a third-order lowpass Butterworth filter and a sinc filter (with the same cutoff frequency). The PACT method determines an upper cutoff frequency based on the distance between the reconstruction location and the transducer array being used in the PACT system. If the transducer array has a circular geometry, the upper cutoff frequency may be calculated as

$f_c^{\prime }=\frac{\mathrm{Nc}}{4{\pi r}^{\prime }}$ where N is the number of elements in the transducer array, c is the speed of sound, is the distance between the reconstruction location and the center of the transducer array. If the transducer array has a linear geometry, the upper cutoff frequency may be calculated using

$f_c^{\prime }=\frac{c}{2d}\sqrt{1+{\left(\frac{z}{\max \left\{\begin{array}{c}y-\frac{N-2}{2}d,\\ y+\frac{N-2}{2}d\end{array}\right\}}\right)}^2},$ where d is the pitch of the transducer array, c is the speed of sound, and ƒ_c′ is the upper cutoff frequency. Temporal filtering may be applied to every pixel/voxel location. In one implementation, location-based temporal filtering is applied to one or more groups of pixels/voxels with similar cutoff frequencies. The signals for these groups of pixels/voxels are filtered with different levels of cutoff frequencies before image reconstruction. During reconstruction, location-dependent cutoff frequencies are calculated, and the filtered signals are retrieved directly.

At operation 730, the PACT system performs image reconstruction. Image reconstruction may include (i) reconstructing one or more two-dimensional images at different positions of the ultrasonic transducer array and/or (ii) reconstructing a volumetric three-dimensional image for a depth scanned by the ultrasonic transducer array. Image reconstruction includes, at least in part, implementing an inverse reconstruction algorithm. Some examples of inverse reconstruction algorithms that can be used include: (i) forward-model-based iterative methods, (ii) time-reversal methods, and (iii) universal back-projection (UBP) method. For example, a 3D back projection algorithm can be used to reconstruct a 3D volumetric image or a 2D back projection algorithm can be used to reconstruct a 2D image. An example of a universal back-projection process can be found in Xu, M. And Wang, L., “Universal back-projection algorithm for photoacoustic computed tomography,” Physical Review E 71, 016706 (2005), which is hereby incorporated by reference in its entirety. Another example of a back-projection process can be found in Anastasio, M. A. et al., “Half-time image reconstruction in thermoacoustic tomography,” IEEE Trans., Med. Imaging 24, pp 199-210 (2005), which is hereby incorporated by reference in its entirety. In another aspect, a dual-speed-of sound (dual-SOS) photoacoustic reconstruction process may be used. Additional examples of inverse reconstruction algorithms are provided in Section I. An example of a single-impulse panoramic photoacoustic computed tomography system that employs a dual-SOS photoacoustic reconstruction process is described in U.S patent application 2019/0307334, titled “SINGLE-IMPULSE PANORAMIC PHOTOACOUSTIC COMPUTED TOMOGRAPHY” and filed on May 29, 2019, which is hereby incorporated by reference in its entirety.

An example of a UBP method that can be implemented in the image reconstruction operation 730 is described with respect to the flowchart shown in FIGS. 9A and 9B. FIG. 9A is a flowchart of operations of a universal back-projection process that can be used to reconstruct either a 2D image or a 3D image, according to an implementation. FIG. 9B is a flowchart of additional operations of the universal back-projection process in FIG. 9A as used for the 3D image, according to an implementation. More details for an example of this process can be found in Anastasio, M. A. et al., “Half-time image reconstruction in thermoacoustic tomography,” IEEE Trans. Med. Imaging 24, 199-210 (2005), which is hereby incorporated by reference in its entirety.

At operation 930, the time delay is calculated for a pixel (voxel) and transducer element location. At operation 940, an acquired photoacoustic signal at the time delay is used to calculate a back-projection term and this is added to the pixel (voxel) value. At operation 950, the method returns to repeat operations 930 and 940 for all combinations of pixel (voxel) and transducer element locations. At operation 960, a 2D (3D) image is formed of all the pixel (voxel) values.

In FIG. 9B, a three-dimensional image output from operation 960 in FIG. 9A is received. At operation 1060, in the elevational direction of each filtered volumetric 3D image, voxels are selected with the largest photoacoustic amplitudes to form a maximum amplitude projection (MAP) 2D image in grayscale. At operation 1070, the depths of the voxels along with their largest photoacoustic amplitude values from operation 1060 are used to form a depth map in grayscale. At operation 1080, the grayscale depth map is transferred to a colorful depth image. At operation 1090, the MAP image is used to modulate the colorful depth image in three channels (RGB), resulting in a color-encoded MAP image. A median filtration with a window size of 3×3 pixels to the depth image. Another median filtration with a window size of 6×6 pixels may be further applied. Different RGB (red, green, blue) color values were assigned to discrete depths. The 2D depth-resolved color-encoded image is multiplied by the MAP image pixel by pixel to represent the maximum amplitudes. In one aspect, to further reduce noise and improve image quality, the parameters are tuned in 2D slices at different depths.

FIG. 8 is a flowchart depicting operations of a PACT method with spatiotemporal antialiasing that applies temporal filtering with a location-dependent antialiasing filter to the temporal signals before applying spatial interpolation. The location-dependent temporal filtering includes at least operations 820, 830, and 840.

At operation 810, one or more control signals are sent to component(s) of the PACT system to perform data acquisition or the photoacoustic data is received or retrieved from a component of the PACT system. For example, photoacoustic data may be retrieved or received from a computer readable media such as, e.g., an onboard buffer.

At operation 820, a distance between the reconstruction location of the pixel (voxel) and the transducer array is determined. At operation 830, it is determined whether the distance determined is places the reconstruction location within the present alias-free zone. In one implementation, the alias-free zone S₂ for a transducer array with a circular geometry may be determined from Eqn. 22. In one implementation, the alias-free zone S₂ for a transducer array with linear geometry may be determined from Eqn. 34. If it is determined that the distance determined places the reconstruction location within the present alias-free zone, the method continues to operation 850 to perform spatial interpolation. If it determined that the distance determined places the reconstruction location outside the present alias-free zone, the method continues to operation 840.

At operation 840, a location-dependent antialiasing filter is applied with an upper cutoff frequency that extends the one-way Nyquist zone (no aliasing in spatial sampling) S₁ to an adjusted one-way Nyquist zone S₁′ that includes the reconstruction location. The antialiasing filter includes one or more low pass filters that can remove components with frequencies higher than the upper cutoff frequency. The antialiasing filter may include a single lowpass filter or a combination of the same type or different types of lowpass filters. Some examples of lowpass filters that be implemented include a Butterworth filter, a Cheyshev filter, and a Bessel filter. In one example, the antialiasing filter includes at least in part a third-order lowpass Butterworth filter and a sinc filter (with the same cutoff frequency). The PACT method determines an upper cutoff frequency based on the distance between the reconstruction location and the transducer array. If the transducer array has a circular geometry, the upper cutoff frequency may be calculated as

$f_c^{\prime }=\frac{Nc}{4{\mathrm{nr}}^{\prime }}$ where N is the number of elements in the transducer array, c is the speed of sound, r′ is the distance between the reconstruction location and the center of the transducer array. If the transducer array has a linear geometry, the upper cutoff frequency may be calculated using

$f_c^{\prime }=\frac{c}{2d}\sqrt{1+{\left(\frac{z}{\max \left\{\begin{array}{c}y-\frac{N-2}{2}d,\\ y+\frac{N-2}{2}d\end{array}\right\}}\right)}^2},$ where d is the pitch of the transducer array, c is the speed of sound, and ƒ_c′ is the upper cutoff frequency. Temporal filtering may be applied to every pixel/voxel location. In one implementation, location-based temporal filtering is applied to one or more groups of pixels/voxels with similar cutoff frequencies. The signals for these groups of pixels/voxels are filtered with different levels of cutoff frequencies before image reconstruction. During reconstruction, location-dependent cutoff frequencies are calculated, and the filtered signals are retrieved directly.

At operation 850, spatial interpolation is applied. The spatial interpolation operation may remove aliasing in image reconstruction to extend the present alias-free two-way Nyquist zone to an adjusted second two-way Nyquist zone that includes the reconstruction location. An example of a spatial interpolation technique that may be used is the frequency domain implementation of Whittaker-Shannon interpolation. Given a time t, signals from all the detection elements formed a vector with length N. A fast Fourier transformation (FFT) is applied to the vector and zeros were padded behind the highest frequency components to double the vector length to 2N in a new vector. The inverse FFT is then applied to the new vector to spatially interpolate the data. In this example, spatial interpolation numerically doubled the number of detection elements.

At optional (denoted by dashed line) operation 860, the method returns to repeat operations 830, 840, and 950 for all combinations of pixel (voxel) locations. In another implementation, location-based temporal filtering may be applied to one or more groups of pixels/voxels with similar cutoff frequencies. The signals for these groups of pixels/voxels are filtered with different levels of cutoff frequencies before image reconstruction. In this implementation, before reconstruction, photoacoustic signals are first filtered with different levels of cutoff frequencies associated with these groups. During reconstruction, location-dependent cutoff frequencies are calculated, and the filtered signals are retrieved directly, which may reduce time in data filtering.

At operation 870, image reconstruction operations are performed including at least on in part implementing an inverse reconstruction algorithm. Image reconstruction may include (i) reconstructing one or more two-dimensional images at different positions of the ultrasonic transducer array and/or (ii) reconstructing a volumetric three-dimensional image for a depth scanned by the ultrasonic transducer array. Some examples of inverse reconstruction algorithms that can be used include: (i) forward-model-based iterative methods, (ii) time-reversal methods, and (iii) universal back-projection (UBP) method. For example, a 3D back projection algorithm can be used to reconstruct a 3D volumetric image or a 2D back projection algorithm can be used to reconstruct a 2D image. An example of a universal back-projection process can be found in Xu, M. And Wang, L., “Universal back-projection algorithm for photoacoustic computed tomography,” Physical Review E 71, 016706 (2005), which is hereby incorporated by reference in its entirety. Another example of a back-projection process can be found in Anastasio, M. A. et al., “Half-time image reconstruction in thermoacoustic tomography,” IEEE Trans., Med. Imaging 24, pp 199-210 (2005), which is hereby incorporated by reference in its entirety. In another aspect, a dual-speed-of sound (dual-SOS) photoacoustic reconstruction process may be used. Additional examples of inverse reconstruction algorithms are provided in Section I. An example of a single-impulse panoramic photoacoustic computed tomography system that employs a dual-SOS photoacoustic reconstruction process is described in U.S patent application 2019/0307334, titled “SINGLE-IMPULSE PANORAMIC PHOTOACOUSTIC COMPUTED TOMOGRAPHY” and filed on May 29, 2019, which is hereby incorporated by reference in its entirety.

IV. Numerical Simulations

Numerical simulations of examples of three-dimensional and two-dimensional detection geometry were performed to evaluate the effects of spatial interpolation and/or location-dependent temporal filtering on mitigating spatial aliasing in image reconstruction.

A. Three-Dimensional Examples

To visualize artifacts caused by spatial aliasing without spatiotemporal antialiasing, a numerical simulation of a three-dimensional example of a hemispherical detection geometry was used. The k-wave toolbox discussed in Treeby B. E. and Cox, B. T., “k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields,” J. Biomed. Opt., vol. 15, no. 2, p. 021314 (2010) was used for the forward problem. The parameters used in the simulation were: the frequency range of the ultrasonic transducer is from 0.1 MHz to 4.5 MHz (2.3-MHz central frequency, 191% bandwidth, ƒ_c=4.5 MHz), the speed of sound (SOS) c=1.5 mm·μ⁻¹, the number of detection elements N=651 (evenly distributed on the hemisphere based on the method in Deserno, M., “How to generate equidistributed points on the surface of a sphere,” Polym. Ed, p. 99 (2004), which is hereby incorporated by reference in its entirety), the radius of the hemisphere R=30 mm, and the simulation grid size is 0.1×0.1×0.1 mm³).

FIG. 10A is an illustration of a numerical simulation of a transducer array with 651 detection elements evenly distributed on a hemisphere and a simple numerical phantom covered by the array, according to an implementation. FIG. 10B is an illustration of a numerical simulation of the transducer array in FIG. 10A with a complex numerical phantom, according to an implementation. FIG. 11 depicts ground truth slices 1102, 1104, 1106, 1108, 1110, and 1112, according to an implementation. FIG. 12 depicts reconstructed images 1202, 1204, 1206, 1208, 1210, and 1212, according to an implementation. Ground truth slice 1102 and reconstructed image 1202 are of a simple phantom at z=0. Ground truth slice 1104 and reconstructed image 1204 are of a simple phantom at z=6.8 mm. Ground truth slice 1106, 1108, 1110, and 1112 and reconstructed image 1206, 1208, 1210, and 1212 are of a complex phantom at z=0, 3.4, 6.8, and 10.2 mm respectively.

In FIG. 10A, the transducer elements are shown as a hemisphere of dots. A simple numerical phantom with nonzero initial pressure from two layers with distances of 0 and 6.8 mm from the xy plane was analyzed first. The phantom is shown as a line of dots in FIG. 10A. The ground-truth images of the two layers are 1102 and 1104, respectively, and the reconstructions are 1202 and 1204, respectively. As can be seen, artifacts appear in both layers and seem stronger in the layer further from the xy plane.

Using the same array, a complex phantom was simulated with nonzero initial pressure from four layers with distances of 0, 3.4, 6.8, and 10.2 mm from the xy plane, shown as a center portion of dots in FIG. 10B. The ground-truth images of the four layers are 1106, 1108, 1110, and 1112 in FIG. 11, respectively, and the reconstructions are 1206, 1208, 1210, and 1212 in FIG. 12 respectively. Strong aliasing artifacts appear in all reconstructed images, and the artifacts tend to be more visible in layers further from the origin. Regions of interest A-G (1.5×1.5 mm²) with increasing distances from the origin were picked from the four layers at locations with zero initial pressure. The standard deviations of the pixel values inside these regions of interest were calculated to quantify the aliasing artifacts.

FIG. 13 is a graph of standard deviations of pixel values in the regions of interest outlined in the boxes shown in reconstructed images 1206, 1208, 1210, and 1212 in FIG. 11, according to an implementation. As can be seen in FIG. 13, artifacts become stronger as the region of interest moves away from the origin.

The three-dimensional (3D) simulation provide a direct observation of spatial aliasing artifacts. In addition, the 3D simulation provides a simplified estimation of the one-way Nyquist zone of the hemisphere geometry. On each of the planes crossing the center of the hemisphere, the one-way Nyquist zone for the 2D case is calculated. All the detection elements for each plane lie on its intersection with the hemisphere, which is a semicircle with length πR (subset of a full circle with length 2πR). Based on the above analysis, one only need to estimate the number of detection elements for this plane. Given the area of the hemisphere 2πR² and the number of elements N, the distance between two neighboring elements is approximately

$\sqrt{\frac{2\pi R^2}{N}}.$ Thus, the number of elements (in a full circle) for this plane is approximately:

$\begin{array}{cc}{N}^{\prime }=\frac{2\pi R}{\sqrt{\frac{2\pi R^2}{N}}}=\sqrt{2\pi N}.& \left(\mathrm{Eqn}35\right)\end{array}$ Using Eqn. 15, the one-way Nyquist zone on this plane is:

$\begin{array}{cc}{S}_1^{\prime }=\left\{r^{\prime }|r\text{'}\le \frac{N^{\prime }{\lambda }_c}{4\pi }\right\}=\left\{r^{\prime }|r^{\prime }\le \frac{c}{4f_c}\sqrt{\frac{2N}{\pi }}\right\}.& \left(\mathrm{Eqn}.36\right)\end{array}$

For a 3D simulation, the plane's normal vector is varied: All the 2D regions S₁′ form a half ball with radius

$\frac{c}{4f_c}\sqrt{\frac{2N}{\pi }}\approx 1.70\mathrm{mm},$ which is much smaller than the radius R=30 mm. This large difference explains the prevalence of artifacts in the reconstructed 3D image.B. Two-Dimensional Examples

In a two-dimensional example, a full-ring transducer array with a radius R=30 mm was used. The frequency range of the full-ring transducer array was from about 0.1 MHz to about 4.5 MHz. The number of detection elements was N=512 and the speed of sound is c=1.5 mm·μs⁻¹. The radius of the one-way Nyquist zone S₁ is thus

$r=\frac{Nc}{4\pi f_c}\approx 13.6\mathrm{mm}.$ In the numerical simulation, a 0.1×0.1 mm² grid size was used. A simple initial pressure distribution was simulated.

Simple Phantom Example

FIG. 14 is an illustration of ground truth of a simple initial pressure p₀ distribution for a simple phantom, according to an implementation. A full-ring transducer array and the circular boundaries of S₁ and S₂ are depicted as circles, respectively. The illustration includes blown up illustrations of boxed regions 1410, 1420, and 1430 in the imaging domain of the full-ring transducer array. The boxed region 1410 lies mainly within zone S₂. The boxed region 1420 lies outside zone S₁ and does not include the simple phantom. The boxed region 1430 lies outside zone S₁ and includes a portion of the simple phantom.

The image of the simple phantom object was reconstructed using a universal back projection algorithm (sometimes referred to as “universal back projection reconstruction”). FIG. 15A is a reconstructed image of the simple phantom object shown in FIG. 14 using a first antialiasing method including universal back projection reconstruction, according to an implementation. Despite the clean ground-truth background in FIG. 14, spatial aliasing-induced artifacts appear in the reconstructed image outside zone S₁, but not inside the box 1410 mainly within zone S₁. The region 1430 is outside zone S₁, but does not show strong artifacts. Due to the specific distribution of the source, the artifacts turn out to be much stronger in the region 1420 than in region 1430.

In a second method, spatial interpolation is used to remove the aliasing in image reconstruction and the image is reconstructed using a universal back projection algorithm. Given a time t, signals from all the detection elements formed a vector with length N. The fast Fourier transformation (FFT) was then applied to the vector, and zeros were padded behind the highest frequency components to double the vector length to 2N. Finally, the inverse FFT was applied to the new vector to interpolate the data. This process is the frequency domain implementation of Whittaker-Shannon interpolation. Spatial interpolation numerically doubled the number of detection elements. FIG. 15B is a reconstructed image of the simple phantom object in FIG. 14 using the second method including UBP reconstruction and spatial interpolation, according to an implementation.

In a third method, to remove the spatial aliasing during spatial sampling, a location-dependent antialiasing filter with one or more lowpass filters was applied to the temporal signals before spatial interpolation. In the simulation, the antialiasing filter included a third-order lowpass Butterworth filter and a sinc filter (with the same cutoff frequency) was implemented for location-dependent temporal filtering. FIG. 15C is a reconstructed image of the simple phantom object in FIG. 14 using a third method including UBP reconstruction, spatial interpolation, and location-dependent temporal filtering, according to an implementation.

To quantify the amplitude of artifacts, regions of interest A-E marked with small boxes in FIG. 15, each region of interest having an area of 1.2×1.2 mm², were selected at locations with zero initial pressure. The standard deviations of the pixel values inside these regions was calculated. FIG. 16 is a graph with a plot of standard deviations of the regions of interest A-E using three methods: (1) a first method with UBP reconstruction alone, (2) a second method with UBP reconstruction and spatial interpolation, and (3) a third method with UBP reconstruction, spatial interpolation (SI), and location-dependent temporal filtering (TF). As shown in FIG. 16, spatial interpolation mitigates the aliasing artifacts significantly, while adding location-dependent temporal filtering before spatial interpolation further diminishes the artifacts. This agrees with the above discussion that both spatial interpolation and location-dependent temporal filtering extend the alias-free region.

To quantify the impact of the antialiasing methods on image resolution, two lines at P and Q, respectively, were selected in each reconstructed image, and their profiles compared in FIGS. 17 and 18 respectively. FIG. 17 is a graph with a plot of profiles of line P using the three methods: (1) a first method with UBP reconstruction alone, (2) a second method with UBP reconstruction and spatial interpolation, and (3) a third method with UBP reconstruction, spatial interpolation (SI), and location-dependent temporal filtering (TF). FIG. 18 is a graph with a plot of profiles of line Q using the three methods. For signals from source points inside zone S₁, there is no spatial aliasing in spatial sampling, and spatial interpolation is accurate at the interpolation points. Moreover, location-dependent temporal filtering does not affect the signals when reconstructing inside zone S₁. Thus, the antialiasing techniques had little impact on the profile of line P (inside zone S₁) in this example, as shown in FIG. 17. For signals from source points outside zone S₁, the spatial interpolation is inaccurate due to spatial aliasing. Thus, directly applying spatial interpolation affects the profile of line Q (outside S₁). Location-dependent temporal filtering smooths the temporal signals before reconstructing outside zone S₁, thus it further smooths the profile of line Q, as shown in FIG. 18.

Complex Phantom Examples

FIG. 19 is an illustration of ground truth of a simple initial pressure p₀ distribution for a first complex phantom, according to an implementation. FIG. 20A is a reconstructed image of the first complex phantom object shown in FIG. 19 using a first antialiasing method including universal back projection reconstruction, according to an implementation. FIG. 20B is a reconstructed image of the first complex phantom object in FIG. 19 using the second method including UBP reconstruction and spatial interpolation, according to an implementation. FIG. 20C is a reconstructed image of the first complex phantom object in FIG. 19 using a third method including UBP reconstruction, spatial interpolation, and location-dependent temporal filtering, according to an implementation. The artifacts in red-boxed region are caused by spatial aliasing in IR, and they are mainly mitigated by SI.

FIG. 21 is a graph with a plot of standard deviations of the regions of interest A-C using three methods: (1) a first method with UBP reconstruction alone, (2) a second method with UBP reconstruction and spatial interpolation, and (3) a third method with UBP reconstruction, spatial interpolation (SI), and location-dependent temporal filtering (TF). FIG. 22 is a graph with a plot of profiles of line P using the three methods. FIG. 23 is a graph with a plot of profiles of line Q using the three methods. The artifacts in the boxed region in are caused by spatial aliasing in SS and IR, and the artifacts are mitigated by TF and SI.

For the first complex phantom case, the object is completely within zone S₁, as shown in FIG. 19. The reconstructed images (reconstructions) of the object in FIG. 19 using the three methods: (1) a first method with UBP reconstruction alone, (2) a second method with UBP reconstruction and spatial interpolation, and (3) a third method with UBP reconstruction, spatial interpolation (SI), and location-dependent temporal filtering (TF) are shown in FIGS. 20A-20C respectively. Regions of interest A-C, each with an area of 1.2×1.2 mm², were chosen at locations with zero initial pressure. The standard deviations were calculated and compared, as shown in FIG. 21. Profiles of lines P and Q are shown in FIG. 22 and FIG. 23, respectively.

FIG. 24 is an illustration of ground truth of a simple initial pressure p₀ distribution for a second complex phantom, according to an implementation. FIG. 25A is a reconstructed image of the second complex phantom object shown in FIG. 24 using a first antialiasing method including universal back projection reconstruction, according to an implementation. FIG. 25B is a reconstructed image of the second complex phantom object in FIG. 24 using the second method including UBP reconstruction and spatial interpolation, according to an implementation. FIG. 25C is a reconstructed image of the second complex phantom object in FIG. 24 using a third method including UBP reconstruction, spatial interpolation, and location-dependent temporal filtering, according to an implementation.

FIG. 26 is a graph with a plot of standard deviations of the regions of interest A-C using three methods: (1) a first method with UBP reconstruction alone, (2) a second method with UBP reconstruction and spatial interpolation, and (3) a third method with UBP reconstruction, spatial interpolation (SI), and location-dependent temporal filtering (TF). FIG. 27 is a graph with a plot of profiles of line P using the three methods. FIG. 28 is a graph with a plot of profiles of line Q using the three methods. The FWHM of the main lobe at was increased from 0.35 mm to 0.48 mm by temporal filtering, while the amplitude was changed from 0.90 to 0.56, respectively.

For the second complex phantom case, the object is outside of zone S₁, and covers most of the area inside the zone S₀, inside the full-ring transducer array, as shown in FIG. 24. The reconstructions of the object in FIG. 24 using the three methods are shown in FIGS. 25A-25C, respectively. The standard deviations of regions of interest A-C are compared in FIG. 26, while the profiles of lines P and Q are shown in FIG. 27 and FIG. 28, respectively. Although spatial interpolation mitigates the aliasing artifacts in both FIG. 20A and FIG. 25A, FIG. 20B and FIG. 25B, visible artifacts remain in the boxed region in FIG. 25C, but not in FIG. 20C.

It can also be seen in FIG. 21 and FIG. 26 that spatial interpolation shows more antialiasing effects on FIG. 20A than on FIG. 25A, while temporal filtering's antialiasing effect on FIG. 25A is more than on FIG. 20A. In fact, the aliasing artifacts in FIG. 20A are solely from the image reconstruction, while those in FIG. 25A are from both the spatial sampling and the image reconstruction. Thus, spatial interpolation works well in antialiasing for FIG. 20A, for which temporal filtering's smoothing effect slightly helps; but not as well for FIG. 25A due to the spatial aliasing in spatial sampling, for which temporal filtering may be necessary. In general, the aliasing artifacts are mitigated by spatial interpolation and further diminished by temporal filtering, as shown in both FIG. 21 and FIG. 26. The antialiasing methods maintain the image resolution well inside zone S₁, as shown in FIG. 22, FIG. 23, and FIG. 27. Due to spatial aliasing in spatial sampling for objects outside zone S₁, the profile of line Q is affected by spatial interpolation. Adding location-dependent temporal filtering further smooths the profile. The full width at half maximum (FWHM) of the main lobe in line Q's profile was increased by temporal filtering while the amplitude was reduced, as shown in FIG. 28. All these observations regarding the complex phantoms agree with the discussions about the simple phantom in FIG. 14.

V. In Vivo Experimental Results

Antialising techniques including spatial interpolation and location-dependent temporal filtering were applied to human breast imaging using a PACT system with spatiotermporal antialiasing. The PACT system employed a 512-element full-ring transducer array (e.g., the 512-element full-ring transducer array with a 110-mm radius, 2.25-MHz central frequency, 95% one-way bandwidth as sold by Imasonic, Inc.). Certain components of the PACT system were similar to those described in Lin, L. Hu, P., Shi, J., Appleton, C. M., Maslov, K., Li, L., Zhang, R., and Wang, L. V., “Single-breath-hold photoacoustic computed tomography of the breast,” Nat. Commun., vol. 9, no. 1, p. 2352, 2018, which is hereby incorporated by reference in its entirety. Based on the point source measurements discussed in Section VII, the cutoff frequency, ƒ_c, is estimated to be about 3.80 Mhz. The acquired signals were filtered by an antialiasing filter including a third-order lowpass Butterworth filter and a sinc filter, both with cutoff frequency 3.80 MHz. Thus, the one-way Nyquist zone S₁ had a radius

$r=\frac{Nc}{4\pi f_c}\approx 16.0\mathrm{mm},$ while the two-way Nyquist zone S₂ has a radius of 8.0 mm. A speed of sound of c=1.49 mm·μs⁻¹ was used in these determinations.

FIG. 29A is a reconstructed image using universal back propagation without either spatial interpolation or location-based temporal filtering. Boundaries of zones S₁ and S₂ are shown as dashed and solid circles, respectively. FIG. 29B is a reconstructed image of the same region in FIG. 29A using universal back propagation with spatial interpolation. FIG. 29C is a reconstructed image of the same region in FIG. 29A using universal back propagation with spatial interpolation and location-based temporal filtering. FIG. 30 is a graph with a plot of profiles of line P using the three methods. FIG. 31 is a graph with a plot of profiles of line Q using the three methods. The FWHM of the main lobe at as increased from 0.44 mm to 0.55 mm by temporal filtering, while the amplitude was changed from 0.64 to 0.38, respectively.

Using universal back propagation, a cross-sectional image of a breast was reconstructed as shown in FIG. 29A. The aliasing artifacts are shown in the peripheral regions in the closeup subset in a boxed region. After applying spatial interpolation of the raw data, the reconstructed image from using universal back propagation is shown in FIG. 29B. Applying location-based temporal filtering and spatial interpolation, the reconstructed image from using universal back propagation is shown in FIG. 29C. As can be seen from a comparison of closeup subsets, the image quality is improved by spatial interpolation, and the aliasing artifacts are further mitigated by location-based temporal filtering. For comparison, the profiles of lines P and Q for the three images are shown in FIGS. 30 and 31 respectively. As shown by the numerical simulation, image resolution outside zone S₁ is compromised by both spatial interpolation and temporal filtering. Temporal filtering smooths the profiles, as shown in FIG. 30 and FIG. 31. Quantitatively, as shown in FIG. 31, temporal filtering increases the FWHM of the main lobe of line Q's profile and reduces the amplitude.

In certain implementations, spatial interpolation maintains the resolution in the one-way Nyquist zone S₁ while mitigating the artifacts caused by aliasing in image reconstruction. Spatial interpolation extends the aliasing-free zone from S₂ to S₁. For objects outside S₁, spatial interpolation may not be useful due to spatial aliasing in spatial sampling. Adding location-dependent temporal filtering does not affect the resolution inside S₁. For objects outside S₁, temporal filtering suppresses high-frequency signals to satisfy the temporal Nyquist sampling requirement. Although reducing the spatial resolution in the affected regions, temporal filtering mitigates aliasing in spatial sampling and makes the spatial interpolation accurate, thus further extends the aliasing-free zone.

VI. Spectrum Analysis

In one aspect, a PACT method includes operations that perform a spectrum analysis of

$\frac{r^{″}-r_n}{c}{h}_e^{″}\left(\frac{r^{″}-r_n}{c}-\frac{r_m^{\prime }-r_n}{c}\right)$ in Eqn. 17. To analyze spatial aliasing in image reconstruction, expression

$\frac{r^{″}-r_n}{c}{h}_e^{″}\left(\frac{r^{″}-r_n}{c}-\frac{r_m^{\prime }-r_n}{c}\right)$ is analyzed, which is a multiplication of a fast-change variable

$h_e^{″}\left(\frac{r^{″}-r_n}{c}-\frac{r_m^{\prime }-r_n}{c}\right)$ and a slow-change variable

$\frac{r^{″}-r_n}{c}.$ Considering that the multiplication is equivalent to a convolution in the frequency domain, the multiplication of

$\frac{r^{″}-r_n}{c}$ causes spectrum change.

To simplify, only the case with

$r_m^{\prime }=\left(r^{\prime },0\right)\mathrm{and}$ $r^{″}=\left(r^{\prime }\cos \left(2\arccos \frac{r^{\prime }}{R}\right),r^{\prime }\sin \left(2\arccos \frac{r^{\prime }}{R}\right)\right)$ is considered. In this case,

$\frac{r^{″}-r_n}{c}-\frac{r_m^{\prime }-r_n}{c}$ can

achieve the maximum sampling step size

$\frac{4\pi r^{\prime }}{N}$ as n varies Eqn. 19, where spatial aliasing is the most severe. Given an element location r(θ)=R (cos θ, sin θ), one defines

$f_1\left(\theta \right)=\frac{r^{″}-r\left(\theta \right)}{c}\mathrm{and}{f}_2\left(\theta \right)=h_e^{″}\left(\frac{r^{″}-r\left(\theta \right)}{c}-\frac{r_m^{\prime }-r\left(\theta \right)}{c}\right).$ Then, the continuous form of

$\frac{r^{″}-r_n}{c}{h}_e^{″}\left(\frac{r^{″}-r_n}{c}-\frac{r_m^{\prime }-r_n}{c}\right)$ can be expressed as ƒ₁(θ) ƒ₂(θ). Applying the Fourier transformation to ƒ₁(ƒ)ƒ₂(θ) and ƒ₂(θ), F(ƒ₁ƒ₂)(ƒ) and F(ƒ₂)(ƒ) are obtained, respectively. F denotes the Fourier transformation operator. The difference between the two normalized spectra is expressed as

$\underset{f}{\max }❘\frac{❘F\left(f_1{f}_2\right)\left(f\right)❘}{\underset{f}{\max }❘F\left(f_1{f}_2\right)\left(f\right)❘}-\frac{❘F\left(f_2\right)\left(f\right)❘}{\underset{f}{\max }❘F\left(f_2\right)\left(f\right)❘}❘,$ which is a function of r′∈[0, R). The difference is calculated between the two normalized spectra for a full-ring array geometry with a radius R=110 mm, and a point source response h_e′ measured through experiments where the speed of sound used was c=1.49 mm·μs⁻¹.

FIG. 32 is a graph with a plot of the difference between the normalized spectra

$\frac{❘F\left(f_1{f}_2\right)\left(f\right)❘}{\underset{f}{\max }❘F\left(f_1{f}_2\right)\left(f\right)❘}\mathrm{and}\frac{❘F\left(f_2\right)\left(f\right)❘}{\underset{f}{\max }❘F\left(f_2\right)\left(f\right)❘}$ for −r′∈[0,0.98R], according to one implementation. As shown in FIG. 32, for −r′∈[0,0.98R], one finds a

$\underset{f}{\max }❘\frac{❘F\left(f_1{f}_2\right)\left(f\right)❘}{\underset{f}{\max }❘F\left(f_1{f}_2\right)\left(f\right)❘}-\frac{❘F\left(f_2\right)\left(f\right)❘}{\underset{f}{\max }❘F\left(f_2\right)\left(f\right)❘}❘<4.5\times 10^{-3}.$ When r′ approaches R, singularity occurs. For a system with a frequency-dependent SNR smaller than

$\frac{1}{4.5\times 10^{-3}}\approx 222,$ which may almost always true for experimental cases discussed herein, the difference between the two normalized spectra is negligible. Therefore, to simplify the problem, the spatial aliasing in

$h_e^{″}\left(\frac{r^{″}-r_n}{c}-\frac{r_m^{\prime }-r_n}{c}\right)$ may represent the spatial aliasing in

$\frac{r^{″}-r_n}{c}{h}_e^{″}\left(\frac{r^{″}-r_n}{c}-\frac{r_m^{\prime }-r_n}{c}\right).$

VII. Estimation of Upper Cutoff Frequency of an Ultrasonic Transducer Array.

In one aspect, a PACT method includes operations that use point source measurements to estimate the upper cutoff frequency ƒ_c of an ultrasonic transducer array. In one operation, an estimation of the point source response is obtained from the point source measurements. In another operation, the response's amplitude and noise level in the frequency domain is quantified and the frequency-dependent signal-to-noise ratio (SNR) is calculated.

The upper cutoff frequency ƒ_c, higher than the central frequency, is selected where the frequency-dependent SNR decreases to one for the first time. The full-ring transducer array is used to acquire the photoacoustic signals generated by a point source located at the array center for J repetitions (e.g., J=100). From each acquisition, N measurements are obtained from N transducer elements (e.g., N=512). Based on Eqn. 5, by ignoring a constant factor, the point source response can be expressed as h_e′(t). The measurement of h_e′(t) by the n-th element in the j-th acquisition is denoted as h_e,j,n′(t). In each acquisition, the mean response is denoted as:

$\begin{array}{cc}{h}_{e,j}^{\prime }\left(t\right)=\frac{1}{N}\sum _{n=1}^N{h}_{e,j,n}^{\prime }\left(t\right),t>0,j=1,2,\dots ,J.& \left(\mathrm{Eqn}.37\right)\end{array}$

The point source response h_e′(t) can be further estimated using

$\begin{array}{cc}\stackrel{\_}{h_e^{\prime }}\left(t\right)=\frac{1}{J}\sum _{j=1}^J{h}_{e,j}^{\prime }\left(t\right),t>0.& \left(\mathrm{Eqn}.38\right)\end{array}$

The normalized value of h_e′(t) is shown in FIG. 33A. By applying the Fourier transformation to h_e,j′(t) and h_e′(t), F(h_e,j′)(ƒ) and F(h_e′(ƒ), can be respectively obtained. The normalized amplitude of F(h_e′)(ƒ) is shown in FIG. 33B.

FIGS. 33A-E are graphs showing the results of operations of a method of estimating an upper cutoff frequency, according to one implementation. FIG. 33A is a graph of a Normalized estimated point source response

$\left(\frac{\stackrel{\_}{h_e^{\prime }}\left(t\right)}{\max \stackrel{\_}{h_e^{\prime }}\left(t\right)}\right),$ according to one implementation. FIG. 33B is a graph of a normalized frequency spectrum of the estimated point source response

$\left(\frac{❘F\left(\stackrel{\_}{h_e^{\prime }}\right)\left(f\right)❘}{\underset{f}{\max }❘F\left(\stackrel{\_}{h_e^{\prime }}\right)\left(f\right)❘}\right),$ according to an implementation. FIG. 33C is a graph of the normalized frequency spectrum of the noise

$\mathrm{STD}\left(\frac{{\sigma }_1\left(f\right)}{\underset{f}{\max }❘F\left(\stackrel{\_}{h_e^{\prime }}\right)\left(f\right)❘}\right),$ according to an implementation. FIG. 33D is a graph of a frequency-dependent signal-to-noise ratio

$\left(SN{R}_1\left(f\right)=\frac{❘F\left(\stackrel{\_}{h_e^{\prime }}\right)\left(f\right)❘}{{\sigma }_1\left(f\right)}\right).$ The upper cutoff frequency is estimated to be ƒ_c=3.80 MHz in this example where the SNR equals to one. FIG. 33E is a graph of frequency-dependent SNR of the derivative of the point source response

$\left(SN{R}_2\left(f\right)=\frac{❘F\left(\stackrel{\_}{h_e^{″}}\right)\left(f\right)❘}{{\sigma }_2\left(f\right)}\right).$ The upper cutoff frequency is also estimated to be ƒ_c=3.80 MHz. FIG. 33F is a graph of Normalized frequency spectra of h_e′ (frequency components with frequencies higher than 3.80 MHz are removed) and h_e″

$\left(\frac{❘F\left(\stackrel{\_}{h_e^{\prime }}\right)\left(f\right)❘}{\underset{f}{\max }❘F\left(h_e^{\prime }\right)\left(f\right)❘}\mathrm{and}\frac{❘F\left(\stackrel{\_}{h_e^{″}}\right)\left(f\right)❘}{\underset{f}{\max }❘F\left(h_e^{″}\right)\left(f\right)❘},\mathrm{respectively}\right).$

Since each pixel value in a reconstructed image is obtained by a weighted summation of the signals from N elements, the noise in h_e,j′(t) (rather than h_e,j,n′(t)) can be used to approximate the noise in the reconstructed image. At a frequency of ƒ, the noise STD in h_e,j′(t) can be estimated as:

$\begin{array}{cc}{\sigma }_1\left(f\right)\approx \sqrt{\frac{1}{J-1}\sum _{j=1}^J{❘F\left(h_{e,j}^{\prime }\right)\left(f\right)-F\left(\stackrel{\_}{h_e^{\prime }}\right)\left(f\right)❘}^2}.& \left(\mathrm{Eqn}.39\right)\end{array}$ Thus, the frequency-dependent SNR can be defined as:

$\begin{array}{cc}{\mathrm{SNR}}_1\left(f\right)=\frac{❘F\left(\stackrel{\_}{h_e^{\prime }}\right)\left(f\right)❘}{{\sigma }_1\left(f\right)}.& \left(\mathrm{Eqn}.40\right)\end{array}$

The normalized value of σ₁(ƒ) and the value of SNR₁(ƒ) are shown in FIGS. 33C and 33D, respectively. The cutoff frequency ƒ_c to be 3.80 MHz (higher than the central frequency 2.25 MHz) was selected where SNR(ƒ) decreases to one for the first time.

To observe the temporal differentiation effect on the spectrum of h_e′, h_e′ in Eqns. 39 and 40 is replaced with h_e″, and obtain σ₂ (ƒ) and SNR₂(ƒ), respectively. As shown in FIG. 33E, the upper cutoff frequency obtained from SNR₂(ƒ) is the same as that from SNR₁(ƒ). In practice, before the UBP reconstruction, the acquired signals may be filtered with an antialiasing filter, e.g., a third-order lowpass Butterworth filter and a sinc filter (both with a cutoff frequency of ƒ_c). Thus, the frequency components with frequencies higher than 3.80 MHz are removed from h_e′. The normalized frequency spectrum of h_e and h_e″ are shown in FIG. 33F. As can be seen, although the spectrum is positively shifted by the temporal differentiation for ƒ<3.80 MHz, the cutoff frequency 3.80 MHz doesn't change.

VIII. Sampling Step Size Approximation

In one aspect, a PACT method includes one or more operations that approximate sampling step size. An approximation of the sampling step size is used in in Eqn. 9, particularly its maximum value

$\frac{2\pi r^{\prime }}{N}.$ The accuracy or this approximation is shown by expressing the sampling step size as a Taylor expansion with the Lagrange remainder. The first-order term is the approximation being used. By analyzing the higher-order terms, the differences are negligible.

FIG. 34 is a schematic diagram of a full-ring transducer array with a radius R, where two adjacent detection element locations r₁ and r₂, and a source point location r′ are marked. These and the following locations are also regarded as vectors from the origin O to them. Vectors r₁ and r₂ form an angle γ, whose bisector intersects with the ring

$r_{\frac{3}{2}}.$ Vector

$r_{\frac{3}{2}}-r^{\prime }$ forms an angle α′ with the tangential dotted line that is perpendicular to vector

$r_{\frac{3}{2}}.$ Vectors

$r_{\frac{3}{2}}$ and r′ form an angle φ, while vectors

$r_{\frac{3}{2}}-r^{\prime }$ and −r′ form an angle β. The angle formed by vectors

$-r_{\frac{3}{2}}$ and

$r^{\prime }-r_{\frac{3}{2}}$ can be expressed as

${\alpha }^{\prime }-\frac{\pi }{2}.$ This graph is used to estimate the sampling step size |∥r′−r₁∥-∥r′−r₂∥|.

FIG. 35 is a graph with a plot of examples of errors calculated using

$r^{\prime }\gamma =\frac{2\pi r^{\prime }}{N}$ to approximate s(r′, θ) for −r′=0.95r₁′≈15.2 mm, according to an implementation. FIG. 35 is a graph with a plot of errors calculated using

$r^{\prime }\gamma =\frac{2\pi r^{\prime }}{N}$ to approximate s(r′, θ) for −r′=6.75r₁′≈107.8 mm, according to an implementation.

For the full-ring transducer array with radius R and centered at the origin O, one considers a source point at r′ and two adjacent detection element locations r₁ and r₂ as shown in FIG. 34. The bisector of the angle formed by vectors r₁ and r₂ intersects with the ring at

$r_{\frac{3}{2}}.$ Vectors r₁ and r₂ form an angle γ; vectors

$r_{\frac{3}{2}}$ and r′ form an angle φ; while vectors

$r_{\frac{3}{2}}-r^{\prime }$ and −r′ form an angle β. Vector

$r_{\frac{3}{2}}-r^{\prime }$ forms an angle α with the tangential clotted line crossing point

$-r_{\frac{3}{2}}$ Thus, the angle formed by vectors

$r_{\frac{3}{2}}.$ and

$r^{\prime }-r_{\frac{3}{2}}$ can be expressed as

${\alpha }^{\prime }-\frac{\pi }{2}.$ Based on the Law of Cosines in triangles Or₁r′ and Or₂r′:

$\begin{array}{cc}r^{\prime }-r_1=\sqrt{R^2+r^{\prime 2}-2{\mathrm{Rr}}^{\prime }\cos \left(\phi +\frac{\gamma }{2}\right)},\mathrm{and}& \left(\mathrm{Eqn}.41\right)\end{array}$ $\begin{array}{cc}r^{\prime }-r_2=\sqrt{R^2+r^{\prime 2}-2{\mathrm{Rr}}^{\prime }\cos \left(\phi -\frac{\gamma }{2}\right)},& \left(\mathrm{Eqn}.42\right)\end{array}$ respectively. To simplify the following expression, functions are defined as:

$\begin{array}{cc}{f}_{k,\phi }\left(\gamma \right)=\frac{r^{\prime }-r_1}{2}=\sqrt{1+k^2-2k\cos \left(\phi +\frac{\gamma }{2}\right)},\mathrm{and}& \left(\mathrm{Eqn}.43\right)\end{array}$ $\begin{array}{cc}g\left(k,\phi \right)={k\left(1+k^2-2k\cos \phi \right)}^{-\frac{1}{2}}\sin \phi ,& \left(\mathrm{Eqn}.44\right)\end{array}$ $\begin{array}{cc}q\left(k,\phi \right)=\frac{1-k^2}{1+k^2-2k\cos \phi }\in \left[\frac{1-k}{1+k},\frac{1+k}{1-k}\right].& \left(\mathrm{Eqn}.45\right)\end{array}$ Here, letting

$k=\frac{r^{\prime }}{R}.$ using the Law or Sines in triangle Or_3/2r′

$\begin{array}{cc}\frac{\sin \phi }{R\sqrt{1+k^2-2k\cos \phi }}=\frac{\sin \beta }{R}=\frac{\sin \left({\alpha }^{\prime }-\frac{\pi }{2}\right)}{r^{\prime }}=\frac{-\cos {\alpha }^{\prime }}{r^{\prime }},& \left(\mathrm{Eqn}.46\right)\end{array}$ One hasg(k,φ)=−cos α′,  (Eqn. 47)and|g(k,φ)|=k sin β|≤k,k∈[0,1),φ∈[0,2π).  (Eqn. 48)One can prove that

$\begin{array}{cc}{f}_{k,\phi }\left(-\gamma \right)=\frac{r^{\prime }-r_2}{R},& \left(\mathrm{Eqn}.49\right)\end{array}$ $\begin{array}{cc}{f^{\prime }}_{k,\phi }\left(\gamma \right)=\frac{1}{2}g\left(k,\phi +\frac{\gamma }{2}\right),& \left(\mathrm{Eqn}.50\right)\end{array}$ $\mathrm{and}$ $\begin{array}{cc}{f}_{k,\phi }^{\mathrm{″\prime }}\left(\gamma \right)=\frac{g\left(k,\phi +\frac{\gamma }{2}\right)}{8}\left(\frac{1}{4}+\frac{1}{3}{q}^2\left(k,\phi +\frac{\gamma }{2}\right)\right).& \left(\mathrm{Eqn}.51\right)\end{array}$ Then, based on the Taylor expansion of ƒ_k,φ(γ)−ƒ_k,φ, (−γ), sampling step size can be expressed as:

$\begin{array}{cc}❘r^{\prime }-r_1-r^{\prime }-r_2❘=R❘f_{k,\phi }\left(\gamma \right)-f_{k,\phi }\left(-\gamma \right)❘=R❘2f_{k,\phi }^{\prime }\left(0\right)\gamma +\frac{f_{k,\phi }^{\mathrm{″\prime }}\left(\gamma \right)+f_{k,\phi }^{\mathrm{″\prime }}\left(-{\gamma }^{\prime }\right)}{6}{\gamma }^3❘=\mathrm{sgn}\left(\sin \phi \right)g\left(k,\phi \right)\mathrm{R\gamma }-\frac{\mathrm{sgn}\left(\sin \phi \right)g\left(k,\phi +\frac{{\gamma }^{\prime }}{2}\right)}{48}\left(\frac{1}{4}+\frac{3}{4}{q}^2\left(k,\phi +\frac{{\gamma }^{\prime }}{2}\right)\right)R{\gamma }^3-\frac{\mathrm{sgn}\left(\sin \phi \right)g\left(k,\phi -\frac{{\gamma }^{\prime }}{2}\right)}{48}\left(\frac{1}{4}+\frac{3}{4}{q}^2\left(k,\phi -\frac{{\gamma }^{\prime }}{2}\right)\right)R{\gamma }^3,\mathrm{for}\mathrm{some}{\gamma }^{\prime }\in \left(0,\gamma \right).& \left(\mathrm{Eqn}.52\right)\end{array}$ The following sign function is used:

$\begin{array}{cc}\mathrm{sgn}\left(x\right)=\{\begin{array}{c}-1,x<0,\\ 0,x=0,\\ 1,x>0.\end{array}& \left(\mathrm{Eqn}.53\right)\end{array}$

In practice, it is the maximum sampling step size that affects spatial aliasing. To have a finer estimation of the maximum sampling step size, first the upper bound of ∥r′−r₁∥-∥r′−r₂∥ is estimated. For

$\phi \in \left[\frac{\gamma }{2},\pi -\frac{\gamma }{2}\right]\bigcup \left[\pi +\frac{\gamma }{2},2\pi -\frac{\gamma }{2}\right],\mathrm{sgn}\left(\mathrm{sin\phi }\right)g\left(k,\phi \pm \frac{{\gamma }^{\prime }}{2}\right)\ge 0.$ Thus, the high order terms in Equation (44) are nonpositive, and one has:

$\begin{array}{cc}||r^{\prime }-r_1||-||r^{\prime }-r_2||\le \left|g\left(k,\phi \right)\right|R\gamma =\left|{\cos \alpha }^{\prime }\right|R\gamma ,k\in \left[0,1\right),\phi \in \left[\frac{\gamma }{2},\pi -\frac{\gamma }{2}\right]\bigcup \left[\pi +\frac{\gamma }{2},2\pi -\frac{\gamma }{2}\right],\gamma =\frac{2\pi }{N}\le \frac{\pi }{4}.& \left(\mathrm{Eqn}.54\right)\end{array}$ N≥8 is assumed. For

$\phi \in \left[-\frac{\gamma }{2},\frac{\gamma }{2}\right],\phi +\frac{\gamma }{2}\in \left[0,\gamma \right]\mathrm{and}\phi -\frac{Y}{2}\in \left[-\gamma ,0\right],$ which means that both ƒ_k,φ(γ) and ƒ_k,φ(−γ) belong to [√{square root over (1+k²−2k)}, √{square root over (1+k²−2k cos γ)}]. Thus:

$\begin{array}{cc}||r^{\prime }-r_1||-\hspace{1em}\hspace{1em}\hspace{1em}||r^{\prime }-r_2||\le \hspace{1em}R\left|\sqrt{1+k^2-2k}-\sqrt{1+k^2-2k\mathrm{cos\gamma }}\right|=R\left|f_{k,\frac{\gamma }{2}}\left(-\gamma \right)-f_{k,\frac{\gamma }{2}}\left(\gamma \right)\right|\le \left|g\left(k,\frac{\gamma }{2}\right)\right|R\gamma \le r^{\prime }\gamma .& \left(\mathrm{Eqn}.55\right)\end{array}$ Here, using Inequality in Eqn. 54 with

$\phi =\frac{\gamma }{2}.$ Similarly, for

$\phi \in \left[\pi -\frac{\gamma }{2},\pi +\frac{\gamma }{2}\right],$ gives us:

$\begin{array}{cc}||r^{\prime }-r_1||-||r^{\prime }-r_2||\le \left|g\left(k,\pi -\frac{\gamma }{2}\right)\right|R\gamma \le r^{\prime }\gamma .& \left(\mathrm{Eqn}.56\right)\end{array}$ Combining Inequalities in Eqns. 54-56, the upper bound of the sampling step size is:

$\begin{array}{cc}||r^{\prime }-r_1||-||r^{\prime }-r_2||\le r^{\prime }\gamma =\frac{2{\pi r}^{\prime }}{N},N\ge 2,r^{\prime }\in \left[0,R\right),\phi \in \left[0,2\pi \right).& \left(\mathrm{Eqn}.57\right)\end{array}$

Next, the lower bound of ∥r′−r₁∥-∥r′r₂∥ is estimated for

$\beta =\frac{\pi }{2}.$ In fact, for each r′>0, there exist N locations of r′ evenly distributed on the circle ∥r′∥=r′ such that

$\beta =\frac{\pi }{2}.$ For each location, φ=arccos k, sgn(sin φ)=1, and g(k, arccos k)=k. In general, there may be singularities in

$q\left(k,\phi +\frac{{\gamma }^{\prime }}{2}\right).$ To avoid the singularities, it is assumed that k≤cos γ, which means that φ=arccos

$k\in \left[\gamma ,\frac{\pi }{2}\right].$ Thus,

$\cos \left(\phi +\frac{{\gamma }^{\prime }}{2}\right)<\mathrm{cos\phi }\mathrm{and}\cos \left(\phi -\frac{{\gamma }^{\prime }}{2}\right)<\cos \left(\phi -\frac{\gamma }{2}\right),$ leading to:

$\begin{array}{cc}q\left(k,\phi +\frac{{\gamma }^{\prime }}{2}\right)<q\left(k,\mathrm{arcos}k\right)=1,& \left(\mathrm{Eqn}.58\right)\\ \mathrm{and}& \\ q\left(k,\phi -\frac{{\gamma }^{\prime }}{2}\right)<q\left(k,\mathrm{arcos}k-\frac{\gamma }{2}\right)=\frac{1-k^2}{1+k^2-2k^2\cos \frac{\gamma }{2}-2k\sqrt{1-k^2}\sin \frac{\gamma }{2}},& \left(\mathrm{Eqn}.59\right)\end{array}$ respectively. One can validate that

$\frac{\partial q\left(k,\arccos k-\frac{\gamma }{2}\right)}{\partial k}\ge 0$ is equivalent to

$k^2\le \frac{1+\cos \frac{\gamma }{2}}{5-3\cos \frac{\gamma }{2}}.$ Thus

$q\left(k,\arccos k-\frac{\gamma }{2}\right),$ as a function of k, is monotonically increasing on

$\left[0,\sqrt{\frac{1+\cos \frac{\gamma }{2}}{5-3\cos \frac{\gamma }{2}}}\right).$ Note that

${\cos }^2\frac{\gamma }{2}<\frac{1+\cos \frac{\gamma }{2}}{5-3\cos \frac{\gamma }{2}}$ is equivalent to

${\left(\cos \frac{\gamma }{2}-1\right)}^2\left(3\cos \frac{\gamma }{2}+1\right)\ge 1,$ which is valid for any

$\gamma =\frac{2\pi }{N}\le \frac{\pi }{4}.$ Thus, for any

$k\le \cos \gamma ,k^2\le {\cos }^2\gamma <{\cos }^2\frac{\gamma }{2}\le \frac{1+\cos \frac{\gamma }{2}}{5-3\cos \frac{\gamma }{2}}.$ Based on the monotonicity of

$q\left(k,\arccos k-\frac{\gamma }{2}\right)$ as a function or k:

$\begin{array}{cc}q\left(k,\arccos k-\frac{\gamma }{2}\right)\le q\left(\cos \gamma ,\arccos \left(\cos \gamma \right)-\frac{\gamma }{2}\right)=\frac{1-{\cos }^2\gamma }{1+{\cos }^2\gamma -2\cos \mathrm{\gamma cos}\frac{\gamma }{2}}.& \left(\mathrm{Eqn}.60\right)\end{array}$ Further, it can be proven that

$q\left(\cos \gamma ,\frac{\gamma }{2}\right)<4$ is equivalent to

${\left(\cos \frac{\gamma }{2}-1\right)}^2\left(5{\cos }^2\frac{\gamma }{2}+6\cos \frac{\gamma }{2}+2\right)>0,$ which is valid for any

$\gamma =\frac{2\pi }{N}\le \frac{\pi }{4}.$ In summary, for any k≤cos γ: and

$\begin{array}{cc}0<q\left(k,\arccos k+\frac{{\gamma }^{\prime }}{2}\right)<1,\mathrm{and}& \left(\mathrm{Eqn}.61\right)\\ 0<q\left(k,\arccos k-\frac{{\gamma }^{\prime }}{2}\right)<4.& \left(\mathrm{Eqn}.62\right)\end{array}$ Combining Inequalities in Eqn. 61 and Eqn. 62 with Eqn. 52, the following it obtained:

$\begin{array}{cc}r^{\prime }-r_1-r^{\prime }-r_2\ge r^{\prime }\gamma -\frac{53r^{\prime }}{192}{\gamma }^3=\frac{2\pi r^{\prime }}{N}\frac{53r^{\prime }}{192}{\left(\frac{2\pi }{N}\right)}^3,N\ge 8,r^{\prime }\le R\cos \gamma ,\phi =\arccos \frac{r^{\prime }}{R}.& \left(\mathrm{Eqn}.63\right)\end{array}$ For any source point=(r′ cos θ, r′ sin θ), the ma sampling step size can be expressed as:

$\begin{array}{cc}s\left(r^{\prime },\theta \right)=\underset{\phi }{\max }r^{\prime }-r_1-r^{\prime }-r_2.& \left(\mathrm{Eqn}.64\right)\end{array}$ From Inequality in Eqn. 57:

$\begin{array}{cc}s\left(r^{\prime },\theta \right)\le \frac{2\pi r^{\prime }}{N},N\ge 8,r^{\prime }\in \left[0,R\right).& \left(\mathrm{Eqn}.65\right)\end{array}$

According to Inequality in Eqn. 63, there exist at least N values of θ evenly distributed in [0,2π) such that:

$\begin{array}{cc}s\left(r^{\prime },\theta \right)\ge \frac{2\pi r^{\prime }}{N}-\frac{53r^{\prime }}{192}{\left(\frac{2\pi }{N}\right)}^3,N\ge 8,r^{\prime }\le R\cos \gamma .& \left(\mathrm{Eqn}.66\right)\end{array}$

As can be seen in Eqn. 65, using

$r^{\prime }\gamma =\frac{2\pi r^{\prime }}{N}$ to approximate s(r′, θ) is sufficient for the antialiasing analysis; while

$\frac{53r^{\prime }}{192}{\left(\frac{2\pi }{N}\right)}^3{\left(\frac{{\lambda }_c}{2}\right)}^{-1}$ in Eqn. 66 can be used to estimate the necessity to further reduce the upper bound of s(r′, θ) for the N specific values of θ. For general values of θ, numerical simulation can be used to estimate.

Next, these estimations are evaluated in numerical simulations and in vivo experiments. In both cases, N=512. The radius constraint, from inequality (58),

$\frac{r^{\prime }}{R}\le \cos \gamma =\cos \frac{2\pi }{N}\approx 1-7.5\times {10}^{-5}$ is substantially close to the intrinsic constraint

$\left(\frac{r^{\prime }}{R}<1\right).$ For the numerical simulations, c=1.5 mm·μs⁻¹, R=30 mm, and ƒ_c=4.5 MHz, thus

$\frac{53r^{\prime }}{192}{\left(\frac{2\pi }{N}\right)}^3{\left(\frac{{\lambda }_c}{2}\right)}^{-1}<\frac{53R}{192}{\left(\frac{2\pi }{N}\right)}^3\left(\frac{2f_c}{c}\right)<9.2\times {10}^{-5}.$ For in vivo experiments, c=1.49 mm·μs⁻¹, R=110 mm, and ƒ_c=3.80 MHz, then

$\frac{53r^{\prime }}{192}{\left(\frac{2\pi }{N}\right)}^3{\left(\frac{{\lambda }_c}{2}\right)}^{-1}<2.9\times {10}^{-4}.$ In both cases, for the N specific values of θ, using

$\frac{2\pi r^{\prime }}{N}$ to approximate s(r′, θ) is substantially accurate. For general values of θ, numerical simulation are used to observe the approximation error for source locations on two circles (in vivo experiment cases): one with a radius of

$0.95{r^{\prime }}_1=0.95\frac{N{\lambda }_c}{4\pi }\approx 15.2\mathrm{mm}$ (r₁′ is the radius of S₁), the other with a radius of 6.75r₁′≈107.8 mm (98% of R).

FIG. 35 is a graph with plot of

$e\left(r^{\prime },\theta \right)=\left(r^{\prime }\gamma -s\left(r^{\prime },\theta \right)\right){\left(\frac{{\lambda }_c}{2}\right)}^{-1}$ against

$\theta \in \left[0,\frac{\pi }{60}\right]$ for r′=0.95r₁′, according to an implementation. FIG. 36 is a graph with plot of e(r′, θ)=

$\left(r^{\prime }\gamma \right)-s\left(r^{\prime },\theta \right)){\left(\frac{{\lambda }_c}{2}\right)}^{-1}$ against

$\theta \in \left[0,\frac{\pi }{60}\right]\mathrm{for}{r}^{\prime }=6.75r_1^{\prime },$ according to an implementation.

Each bar marks a choice of θ assumed in inequality in Eqn. 66. One can see from FIGS. 35 and 36 that, for general values of θ, to further reduce the upper bound of s(r′, θ) is unnecessary.

Modifications, additions, or omissions may be made to any of the above-described embodiments without departing from the scope of the disclosure. Any of the embodiments described above may include more, fewer, or other features without departing from the scope of the disclosure. Additionally, the steps of described features may be performed in any suitable order without departing from the scope of the disclosure. Also, one or more features from any embodiment may be combined with one or more features of any other embodiment without departing from the scope of the disclosure. The components of any embodiment may be integrated or separated according to particular needs without departing from the scope of the disclosure.

It should be understood that certain aspects described above can be implemented in the form of logic using computer software in a modular or integrated manner. Based on the disclosure and teachings provided herein, a person of ordinary skill in the art will know and appreciate other ways and/or methods to implement the present invention using hardware and a combination of hardware and software.

Any of the software components or functions described in this application, may be implemented as software code using any suitable computer language and/or computational software such as, for example, Java, C, C#, C++ or Python, LabVIEW, Mathematica, or other suitable language/computational software, including low level code, including code written for field programmable gate arrays, for example in VHDL. The code may include software libraries for functions like data acquisition and control, motion control, image acquisition and display, etc. Some or all of the code may also run on a personal computer, single board computer, embedded controller, microcontroller, digital signal processor, field programmable gate array and/or any combination thereof or any similar computation device and/or logic device(s). The software code may be stored as a series of instructions, or commands on a CRM such as a random access memory (RAM), a read only memory (ROM), a magnetic media such as a hard-drive or a floppy disk, or an optical media such as a CD-ROM, or solid stage storage such as a solid state hard drive or removable flash memory device or any suitable storage device. Any such CRM may reside on or within a single computational apparatus, and may be present on or within different computational apparatuses within a system or network. Although the foregoing disclosed embodiments have been described in some detail to facilitate understanding, the described embodiments are to be considered illustrative and not limiting. It will be apparent to one of ordinary skill in the art that certain changes and modifications can be practiced within the scope of the appended claims.

The terms “comprise,” “have” and “include” are open-ended linking verbs. Any forms or tenses of one or more of these verbs, such as “comprises,” “comprising,” “has,” “having,” “includes” and “including,” are also open-ended. For example, any method that “comprises,” “has” or “includes” one or more steps is not limited to possessing only those one or more steps and can also cover other unlisted steps. Similarly, any composition or device that “comprises,” “has” or “includes” one or more features is not limited to possessing only those one or more features and can cover other unlisted features.

All methods described herein can be performed in any suitable order unless otherwise indicated herein or otherwise clearly contradicted by context. The use of any and all examples, or exemplary language (e.g. “such as”) provided with respect to certain embodiments herein is intended merely to better illuminate the present disclosure and does not pose a limitation on the scope of the present disclosure otherwise claimed. No language in the specification should be construed as indicating any non-claimed element essential to the practice of the present disclosure.

Groupings of alternative elements or embodiments of the present disclosure disclosed herein are not to be construed as limitations. Each group member can be referred to and claimed individually or in any combination with other members of the group or other elements found herein. One or more members of a group can be included in, or deleted from, a group for reasons of convenience or patentability. When any such inclusion or deletion occurs, the specification is herein deemed to contain the group as modified thus fulfilling the written description of all Markush groups used in the appended claims.

Claims

1. A photoacoustic computed tomography method, comprising: acquiring photoacoustic data recorded by one or more data acquisition devices of a photoacoustic computed tomography system;applying location-based temporal filtering to the photoacoustic data acquired, wherein applying the location-based temporal filtering comprises determining an upper cutoff frequency based at least in part on geometry of an ultrasonic transducer array of the photoacoustic computed tomography system andapplying one or more lowpass filters with the upper cutoff frequency; andreconstructing one or more photoacoustic images from the filtered photoacoustic data.

2. The photoacoustic computed tomography method of claim 1, wherein the upper cutoff frequency is determined based on a distance from a reconstruction location to a center of the ultrasonic transducer array.

3. The photoacoustic computed tomography method of claim 1, further comprising applying spatial interpolation after the application of the location-based temporal filtering.

4. The photoacoustic computed tomography method of claim 3, wherein universal back projection is used to reconstruct the one or more photoacoustic images.

5. A photoacoustic computed tomography method, comprising: acquiring photoacoustic data;determining whether a reconstruction location is (a) at or inside, or (b) outside an alias-free region; andif it is determined that the reconstruction location is outside the alias-free region,(i) applying location-based temporal filtering to the photoacoustic data acquired, wherein the application of the location-based temporal filtering comprises applying one or more lowpass filters with an upper cutoff frequency based at least in part on geometry of an ultrasonic transducer array; and(ii) reconstructing one or more photoacoustic images from the filtered photoacoustic data.

6. The photoacoustic computed tomography method of claim 5, wherein the determination of whether the reconstruction location is (a) at or inside, or (b) outside the alias-free region is based on Nyquist sampling criterion.

7. The photoacoustic computed tomography method of claim 5, wherein the method further comprises applying spatial interpolation after the application of the location-based temporal filtering if it is determined that the reconstruction location is outside the alias-free region.

8. The photoacoustic computed tomography method of claim 5, wherein the one or more photoacoustic images are reconstructed from the filtered photoacoustic data using universal back propagation.

9. A non-transitory computer readable media for generating one or more photoacoustic images of a specimen being imaged, the non-transitory computer readable media, when read by one or more processors, is configured to perform operations comprising: acquiring photoacoustic data recorded by one or more data acquisition devices of a photoacoustic computed tomography system;applying location-based temporal filtering to the photoacoustic data acquired, wherein applying the location-based temporal filtering comprisesdetermining an upper cutoff frequency based at least in part on geometry of an ultrasonic transducer array of the photoacoustic computed tomography system andapplying one or more lowpass filters with the upper cutoff frequency; andreconstructing one or more photoacoustic images from the filtered photoacoustic data.

10. The non-transitory computer readable media of claim 9, wherein the upper cutoff frequency is determined based on a distance from a reconstruction location to a center of the ultrasonic transducer array.

11. The non-transitory computer readable media of claim 9, wherein the operations further comprise applying spatial interpolation after the application of the location-based temporal filtering.

12. The non-transitory computer readable media of claim 9, wherein universal back propagation is used to reconstruct the one or more photoacoustic images.

13. A non-transitory computer readable media for generating one or more photoacoustic images of a specimen being imaged, the non-transitory computer readable media, when read by one or more processors, is configured to perform operations comprising: acquiring photoacoustic data;determining whether a reconstruction location is (a) at or inside, or (b) outside an alias-free region; andif it is determined that the reconstruction location is outside the alias-free region,(i) applying location-based temporal filtering to the photoacoustic data acquired, wherein the application of the location-based temporal filtering comprises applying one or more lowpass filters with an upper cutoff frequency based at least in part on geometry of an ultrasonic transducer array; and(ii) reconstructing one or more photoacoustic images from the filtered photoacoustic data.

14. The non-transitory computer readable media of claim 13, wherein the determination of whether the reconstruction location is (a) at or inside, or (b) outside the alias-free region is based on Nyquist sampling criterion.

15. The non-transitory computer readable media of claim 13, wherein the operations further comprise applying spatial interpolation after the application of the location-based temporal filtering if it is determined that the reconstruction location is outside the alias-free region.

16. The non-transitory computer readable media of claim 13, wherein the one or more photoacoustic images are reconstructed using universal back propagation.