An x-ray thermoacoustic imaging method based on diffusion model and sound velocity compensation

Abstract

The application discloses an X-ray thermoacoustic imaging method based on a diffusion model and sound velocity compensation and relates to the technical fields of X-ray thermoacoustic imaging and deep learning imaging reconstruction. The application fully considers the influence of the non-uniformity of sound velocity on the thermoacoustic wave propagation path and time delay, and proposes a new method of combining sound velocity compensation and diffusion model reconstruction; under the premise of not increasing the hardware complexity of the imaging system, the sound velocity modeling is used to realize the accurate correction of the wave propagation path, and the imaging artifacts and spatial positioning errors are effectively avoided; meanwhile, the image inversion framework based on the conditional diffusion model has excellent noise suppression and high-frequency detail recovery capability, and can generate high-quality thermoacoustic images with clear structure and significant contrast under imaging conditions such as low signal-to-noise ratio and low energy dose.

Description

Technical Field

[0001] This invention relates to the field of X-ray thermoacoustic imaging and deep learning imaging reconstruction technology, specifically to an X-ray thermoacoustic imaging method based on a diffusion model and sound velocity compensation. Background Technology

[0002] X-ray-induced thermoacoustic imaging is a hybrid imaging method that combines the penetrability of X-ray imaging with the high resolution of ultrasound imaging. Its basic principle is to irradiate tissue with short-pulse X-rays. The thermoelastic expansion generated by the interaction between the X-rays and the tissue induces localized ultrasound signals. These signals are received by an external ultrasound sensor array, and image reconstruction is performed, thereby enabling visualization of the internal structure of the tissue.

[0003] This method combines the high penetrability of X-rays with the high-contrast imaging capability of ultrasound for soft tissues, thus showing potential in areas such as tissue lesion detection and combined imaging of bone and soft tissues. However, current X-ray thermoacoustic imaging still faces many challenges in terms of imaging resolution, contrast, and reconstruction accuracy.

[0004] On the one hand, the spatial distribution of sound velocity in biological tissues is usually not uniform. The resulting sound velocity mismatch leads to deviations in the propagation path of thermoacoustic signals and time delay errors, ultimately causing spatial distortion and blurring of the image. On the other hand, traditional thermoacoustic image reconstruction methods (such as time-delay stacking or model inversion methods) have limited ability to model sound velocity distribution, making it difficult to accurately recover tissue boundaries and detailed structures. In addition, the limited X-ray irradiation dose and low signal-to-noise ratio also affect image quality.

[0005] In recent years, diffusion models have emerged as a breakthrough in generative deep learning models, demonstrating powerful image modeling and reconstruction capabilities. However, their application in the field of thermoacoustic imaging is still in the exploratory stage, particularly in how to integrate them with physical priors (such as sound velocity distribution) to enhance image reconstruction quality, which still lacks systematic research and engineering implementation paths.

[0006] Therefore, a new solution is needed to address the above problems. Summary of the Invention

[0007] The purpose of this invention is to provide an X-ray thermoacoustic imaging method based on a diffusion model and sound velocity compensation, which addresses the problems of image blurring, severe reconstruction artifacts, and insensitivity to changes in tissue sound velocity in existing X-ray thermoacoustic imaging methods.

[0008] To achieve the above objectives, the present invention provides the following technical solution: an X-ray thermoacoustic imaging method based on a diffusion model and sound velocity compensation, comprising at least the following steps:

[0009] S1: X-ray thermoacoustic signal simulation data generated based on joint modeling of MCX and K-Wave;

[0010] S2: Construct a diffusion model for learning the distribution of X-ray thermoacoustic signals, wherein the diffusion model simulates the distribution evolution of thermoacoustic signals through forward and reverse diffusion processes;

[0011] S3: Construct a sound velocity compensation module to incorporate the sound velocity information of the tissue region into the delay summation reconstruction process to correct the sound wave propagation time and improve the accuracy of the reconstructed image;

[0012] S4: Train the joint network of the diffusion model and the sound speed compensation module using the training dataset;

[0013] S5: Perform domain transformation on the X-ray thermoacoustic signals acquired under sparse viewpoints, and use the trained network to perform image completion and enhancement to obtain high-quality X-ray thermoacoustic image reconstruction results.

[0014] Furthermore, S1 includes at least the following steps:

[0015] Based on CT slices or standard anatomical images, a three-dimensional tissue structure model of the target area is constructed, and different tissue types are classified.

[0016] Based on the tissue information in the CT images, optical parameters for each tissue are set. The optical parameters include absorption coefficient μa, scattering coefficient μs, anisotropy coefficient g, refractive index n, and thermoacoustic parameters, including density ρ, Grunison coefficient Γ, and sound velocity c.

[0017] In the tissue model, the incident direction, energy spectrum and photon number of the X-ray source are set, and the energy deposition distribution H(x,y,z) in three-dimensional space is calculated using Monte Carlo simulation (MCX).

[0018] Based on the thermoelastic expansion theory, the energy deposition H(x,y,z) is converted into the initial sound pressure distribution p0(x,y,z), and the conversion formula is as follows:

[0019]

[0020] Where H is the deposition energy per unit volume, and Γ(x,y,z) is the Gronnison coefficient of the tissue;

[0021] Based on the distribution of tissues in CT images, a corresponding sound velocity value is assigned to each voxel to construct a non-uniform sound velocity field, and simulation parameters are configured, including spatial resolution, ultrasonic sensor array arrangement, and sensor parameters.

[0022] The spatial resolution of the image is preset, and the number, shape, and spatial distribution of the sensor array are defined in K-Wave; the key parameters of the ultrasonic sensor are set, including the center frequency, bandwidth, and sampling frequency.

[0023] The initial sound pressure distribution p0(x,y,z) is used as the initial condition for K-Wave acoustic simulation to simulate the propagation of thermoacoustic waves in a non-uniform medium, and noise is added to improve the simulation realism.

[0024] Add simulated noise with statistical distribution characteristics (such as Gaussian white noise) to approximate the signal collected by the real detector;

[0025] Under different sampling conditions (such as the number of viewpoints, sensor density, and bandwidth configuration), the simulation was repeated to generate multiple sets of ultrasonic signal data, and sparse and fully sampled data pairs were constructed to serve as the training input and target for training the image reconstruction model.

[0026] Furthermore, the diffusion model is used for: denoising the received signal and learning the logarithmic gradient of the signal distribution. To optimize the image reconstruction process;

[0027] The forward diffusion process of the diffusion model in S2 satisfies the following stochastic differential equation:

[0028] dx=f(x,t)dt+g(t)dw

[0029] in, It is the drift coefficient; It is the diffusion coefficient; dw represents Brownian motion; dw represents Brownian increment motion.

[0030] The reverse diffusion process of the diffusion model satisfies:

[0031]

[0032] in It is an infinitesimal time step. It is an inverse Brownian motion;

[0033] The diffusion model is trained using a score network. To estimate the fractional function of time step t The parameter optimization objective is:

[0034]

[0035] in, This represents the squared error between the output of the fractional network and the true fractional function. It is time The weighting factor is used to control the impact of different time steps; It is the expectation operation, which represents different random variables. Perform expectation calculation; This is also the desired operation.

[0036] Furthermore, the sound velocity compensation method of the compensation module in S3 is as follows: calculating the propagation path of the sound wave in the non-uniform medium. Calculate time delay using path integral :

[0037]

[0038] Where c(x,y) represents the speed of sound at position (x,y); dl represents the arc length along the propagation path Γ;

[0039] The sound velocity compensation method uses the finite difference method or the finite element method to simulate sound wave propagation in order to calculate the sound wave delay compensation parameters in each tissue layer.

[0040] Furthermore, S5 includes at least the following steps:

[0041] Construct a dataset containing the original images and corresponding thermoacoustic signals;

[0042] An end-to-end training approach is adopted, with the difference between the reconstructed image and the original image as the optimization objective;

[0043] By combining perceptual loss, structural similarity loss, and mean squared error loss, the fractional network parameters and sound velocity compensation parameters are jointly optimized.

[0044] Compared with the prior art, the beneficial effects of the present invention are:

[0045] This invention fully considers the impact of sound velocity non-uniformity on the propagation path and time delay of thermoacoustic waves, and proposes a novel method combining sound velocity compensation and diffusion model reconstruction. Without increasing the hardware complexity of the imaging system, sound velocity modeling is used to accurately correct the wave propagation path, effectively avoiding imaging artifacts and spatial positioning errors. Simultaneously, the image inversion framework based on the conditional diffusion model possesses excellent noise suppression and high-frequency detail recovery capabilities, enabling the generation of high-quality thermoacoustic images with clear structure and significant contrast under imaging conditions such as low signal-to-noise ratio and low energy dose. Furthermore, the method of this invention can integrate multimodal images (such as CT and MRI) as reconstruction conditions, improving tissue structure recognition and enhancing the system's adaptability and robustness. This provides technical support and engineering feasibility for the practical application of X-ray thermoacoustic imaging in early lesion detection, tumor localization, and material testing. This invention will provide technical support for the implementation of X-ray thermoacoustic imaging in clinical disease early screening, tissue damage assessment, and non-destructive testing of materials. Attached Figure Description

[0046] To more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0047] Figure 1 This is a schematic diagram illustrating the comparative effect of the sound speed compensation of the present invention.

[0048] Figure 2 This invention uses USCT to obtain the sound velocity distribution map of the target tissue;

[0049] Figure 3 This is an overall flowchart of the invention based on the diffusion model and sound velocity compensation;

[0050] Figure 4 This is a comparison chart of the iterative reconstruction effects of the method proposed in this invention. Detailed Implementation

[0051] The technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments.

[0052] See Figures 1-4 An X-ray thermoacoustic imaging method based on a diffusion model and sound velocity compensation includes at least the following steps:

[0053] S1: Simulation data of X-ray thermoacoustic signals generated based on joint modeling of MCX and K-Wave:

[0054] The contents described in S1 include the following:

[0055] A three-dimensional tissue structure model of the target area is established based on CT slices or standard anatomical images, and different tissue regions are divided. Based on different tissue types (such as bone, soft tissue, and adipose tissue), specific optical parameters are assigned, including absorption coefficient (μa), scattering coefficient (μs), anisotropy coefficient (g), and refractive index (n). The incident direction, energy spectrum, and photon number of the X-ray source are set in the tissue model. Monte Carlo simulation of the photon propagation path in the tissue is performed using MCX software to obtain the energy deposition distribution map of each voxel in three-dimensional space. Based on the thermoelastic expansion formula, the energy deposition distribution obtained from the MCX simulation is converted into an initial sound pressure distribution map, where the conversion parameters include tissue density and the Glunison coefficient. The sound pressure conversion relationship is as follows:

[0056]

[0057] Where H is the deposition energy per unit volume, and Γ is the Glunison coefficient, which can be obtained by looking up a table or fitting based on tissue-specific parameters. Based on the tissue distribution in the CT image, a corresponding sound velocity value is assigned to each voxel to construct a non-uniform sound velocity field; the spatial resolution of the image is preset, and the number, shape, and spatial distribution of the sensor array are defined in K-Wave; key parameters of the ultrasound sensor are set, including center frequency, bandwidth, and sampling frequency; the initial sound pressure distribution is set as the initial condition in the K-Wave simulation to simulate the propagation process of thermoacoustic waves in a non-uniform medium; simulated noise with statistical distribution characteristics (such as Gaussian white noise) is added to approximate the signal collected by the real detector; the above simulation process is repeated under different sampling conditions (such as the number of viewing angles, sensor density, and bandwidth configuration), outputting multiple ultrasound simulation signal datasets; sparse and non-sparse data pairs are constructed as training inputs and targets for subsequent image reconstruction models.

[0058] S2: Construct a diffusion model for learning the distribution of X-ray thermoacoustic signals, wherein the diffusion model simulates the distribution evolution of thermoacoustic signals through forward and reverse diffusion processes;

[0059] like Figure 3 As shown, in Figure 3 In the upper part, specifically, the initial X-ray thermoacoustic signal sine wave is first processed. Nearest neighbor interpolation is used to reshape the sine wave into a symmetrical size that is more suitable for convolution operations. Then, a mask is applied to the reshaped sine wave to extract valid data. The mask operation is the core of sparse data reconstruction. The mask is a two-dimensional matrix with a value of 0 or 1 and the same shape as the sine wave. The corresponding mask matrix can be defined according to the required sampling angle or sensor array settings. 1 indicates that the data at that position is sampled, and 0 indicates that the data at that position is discarded.

[0060] To ensure consistency between the reconstructed signal and the target data sample, the model employs data consistency and regularization constraints for optimization. During the optimization process, the data consistency term is defined as:

[0061]

[0062] in The extraction mask for the sine wave is used to determine and differentiate sampling strategies. This represents the sine wave of the current reconstruction. This indicates that the sine wave reconstructed through interpolation, when constrained by data consistency during optimization reconstruction, can be measured by the L2 norm. and The error between them is minimized to ensure that the generated sine wave data and the original sparse data remain consistent at the sampling points during the optimization process.

[0063] Since sparse sampling leads to an ill-posed optimization problem, a regularization term needs to be introduced. This can alleviate the pathological nature of sparse reconstruction.

[0064]

[0065] in This is the regularization coefficient, used to balance the weighting of data consistency and regularization constraints. The regularization function is derived from high-quality prior knowledge obtained by the noisy gradient scoring network generated by the diffusion model, used to generate more reliable data. Therefore, the final optimization objective is:

[0066]

[0067] This invention achieves iterative reconstruction of the sine wave through two alternating update operations: prediction and correction. The prediction aims to generate a preliminary estimate of the samples at the current time step, attempting to approach the target distribution using current gradient information and random perturbations. The target signal data is generated from the prior distribution using the following formula:

[0068]

[0069] in, It is a noise scale. It is the number of iterations in the model solution. It is Gaussian white noise that conforms to a standard normal distribution.

[0070] Samples in the current state It will be based on the score network The judgment shifts in the direction of gradient ascent, adjusting to a more probable position. The adjustment range was determined, and after adjustment, some random noise was added. By introducing randomness, the prediction results are made more consistent with the true distribution while exhibiting diversity.

[0071] The prediction operation provides the direction for sample adjustment, but it is not accurate. Therefore, a correction operation is needed to further refine the samples, making them more consistent with the target probability distribution. The Langevin Markov correction algorithm is used for further optimization, as shown in the following equation:

[0072] First, the predicted value obtained from the prediction operation will be... Based on this, use the gradient information provided by the model. Further precise adjustments were made, ε i The adjustment step size is determined, moving closer to regions with higher probability distributions of the target. Furthermore, some random perturbations are incorporated during the adjustment process. Ensure the diversity of generated samples.

[0073] In the prediction and correction phase, according to The design alternately performs replacement fidelity operations. During the iterative reconstruction process, each intermediate result generated is affected by the fidelity item. In the result generated in the current iteration, a portion of the data is replaced by the data in the original sine graph, thereby ensuring data consistency.

[0074] S3: Construct a sound velocity compensation module to incorporate the sound velocity information of the tissue region into the delay summation reconstruction process to correct the sound wave propagation time and improve the accuracy of the reconstructed image;

[0075] like Figure 3 As shown, in Figure 3 In the lower part, specifically, after obtaining the iteratively reconstructed sine wave in S2, the Delayed Summation (DAS) method is used to transform the sine wave data into the image domain. However, the traditional DAS method assumes that the sound velocity c is uniform throughout the medium. This assumption may lead to artifacts and reduced resolution in the reconstruction results when the sound velocity distribution is complex. Considering the practical factors of X-ray thermoacoustic imaging, the sound velocity distribution information c(x,y) is embedded in the DAS method. Non-uniform sound velocity is used to correct the signal propagation time, thereby further improving the accuracy of image reconstruction. In the traditional DAS method, assuming that the sound velocity is uniform, each signal S in the sine wave... m (t m The time delay t) m Based on the detector position r m (x m ,y m ) and reconstruct the target point location r p The Euclidean distance between (x, y) is used to determine:

[0076]

[0077] When the heterogeneity of sound speed is embedded in DAS reconstruction, the sound speed is not uniformly distributed. The sound speed at position (x,y) can then be expressed as c = c(x,y), and the propagation time t... m The calculation can be performed using path integrals:

[0078]

[0079] in Indicates the signal originating from the sound pressure level. propagation to detector location The true transmission path, Indicates the path along the signal propagation path The minute arc length element.

[0080] S4: Train the joint network of the diffusion model and the sound speed compensation module using the training dataset;

[0081] S5: Perform domain transformation on the X-ray thermoacoustic signals acquired under sparse viewpoints, and use the trained network to perform image completion and enhancement to obtain high-quality X-ray thermoacoustic image reconstruction results.

[0082] like Figure 4 As shown, specifically, in the reconstruction process of this method from the 1st iteration to the 500th iteration, with the increase of the number of iterations, the noise in the X-ray thermoacoustic signal sine wave is gradually eliminated, and the reconstructed image gradually becomes clearer. After the 500th iteration, the method of this invention has basically recovered all the X-ray thermoacoustic signal information, and has also obtained better reconstructed image results, further improving the image quality.

[0083] Further clarification of the image reconstruction process includes:

[0084] (1) Map the sparsely sampled thermoacoustic signal to the potential space;

[0085] (2) Image completion and enhancement are performed using a trained backdiffusion network;

[0086] (3) Combine sound speed compensation guidance to make the reconstructed image approximate the real distribution in the latent space;

[0087] (4) Output high-resolution X-ray thermoacoustic images.

[0088] In summary:

[0089] This invention utilizes a generative model to learn prior information from X-ray thermoacoustic signals. By learning the prior distribution of X-ray thermoacoustic signals sampled under non-sparse conditions, it progressively guides and optimizes X-ray thermoacoustic signals from sparse perspectives, effectively overcoming the signal sparsity problem in sparse reconstruction. It applies a domain transformation to the X-ray thermoacoustic signals using the Delayed Summation (DAS) method, mapping the signals to the image space. During the domain transformation using DAS, sound velocity compensation further optimizes the reconstruction quality. An image reconstruction network based on a diffusion probability model is constructed, introducing noise during the diffusion process to simulate signal degradation. A reverse diffusion process then progressively restores the true thermoacoustic image, thereby improving the detail and contrast of the reconstructed image. The diffusion model is trained using a training dataset containing X-ray thermoacoustic images from both full and sparse perspectives, enhancing its ability to recover images under different sampling conditions. This invention enables high-quality reconstruction and enhancement of X-ray thermoacoustic images even under conditions of limited acquisition perspective or poor signal quality, improving the robustness and imaging effect of the imaging system and promoting the practical application of X-ray thermoacoustic imaging technology in biomedical diagnosis, tissue functional imaging, and other fields.

[0090] It will be apparent to those skilled in the art that the present invention is not limited to the details of the exemplary embodiments described above, and that the invention can be implemented in other specific forms without departing from its spirit or essential characteristics. Therefore, the embodiments should be considered in all respects as exemplary and non-limiting, and the scope of the invention is defined by the appended claims rather than the foregoing description. Thus, all variations falling within the meaning and scope of equivalents of the claims are intended to be included within the present invention. No reference numerals in the claims should be construed as limiting the scope of the claims.

Claims

1. An X-ray thermoacoustic imaging method based on a diffusion model and sound velocity compensation, characterized in that: At least the following steps are included: S1: X-ray thermoacoustic signal simulation data generated based on Monte Carlo simulation of MCX and K-Wave joint modeling; S2: Construct a diffusion model for learning the distribution of X-ray thermoacoustic signals, wherein the diffusion model simulates the distribution evolution of thermoacoustic signals through forward and reverse diffusion processes; The diffusion model is used for: denoising the received signal and learning the logarithmic gradient of the signal distribution. To optimize the image reconstruction process; The forward diffusion process of the diffusion model in S2 satisfies the following stochastic differential equation: dx=f(x,t)dt+g(t)dw in, It is the drift coefficient; It is the diffusion coefficient; dw represents Brownian motion; dw represents Brownian increment motion. The reverse diffusion process of the diffusion model satisfies: Where dt is an infinitesimal time step. It is an inverse Brownian motion; The diffusion model is trained using a score network S. θ (x t The fractional function (,t) is used to estimate the time step t. Its parameter optimization objective is: in, λ(t) represents the squared error between the output of the fractional network and the true fractional function; λ(t) is the weighting factor for time t, used to control the effect of different time steps. It is the expectation operation, representing the expectation operation for different random variables x0 and x... t Perform expectation calculation; E t～N(0,T) Also for the expected operation; S3: Construct a sound velocity compensation module to incorporate the sound velocity information of the tissue region into the delay summation reconstruction process to correct the sound wave propagation time and improve the accuracy of the reconstructed image; The sound velocity compensation method of the compensation module in S3 is as follows: calculate the propagation path Γ of the sound wave in the non-uniform medium, and calculate the time delay t through path integral. m : Where c(x,y) represents the speed of sound at position (x,y); dl represents the arc length along the propagation path Γ; The sound velocity compensation method uses the finite difference method or the finite element method to simulate sound wave propagation in order to calculate the sound wave delay compensation parameters in each tissue layer. S4: Train the joint network of the diffusion model and the sound speed compensation module using the training dataset; S5: Perform domain transformation on the X-ray thermoacoustic signals acquired under sparse viewpoints, and use the trained network to perform image completion and enhancement to obtain high-quality X-ray thermoacoustic image reconstruction results.

2. The X-ray thermoacoustic imaging method based on diffusion model and sound velocity compensation according to claim 1, characterized in that: S1 includes at least the following steps Based on CT slices or standard anatomical images, a three-dimensional tissue structure model of the target area is constructed, and different tissue types are classified. Based on the tissue information in the CT images, optical parameters for each tissue are set. The optical parameters include absorption coefficient μa, scattering coefficient μs, anisotropy coefficient g, refractive index n, and thermoacoustic parameters, including density ρ, Grunison coefficient Γ, and sound velocity c. In the tissue model, the incident direction, energy spectrum and photon number of the X-ray source are set, and the energy deposition distribution H(x,y,z) in three-dimensional space is calculated using Monte Carlo simulation MCX. Based on the thermoelastic expansion theory, the energy deposition H(x,y,z) is converted into the initial sound pressure distribution p0(x,y,z), and the conversion formula is as follows: p0(x,y,z)=Γ(x,y,z)·H(x,y,z) Where H is the deposition energy per unit volume, and Γ(x,y,z) is the Gronnison coefficient of the tissue; Based on the distribution of tissues in CT images, a corresponding sound velocity value is assigned to each voxel to construct a non-uniform sound velocity field, and simulation parameters are configured, including spatial resolution, ultrasonic sensor array arrangement, and sensor parameters. The spatial resolution of the image is preset, and the number, shape, and spatial distribution of the sensor array are defined in K-Wave; the key parameters of the ultrasonic sensor are set, including the center frequency, bandwidth, and sampling frequency. The initial sound pressure distribution p0(x,y,z) is used as the initial condition for K-Wave acoustic simulation to simulate the propagation of thermoacoustic waves in a non-uniform medium, and noise is added to improve the simulation realism. Add simulated noise with statistical distribution characteristics to approximate the signal collected by the real detector; Under different sampling conditions, the simulation was repeated to generate multiple sets of ultrasonic signal data, and sparse and fully sampled data pairs were constructed to serve as the training input and target for training the image reconstruction model.

3. The X-ray thermoacoustic imaging method based on diffusion model and sound velocity compensation according to claim 1, characterized in that: The S5 includes at least the following steps: Construct a dataset containing the original images and corresponding thermoacoustic signals; An end-to-end training approach is adopted, with the difference between the reconstructed image and the original image as the optimization objective; By combining perceptual loss, structural similarity loss, and mean squared error loss, the fractional network parameters and sound velocity compensation parameters are jointly optimized.

Images