A rapid three-dimensional photoacoustic imaging system and a reconstruction method thereof

Abstract

This invention discloses a rapid three-dimensional photoacoustic imaging system and its reconstruction method. The system includes: a pulse excitation unit for emitting pulsed laser light to excite biological tissue to generate broadband photoacoustic signals; a three-dimensional calibration motion unit for scanning the entire imaging field of view with a standard point sound source during the system initialization calibration phase to obtain a spatial non-uniform point spread function dictionary characterizing the system response; a sensing and encoding module including an arc-shaped ultrasonic transducer and a non-uniform perforated encoding mask closely attached to the front end of the arc-shaped ultrasonic transducer, the non-uniform perforated encoding mask being used to spatially pre-modulate the incident sound wave, causing the signals from multiple sound sources at different locations to be aliased and output on a small number of receiving channels; and a data acquisition and processing unit for receiving the aliased signals and performing rapid image reconstruction. This invention achieves three-dimensional photoacoustic imaging with extremely simplified hardware, automated calibration, and rapid reconstruction, significantly simplifying the system architecture and reducing manufacturing costs.

Description

Technical Field

[0001] This invention belongs to the field of biomedical imaging technology, specifically relating to a rapid three-dimensional photoacoustic imaging system and its reconstruction method. Background Technology

[0002] Traditional 3D photoacoustic computed tomography (3D-PACT) systems typically require arrays of hundreds or even thousands of channels and extremely complex parallel acquisition circuits to obtain high-quality images. This results in extremely high system costs, complex control logic, and massive data redundancy. Therefore, developing a 3D photoacoustic imaging system with a simple hardware architecture, streamlined control circuitry, fast reconstruction speed, and low cost has significant clinical and research value. Summary of the Invention

[0003] The purpose of this invention is to address the technical problems of existing three-dimensional photoacoustic imaging systems, which rely on high-channel-count transducer arrays and complex parallel acquisition circuits, resulting in high costs, complex control logic, and slow reconstruction speed. The invention provides a fast three-dimensional photoacoustic imaging system and its reconstruction method.

[0004] The objective of this invention is achieved through the following technical solution: A first aspect of this invention provides a fast three-dimensional photoacoustic imaging system, comprising: The pulse excitation unit includes a nanosecond-level pulsed laser for emitting pulsed lasers to excite biological tissues to generate broadband photoacoustic signals. The three-dimensional calibration motion unit includes a three-axis precision displacement stage, which is used to carry a standard point sound source to scan the entire imaging field of view during the system initialization calibration phase in order to obtain a dictionary of spatial non-uniform point diffusion functions characterizing the system response. The sensing and coding module includes an arc-shaped ultrasonic transducer and a non-uniform perforated coding mask that is closely attached to the front end of the arc-shaped ultrasonic transducer. The non-uniform perforated coding mask is used to spatially pre-modulate the incident sound wave so that the signals from multiple sound sources are aliased and output on a small number of receiving channels. The data acquisition and processing unit, including a multi-channel synchronous acquisition board and a host computer, is used to receive aliased signals and perform fast image reconstruction.

[0005] Furthermore, the design method for the non-uniform punch code mask specifically includes: Based on the element spacing of the arc-shaped ultrasonic transducer and the radius of curvature of the probe, calculate its projection resolution unit on the mask sphere. Set a minimum center-to-center distance for the holes so that the acoustic responses produced by different hole positions can be effectively distinguished. A non-uniform perforation density distribution strategy along the axial direction is introduced to make the perforation probability of the two ends of the mask higher than that of the center region, in order to compensate for the decrease in sensitivity of the edge field of view of the arc-shaped ultrasonic transducer caused by the attenuation of directivity. Perform boundary safety checks on all candidate aperture locations and eliminate incomplete apertures that are close to the mask edge.

[0006] Furthermore, when the data acquisition and processing unit performs fast image reconstruction, it employs a fast deconvolution reconstruction algorithm based on a dictionary of measured point spread functions. This algorithm includes: The dictionary of diffusion functions of measured points obtained during the calibration phase is used as the forward model of the system; During the reconstruction process, spatially related gain compensation weights are embedded in the deconvolution operator to dynamically enhance weak edge signals. An accelerated version of the Richardson-Lucy iterative framework is used to perform matrix decoupling operations on the acquired coded signals and the point spread function dictionary to reconstruct the three-dimensional initial pressure distribution of the object.

[0007] Furthermore, the data acquisition and processing unit is also used for: After reconstructing the initial three-dimensional pressure distribution, a three-dimensional adaptive threshold segmentation process is performed. The noise suppression threshold is dynamically set by calculating the standard deviation of global voxels to improve the contrast between biological structures and background tissues. By using coordinate space resampling and interpolation methods, non-uniformly distributed computational voxels are mapped to the standard Cartesian coordinate system; The maximum intensity projection algorithm or volume rendering technique is applied to the processed three-dimensional initial pressure distribution to generate the final three-dimensional photoacoustic image.

[0008] A second aspect of this invention provides a rapid three-dimensional photoacoustic imaging reconstruction method, applied to the aforementioned system, comprising the following steps: (1) System setup and initialization: Configure the pulse excitation unit, sensing and encoding module, data acquisition and processing unit and three-dimensional calibration motion unit to build and initialize the system; (2) Point spread function dictionary calibration: Under the condition of no biological sample, the three-dimensional calibration motion unit carries a standard point sound source to scan the entire imaging field of view and obtain the spatial non-uniform point spread function dictionary characterizing the system response; (3) Encoded signal acquisition: The sample to be tested is placed in the imaging field of view, and the photoacoustic signal is spatially premodulated by the non-uniform perforation encoding mask of the sensing and encoding module through single pulse laser excitation, and the aliasing signal is acquired by the data acquisition and processing unit. (4) Fast 3D reconstruction: Based on the point spread function dictionary and the acquired aliasing signal, the accelerated Richardson-Lucy iterative algorithm with embedded spatial position correlation gain compensation weight is used to decouple and reconstruct the 3D initial pressure distribution of the sample to be tested.

[0009] Furthermore, the point spread function dictionary calibration specifically includes: The standard point sound source is fixed at the end of the three-axis precision displacement stage, and the Z-axis of the arc-shaped ultrasonic transducer is determined according to the physical installation orientation of the three-axis precision displacement stage relative to the arc-shaped ultrasonic transducer. Without installing the encoding mask, the mapping relationship between the motion space of the triaxial precision displacement stage and the acoustic coordinate system of the arc-shaped ultrasonic transducer is calibrated using a transducer array self-localization method, so as to determine the projection of the triaxial precision displacement stage in the acoustic coordinate system. Install a non-uniform perforation coding mask and perform geometric consistency verification on the pose deviation during the installation process; The three-axis precision displacement stage is driven to perform automatic point-by-point scanning according to the preset spatial step size to traverse the entire imaging field of view and record the multi-channel encoded signal after mask modulation corresponding to each spatial position to form a time-domain feature vector. All time-domain feature vectors are indexed by spatial position to construct a dictionary of spread functions for the measured points.

[0010] Furthermore, the method based on transducer array self-localization specifically includes: By utilizing the redundancy of the entire array data, the flight time vector of the sound wave arriving at each array element is extracted from the original multi-channel radio frequency signal. An overdetermined nonlinear equation set is constructed using the spatial coordinates of the point sound source, the inherent electrical delay of the system, and the equivalent sound velocity at the current water temperature as the parameters to be determined. The Levenberg-Marquardt optimization algorithm was used to iteratively solve the overdetermined nonlinear equations to obtain the projected coordinates of the triaxial precision displacement stage in the acoustic coordinate system.

[0011] Furthermore, the iterative update formula for the accelerated Richardson-Lucy iterative algorithm with embedded spatial location-related gain compensation weights is as follows:

[0012] In the formula, This represents the three-dimensional pressure distribution obtained after the k-th iteration. Let f denote the weighted point spread function matrix, the superscript T denotes the transpose of the matrix, and f denotes the preprocessed aliased signal. To prevent constants from being divided by zero, For a full 1 tensor, This indicates element-wise multiplication.

[0013] Furthermore, image post-processing is included after the rapid 3D reconstruction: The reconstructed 3D data volume is subjected to 3D adaptive threshold segmentation to suppress background noise. Data is mapped to the standard Cartesian coordinate system using coordinate space resampling and interpolation methods. The final 3D photoacoustic image is generated using the maximum intensity projection algorithm or volume rendering technology.

[0014] Compared with existing technologies, the beneficial effects of this invention are as follows: This invention utilizes simplified spatial coding hardware to realize a system, mask design process, and corresponding image reconstruction for rapid three-dimensional photoacoustic computed tomography (PACT); by introducing a spatial coding mechanism at the hardware level, this invention replaces the traditional electronic channel stacking with a physical structure, significantly simplifying the system architecture and reducing manufacturing costs; through a three-in-one technical approach of "physical mask coding + measured PSF modeling + lightweight deconvolution," this invention achieves three-dimensional photoacoustic imaging with extremely simplified hardware, automated calibration, and rapid reconstruction, demonstrating outstanding clinical translational potential and engineering application value. Attached Figure Description

[0015] Figure 1 This is a schematic diagram of the overall structure of the fast three-dimensional photoacoustic imaging system of the present invention; Figure 2 This is a design drawing of the non-uniform perforated coding mask of the present invention; Figure 3 This is a flowchart of the rapid three-dimensional photoacoustic imaging reconstruction method of the present invention. Detailed Implementation

[0016] Exemplary embodiments will now be described in detail, examples of which are illustrated in the accompanying drawings. When the following description relates to the drawings, unless otherwise indicated, the same numerals in different drawings denote the same or similar elements. The embodiments described in the following exemplary embodiments do not represent all embodiments consistent with the present invention. Rather, they are merely examples of apparatuses and methods consistent with some aspects of the invention as detailed in the appended claims.

[0017] The terminology used herein is for the purpose of describing particular embodiments only and is not intended to be limiting of the invention. The singular forms “a,” “the,” and “the” used in this invention and the appended claims are also intended to include the plural forms unless the context clearly indicates otherwise. It should also be understood that the term “and / or” as used herein refers to and includes any or all possible combinations of one or more of the associated listed items.

[0018] It should be understood that although the terms first, second, third, etc., may be used in this invention to describe various information, this information should not be limited to these terms. These terms are only used to distinguish information of the same type from one another. For example, first information may also be referred to as second information without departing from the scope of this invention, and similarly, second information may also be referred to as first information. Depending on the context, the word "if" as used herein may be interpreted as "when," "when," or "in response to a determination."

[0019] The present invention will now be described in detail with reference to the accompanying drawings. Unless otherwise specified, the features of the following embodiments and implementations can be combined with each other.

[0020] Example 1: Fast 3D Photoacoustic Imaging System See Figure 1 The rapid three-dimensional photoacoustic imaging system of the present invention includes a pulse excitation unit, a three-dimensional calibration motion unit, a sensing and encoding module, and a data acquisition and processing unit. The pulse excitation unit includes a nanosecond-level pulsed laser for emitting pulsed laser light to excite biological tissue to generate broadband photoacoustic signals. The three-dimensional calibration motion unit includes a triaxial precision displacement stage for scanning the entire imaging field of view with a standard point sound source during the system initialization calibration phase to obtain a spatial non-uniform point spread function (PSF) dictionary characterizing the system response. The sensing and encoding module includes an arc-shaped ultrasonic transducer and a non-uniformly punched coded mask positioned close to the front end of the arc-shaped ultrasonic transducer. This non-uniformly punched coded mask is used to spatially pre-modulate the incident sound wave, causing aliasing of sound source signals from multiple locations on a small number of receiving channels, thereby significantly reducing the requirements for the number of back-end acquisition channels and circuit complexity. The data acquisition and processing unit includes a multi-channel synchronous acquisition board and a host computer for receiving aliased signals and performing rapid image reconstruction.

[0021] It should be understood that the three-axis precision displacement stage is equipped with three precision motors, which achieve precision displacement on the lead screw. The lead screw has a photoelectric gate; the precision motor returns to the photoelectric gate before each displacement to achieve calibration.

[0022] Furthermore, the design method for non-uniform punched coding masks specifically includes: First, based on the element spacing pitch of the arc-shaped ultrasonic transducer and the radius of curvature of the probe... Calculate its projection resolution unit on the mask sphere. :

[0023] In the formula, Let be the radius of curvature of the mask sphere. The "equivalent projection size" representing a single arc-shaped ultrasonic transducer element on the mask sphere is the theoretical lower limit of the mask coding spatial resolution. In this embodiment, Pitch = 0.6mm. , Subsequently, the minimum center-to-center distance of the holes is set to satisfy... (In this embodiment, 0.75mm is used) to effectively distinguish the acoustic responses generated by different aperture positions, thereby reducing the difficulty of signal decoupling at the physical level; wherein This represents the minimum center-to-center distance between any two holes. Its core function is to ensure, at the physical level, that the acoustic responses generated by different holes are distinguishable at the transducer end. The physical diameter of a single hole in the encoding mask is given. Based on this, a non-uniform perforation density distribution strategy along the axial direction is introduced, making the perforation probability higher in the two ends of the mask (Z=±10mm) than in the center region. This compensates for the sensitivity decrease caused by directivity attenuation at the edge of the curved ultrasonic transducer's field of view. Finally, boundary safety checks are performed on all candidate hole positions, eliminating incomplete holes near the mask edge to ensure the mathematical integrity and physical realizability of the encoding operator. Thus, a precise non-uniform perforation encoding mask, strictly aligned with the geometric characteristics of the curved ultrasonic transducer, is obtained, as shown below. Figure 2 As shown.

[0024] Furthermore, when the data acquisition and processing unit performs fast image reconstruction, it employs a fast deconvolution reconstruction algorithm based on a measured point spread function dictionary. This algorithm specifically includes: using the measured point spread function dictionary obtained during the calibration phase as the system's forward model; embedding spatially related gain compensation weights into the deconvolution operator during reconstruction to dynamically enhance weak edge signals; subsequently, using an accelerated Richardson-Lucy iterative framework, performing matrix decoupling operations on the acquired coded signal and the point spread function dictionary, requiring only 5 to 10 iterations to faithfully reconstruct the object's initial 3D pressure distribution. This algorithm abandons computationally intensive neural networks or large-scale iterative optimization methods, has low computational cost, and can achieve near real-time 3D reconstruction on ordinary personal computers, significantly improving imaging efficiency.

[0025] Furthermore, the data acquisition and processing unit is also used to: after reconstructing the three-dimensional initial pressure distribution, perform three-dimensional adaptive threshold segmentation processing, dynamically set the noise suppression threshold by calculating the standard deviation of global voxels to improve the contrast between biological structures and background tissues; map the non-uniformly distributed computational voxels to the standard Cartesian coordinate system through coordinate space resampling and interpolation methods; and perform maximum intensity projection algorithm or volume rendering technology on the processed three-dimensional initial pressure distribution to generate the final three-dimensional photoacoustic image.

[0026] Example 2: A Reconstruction Method for Fast 3D Photoacoustic Imaging The rapid three-dimensional photoacoustic imaging reconstruction method of this invention comprises four stages: system setup and initialization, point spread function dictionary calibration, coded signal acquisition, and rapid three-dimensional reconstruction. The reconstruction process is detailed as follows: Figure 3 As shown below, the operation process of each step is explained in detail using in vivo imaging of blood vessels in the mouse ear as an example.

[0027] (1) System setup and initialization: Configure the pulse excitation unit, sensing and encoding module, data acquisition and processing unit and three-dimensional calibration motion unit to build and initialize the system.

[0028] like Figure 1 As shown, a nanosecond-level pulsed laser is first fixed to an optical platform, and the pulsed laser beam emitted by the nanosecond-level pulsed laser is guided to the top of the imaging tank via a mirror assembly. The nanosecond-level pulsed laser has a wavelength of 532-1100 nm, a repetition rate of 10 Hz, and a pulse width of 8 ns. The mirror assembly is integrated onto the nanosecond-level pulsed laser. Then, an arc-shaped ultrasonic transducer with a center frequency of 7.5 MHz and 256 array elements is installed at the bottom of the tank. The probe curvature radius of the arc-shaped ultrasonic transducer is... It measures 50.7mm and effectively receives solid angles of ±30°. Subsequently, a customized piece of... Figure 2 The non-uniform perforation coding mask shown is tightly attached to the front end of the curved ultrasonic transducer using a 3D-printed bracket; the non-uniform perforation coding mask is made of photosensitive resin, with a thickness of 0.5 mm and a radius of curvature of [missing information]. The array elements of the 256-element arc-shaped ultrasonic transducer are located directly below a non-uniformly perforated coded mask, with the entire structure—the sensing and encoding modules—contained within the water tank. Furthermore, the outputs of all 256 elements of the arc-shaped ultrasonic transducer are connected to a 256-channel synchronous acquisition board. This board has a sampling rate of 40 MS / s and an analog bandwidth of 0-20 MHz. Simultaneously, a three-axis precision stage is positioned above the water tank for subsequent calibration scanning and precise sample positioning; the stage has a travel of 50 mm × 50 mm × 30 mm and a positioning resolution of 1 μm. On the software side, the host computer control program (developed using a hybrid Python and LabVIEW architecture) is launched, loading motion control, data acquisition, and image reconstruction modules. Strict synchronization between the pulsed laser trigger signal and the acquisition is configured to ensure that each laser pulse corresponds to a complete 256-channel signal capture.

[0029] (2) Point spread function dictionary calibration: Under the condition of no biological sample, the three-dimensional calibration motion unit carries a standard point sound source to scan the entire imaging field of view and obtain the spatial non-uniform point spread function dictionary characterizing the system response.

[0030] Further, the point spread function dictionary calibration specifically includes: fixing a standard point sound source at the end of a triaxial precision displacement stage, and determining the Z-axis of the arc-shaped ultrasonic transducer based on the physical installation orientation of the triaxial precision displacement stage relative to the arc-shaped ultrasonic transducer; accurately calibrating the mapping relationship between the motion space of the triaxial precision displacement stage and the acoustic coordinate system of the arc-shaped ultrasonic transducer using a transducer array self-positioning method without installing the coded mask, so as to determine the projection of the triaxial precision displacement stage in the acoustic coordinate system; installing a non-uniformly perforated coded mask, and performing geometric consistency verification on the pose deviation during the installation process; driving the triaxial precision displacement stage to perform point-by-point automatic scanning according to a preset spatial step size to traverse the entire imaging field of view, and recording the multi-channel coded signal after mask modulation corresponding to each spatial position to form a temporal feature vector, and constructing a measured point spread function dictionary by indexing all temporal feature vectors according to spatial position.

[0031] Specifically, an optical fiber with a diameter of approximately 50 μm is fixed to the end of a triaxial precision displacement stage, and carbon particles are attached to the end of the fiber as a standard point sound source. The calibration space volume is set as follows: With a spatial step size of 0.2 mm, a total of 125,000 discrete spatial sampling points were generated. A three-axis precision displacement stage sequentially moved the standard point sound source to each preset position (i.e., a discrete spatial sampling point), and laser irradiation of the carbon fiber generated a localized photoacoustic signal. At this time, the masked modulated sound wave was received by an arc-shaped ultrasonic transducer, and the corresponding aliasing signal was recorded by a 256-channel synchronous acquisition board. ,in This represents the aliased signal at time t in the i-th path. This indicates the number of time sampling points. All acquired aliased signals are indexed according to their corresponding spatial coordinates to construct a five-dimensional tensor form dictionary of measured point spread functions. ,in , and These represent the number of discrete voxels in the imaging region along the x, y, and z directions, respectively. This is the diffusion function dictionary for the measured point. A complete forward model of the entire system, including mask diffraction effect, directional response of arc ultrasonic transducer and electronic link gain, was characterized.

[0032] It is worth noting that the point spread function dictionary calibration process described above only needs to be performed once, and does not need to be repeated when changing to different biological samples. Point spread function dictionary calibration can be used for three-dimensional reconstruction of signals.

[0033] Specifically, the accuracy of the calibration point positions affects the final imaging accuracy when constructing the diffusion function dictionary for measured points. Using a transducer array self-localization method, the mapping relationship between the motion space of the three-axis precision displacement stage and the acoustic coordinate system of the arc-shaped ultrasonic transducer can be accurately established without relying on high-precision mechanical references. The specific operation procedure for constructing the diffusion function dictionary for measured points is as follows: ① Initial alignment of the spatial coordinate system and determination of the vertical transducer axis (Z-axis): Without installing the encoding mask, the carbon fiber standard point sound source is fixed at the end of the triaxial precision displacement stage and placed within the arc-shaped opening area of ​​the arc-shaped ultrasonic transducer. Based on the physical installation orientation of the triaxial precision displacement stage relative to the arc-shaped ultrasonic transducer, the sign of the Z-axis is preset in the software logic (for example, the direction away from the plane of the arc-shaped ultrasonic transducer is set as positive), thereby manually eliminating the multivalued solution space problem caused by the plane symmetry of the arc-shaped ultrasonic transducer in subsequent calculations.

[0034] ② Three-dimensional position calculation of point source based on 256-channel redundant signal: A three-axis precision displacement stage drives the standard point sound source to the scanning start position. Laser illumination of the standard point sound source generates spherical sound waves, which are synchronously received by 256 array elements. This embodiment utilizes the redundancy of the entire array data for high-precision coordinate inversion. Time delay extraction is achieved using segmented upsampling interpolation and cross-correlation algorithms to extract the flight time vectors of the sound waves arriving at each array element from the 256 channels of original radio frequency (RF) signals. The time resolution is better than 10 ns; among which This represents the flight time of the i-th path to the array element. Further, an overdetermined nonlinear equation system is constructed, using the spatial coordinates of the point sound source. The inherent electrical delay of the system Given the equivalent speed of sound c at the current water temperature as the parameter to be determined, an overdetermined nonlinear equation system is constructed, the expression of which is:

[0035] In the formula, Let be the geometric center coordinates of the i-th element, i = 1, 2, ..., 256. The Levenberg-Marquardt (LM) optimization algorithm is then used to iteratively solve the overdetermined nonlinear equations, with preset initial values ​​for the Z-axis sign used in the iterative search. Utilizing up to 256 sets of observation data, this method effectively suppresses signal-to-noise ratio fluctuations in a single channel, improving the absolute positioning accuracy of the point source in the imaging space to the 5μm-15μm level, thereby accurately calibrating the projection of the current physical position of the three-axis precision displacement stage into the acoustic coordinate system.

[0036] ③ Installing the non-uniformly perforated coded mask and performing geometric consistency verification: After completing coordinate system alignment, i.e., determining the projected coordinates of the triaxial precision displacement stage in the acoustic coordinate system, keep the arc-shaped ultrasonic transducer stationary and install the non-uniformly perforated coded mask on the front end of the arc-shaped ultrasonic transducer. To address potential minor pose deviations during coded mask installation, a geometric consistency verification is introduced: by comparing the shift of the signal envelope centroid before and after coded mask installation, the actual installation tilt angle and eccentricity of the coded mask are deduced.

[0037] ④ Automated construction of the PSF dictionary: driving the three-axis precision displacement stage according to preset... Automatic point-by-point scanning is performed with a spatial step size (0.2 mm) to traverse the entire imaging field of view. For each defined spatial position coordinate... Laser excitation records the 256-channel encoded signal after mask modulation, forming a temporal feature vector. Furthermore, all temporal feature vectors are indexed by their spatial location i to construct a five-dimensional tensor form of the measured PSF dictionary. Furthermore, during the acquisition process, abnormal waveforms caused by bubbles or extreme diffraction are automatically eliminated by calculating the phase coherence factor of each channel signal, ensuring that each impulse response function in the dictionary can truly reflect the spatial encoding characteristics of the mask for that point source. In addition, in the PSF dictionary... During the construction process, considering the limited depth of focus characteristics of the 7.5 MHz arc-shaped ultrasonic transducer along the z-axis, this embodiment not only records the original time-domain signal but also performs spatial domain normalization processing on the signal. That is, for PSFs at different depths, the energy attenuation factor is calculated and used as a pre-weighting coefficient. This storage allows the algorithm to automatically compensate for the signal-to-noise ratio loss caused by the point source deviating from the focus center during subsequent rapid 3D reconstruction, ensuring the consistency of imaging quality within the 3D field of view.

[0038] (3) Encoded signal acquisition: The sample to be tested is placed in the imaging field of view and excited by a single pulse laser. The photoacoustic signal is spatially premodulated by the non-uniform punched encoding mask of the sensing and encoding module, and the aliased signal is acquired by the data acquisition and processing unit.

[0039] Specifically, during the signal acquisition phase, the sample to be tested (e.g., the ear of a C57BL / 6 mouse) is carefully fixed on a triaxial precision displacement stage and completely immersed in deaerated water to ensure acoustic coupling. A single laser pulse irradiates the ear tissue, inducing its endogenous hemoglobin to absorb light energy and generate a three-dimensionally distributed initial pressure field. The generated photoacoustic waves, in their propagation path to the arc-shaped ultrasonic transducer, first pass through the non-uniformly perforated coded mask at the front end, where they are spatially coded and modulated by the non-uniform perforation structure. Subsequently, all 256 elements of the arc-shaped ultrasonic transducer receive and output 256 channels of aliased signals. Thanks to the spatial multiplexing properties of the coded mask, three-dimensional information within the entire imaging field of view can be captured with a single laser excitation, without the need for probe rotation or mechanical scanning. The acquired raw signal is then subjected to routine preprocessing, including 1-15MHz bandpass filtering, DC component removal, and time gain compensation (TGC), ultimately forming the input data vector for reconstruction.

[0040] (4) Fast 3D reconstruction: Based on the point spread function dictionary and the acquired aliasing signal, the accelerated Richardson-Lucy iterative algorithm with embedded spatial position correlation gain compensation weight is used to decouple and reconstruct the 3D initial pressure distribution of the sample to be tested.

[0041] Specifically, in the fast 3D image reconstruction stage, the preprocessed signal f is first compared with the aforementioned PSF dictionary. Constructing a linear observation model:

[0042] In the formula, n represents the system noise. To compensate for signal attenuation caused by the edge directivity of the curved ultrasonic transducer, a spatially related gain compensation weight is introduced during the deconvolution process. Based on this, a weighted PSF matrix is ​​constructed. :

[0043]

[0044]

[0045] In the formula, This represents the flattened response vector corresponding to the k-th iteration. express The corresponding matrix, Indicates will Convert to matrix M represents the total number of channels. To prevent the constant from being divided by zero, an improved, accelerated version of the Richardson-Lucy iterative algorithm is then used to solve the above ill-conditioned inverse problem: initializing the three-dimensional pressure distribution. If the tensor is all 1s, then in the k-th iteration, the update formula is:

[0046] In the formula, This represents the three-dimensional pressure distribution obtained after the k-th iteration. Let f denote the weighted point spread function matrix, the superscript T denotes the transpose of the matrix, and f denotes the preprocessed aliased signal. To prevent constants from being divided by zero, For a full 1 tensor, This indicates element-wise multiplication.

[0047] Furthermore, after calculating the initial pressure weight distribution of each voxel within the imaging region using an iterative algorithm, the final three-dimensional photoacoustic image is constructed and output through the following steps: The non-negative weight matrix after iterative convergence is spatially scalar mapped, and then transformed into an initial pressure distribution field reflecting the tissue's light absorption characteristics based on the photoacoustic emission coefficient. That is, by decoupling and reconstructing according to the voxel distribution, the three-dimensional initial pressure distribution is obtained, and the three-dimensional reconstruction of the three-dimensional photoacoustic image is completed.

[0048] Furthermore, after rapid 3D reconstruction, image post-processing is also included: performing 3D adaptive threshold segmentation on the reconstructed 3D data volume to suppress background noise; mapping the data to the standard Cartesian coordinate system through coordinate space resampling and interpolation methods; and generating the final 3D photoacoustic image using the maximum intensity projection algorithm or volume rendering technology.

[0049] Specifically, to eliminate background noise and subtle artifacts from mask diffraction resulting from numerical calculations, the reconstructed 3D data volume undergoes 3D adaptive threshold segmentation. A noise suppression threshold is dynamically set by calculating the standard deviation of global voxels, zeroing out spurious signals below a certain percentage (e.g., 5%-10%) of peak intensity, thereby improving the contrast between vascular structures and background tissue and suppressing background noise. Subsequently, considering that the original calculations were performed in an acoustic coordinate system with the geometric center of the arc-shaped ultrasound transducer as the origin, coordinate space resampling and interpolation methods are used to map the non-uniformly distributed calculation voxels to a standard Cartesian coordinate system to match the observation habits of conventional medical imaging. For in vivo imaging of mouse ear vessels, the system further executes the maximum intensity projection (MIP) algorithm on the processed 3D data volume (i.e., the initial 3D pressure distribution), generating 2D projection views along the three orthogonal axes X, Y, and Z to visually represent the spatial topological connections of the microvascular network. Simultaneously, volume rendering technology is used to perform ray projection processing on the 3D data volume (i.e., the 3D initial pressure distribution). By adjusting the transparency transfer function and pseudo-color mapping table, the sense of depth of the vascular branches is enhanced. Finally, the system imports the processed standard format data (such as DICOM) into a 3D visualization analysis platform, which can not only provide tomographic images of arbitrary sections, but also perform quantitative analysis of vascular diameter, branch angle, and tissue hemoglobin concentration, thus completing 3D reconstruction.

[0050] Using the above method, this invention successfully achieves high-quality three-dimensional photoacoustic imaging of blood vessels in the ear of a live mouse by efficiently spatially encoding the sound field using a front-end physical mask while retaining the full 256-channel signal acquisition capability, combined with a fast deconvolution algorithm based on a measured PSF dictionary. This method verifies the effectiveness of the mask design and fully demonstrates the good balance between system flexibility and imaging performance achieved by this invention.

[0051] The above embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit it. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A rapid three-dimensional photoacoustic imaging system, characterized in that, include: The pulse excitation unit includes a nanosecond-level pulsed laser for emitting pulsed lasers to excite biological tissues to generate broadband photoacoustic signals. The three-dimensional calibration motion unit includes a three-axis precision displacement stage, which is used to carry a standard point sound source to scan the entire imaging field of view during the system initialization calibration phase in order to obtain a dictionary of spatial non-uniform point diffusion functions characterizing the system response. The sensing and coding module includes an arc-shaped ultrasonic transducer and a non-uniform perforated coding mask that is closely attached to the front end of the arc-shaped ultrasonic transducer. The non-uniform perforated coding mask is used to spatially pre-modulate the incident sound wave so that the signals from multiple sound sources are aliased and output on a small number of receiving channels. The data acquisition and processing unit, including a multi-channel synchronous acquisition board and a host computer, is used to receive aliased signals and perform fast image reconstruction.

2. The rapid three-dimensional photoacoustic imaging system according to claim 1, characterized in that, The design method for the non-uniform perforation coding mask specifically includes: Based on the element spacing of the arc-shaped ultrasonic transducer and the radius of curvature of the probe, calculate its projection resolution unit on the mask sphere. Set a minimum center-to-center distance for the holes so that the acoustic responses produced by different hole positions can be effectively distinguished. A non-uniform perforation density distribution strategy along the axial direction is introduced to make the perforation probability of the two ends of the mask higher than that of the center region, in order to compensate for the decrease in sensitivity of the edge field of view of the arc-shaped ultrasonic transducer caused by the attenuation of directivity. Perform boundary safety checks on all candidate aperture locations and eliminate incomplete apertures that are close to the mask edge.

3. The rapid three-dimensional photoacoustic imaging system according to claim 1, characterized in that, When the data acquisition and processing unit performs fast image reconstruction, it employs a fast deconvolution reconstruction algorithm based on a dictionary of measured point spread functions. This algorithm includes: The dictionary of diffusion functions of measured points obtained during the calibration phase is used as the forward model of the system; During the reconstruction process, spatially related gain compensation weights are embedded in the deconvolution operator to dynamically enhance weak edge signals. An accelerated version of the Richardson-Lucy iterative framework is used to perform matrix decoupling operations on the acquired coded signals and the point spread function dictionary to reconstruct the three-dimensional initial pressure distribution of the object.

4. The rapid three-dimensional photoacoustic imaging system according to claim 1, characterized in that, The data acquisition and processing unit is also used for: After reconstructing the initial three-dimensional pressure distribution, a three-dimensional adaptive threshold segmentation process is performed. The noise suppression threshold is dynamically set by calculating the standard deviation of global voxels to improve the contrast between biological structures and background tissues. By using coordinate space resampling and interpolation methods, non-uniformly distributed computational voxels are mapped to the standard Cartesian coordinate system; The maximum intensity projection algorithm or volume rendering technique is applied to the processed three-dimensional initial pressure distribution to generate the final three-dimensional photoacoustic image.

5. A rapid three-dimensional photoacoustic imaging reconstruction method, applied to the system according to any one of claims 1-4, characterized in that, Includes the following steps: (1) System setup and initialization: Configure the pulse excitation unit, sensing and encoding module, data acquisition and processing unit and three-dimensional calibration motion unit to build and initialize the system; (2) Point spread function dictionary calibration: Under the condition of no biological sample, the three-dimensional calibration motion unit carries a standard point sound source to scan the entire imaging field of view and obtain the spatial non-uniform point spread function dictionary characterizing the system response; (3) Encoded signal acquisition: The sample to be tested is placed in the imaging field of view, and the photoacoustic signal is spatially premodulated by the non-uniform perforation encoding mask of the sensing and encoding module through single pulse laser excitation, and the aliasing signal is acquired by the data acquisition and processing unit. (4) Fast 3D reconstruction: Based on the point spread function dictionary and the acquired aliasing signal, the accelerated Richardson-Lucy iterative algorithm with embedded spatial position correlation gain compensation weight is used to decouple and reconstruct the 3D initial pressure distribution of the sample to be tested.

6. The reconstruction method for rapid three-dimensional photoacoustic imaging according to claim 5, characterized in that, The point spread function dictionary calibration specifically includes: The standard point sound source is fixed at the end of the three-axis precision displacement stage, and the Z-axis of the arc-shaped ultrasonic transducer is determined according to the physical installation orientation of the three-axis precision displacement stage relative to the arc-shaped ultrasonic transducer. Without installing the encoding mask, the mapping relationship between the motion space of the triaxial precision displacement stage and the acoustic coordinate system of the arc-shaped ultrasonic transducer is calibrated using a transducer array self-localization method, so as to determine the projection of the triaxial precision displacement stage in the acoustic coordinate system. Install a non-uniform perforation coding mask and perform geometric consistency verification on the pose deviation during the installation process; The three-axis precision displacement stage is driven to perform automatic point-by-point scanning according to the preset spatial step size to traverse the entire imaging field of view and record the multi-channel encoded signal after mask modulation corresponding to each spatial position to form a time-domain feature vector. All time-domain feature vectors are indexed by spatial position to construct a dictionary of spread functions for the measured points.

7. The reconstruction method for rapid three-dimensional photoacoustic imaging according to claim 6, characterized in that, The method based on transducer array self-localization specifically includes: By utilizing the redundancy of the entire array data, the flight time vector of the sound wave arriving at each array element is extracted from the original multi-channel radio frequency signal. An overdetermined nonlinear equation set is constructed using the spatial coordinates of the point sound source, the inherent electrical delay of the system, and the equivalent sound velocity at the current water temperature as the parameters to be determined. The Levenberg-Marquardt optimization algorithm was used to iteratively solve the overdetermined nonlinear equations to obtain the projected coordinates of the triaxial precision displacement stage in the acoustic coordinate system.

8. The reconstruction method for rapid three-dimensional photoacoustic imaging according to claim 5, characterized in that, The iterative update formula for the accelerated Richardson-Lucy iterative algorithm with embedded spatial location-related gain compensation weights is as follows: In the formula, This represents the three-dimensional pressure distribution obtained after the k-th iteration. Let f denote the weighted point spread function matrix, the superscript T denotes the transpose of the matrix, and f denotes the preprocessed aliased signal. To prevent constants from being divided by zero, For a full 1 tensor, This indicates element-wise multiplication.

9. The reconstruction method for rapid three-dimensional photoacoustic imaging according to claim 5, characterized in that, Following the rapid 3D reconstruction, image post-processing is also included: The reconstructed 3D data volume is subjected to 3D adaptive threshold segmentation to suppress background noise. Data is mapped to the standard Cartesian coordinate system using coordinate space resampling and interpolation methods. The final 3D photoacoustic image is generated using the maximum intensity projection algorithm or volume rendering technology.

Images