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

Abstract

The invention discloses an X-ray thermoacoustic imaging method based on a diffusion model and sound velocity compensation, and relates to the technical field of X-ray thermoacoustic imaging and deep learning imaging reconstruction. According to the method, the influence of sound velocity non-uniformity on a thermal sound wave propagation path and time delay is fully considered, and a new method combining sound velocity compensation and diffusion model reconstruction is provided; on the premise that the hardware complexity of an imaging system is not increased, accurate correction of a wave propagation path is achieved through sound velocity modeling, and imaging artifacts and space 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 capabilities, and can generate a high-quality thermoacoustic image with a clear structure and remarkable contrast under the imaging conditions of low signal-to-noise ratio, low energy dose and the like.

Description

Technical Field

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

[0002] X-ray-induced thermoacoustic imaging is a hybrid imaging method that combines the penetrating properties of X-ray imaging with the high resolution of ultrasound imaging. Its basic principle is to illuminate tissue with short pulses of X-rays. The interaction between the X-rays and the tissue produces thermoelastic expansion, which induces localized ultrasonic signals. These signals are received by an external ultrasonic transducer array and image reconstructed, thereby visualizing the internal structure of the tissue.

[0003] This method combines the high penetrability of X-rays with the high-contrast soft tissue imaging capabilities of ultrasound, and therefore has potential for tissue lesion detection and combined bone and soft tissue imaging. However, X-ray thermoacoustic imaging currently faces many challenges in imaging resolution, contrast, and reconstruction accuracy.

[0004] On the one hand, the spatial distribution of sound velocity in biological tissue is typically uneven. The resulting sound velocity mismatch can lead to deviations in the propagation path of the thermoacoustic signal 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-delayed superposition or model inversion) have limited ability to model the sound velocity distribution, making it difficult to accurately restore tissue boundaries and detailed structures. Furthermore, due to the limited X-ray dose, the signal-to-noise ratio is low, which also affects imaging quality.

[0005] In recent years, diffusion models, a breakthrough in generative deep learning models, have demonstrated powerful image modeling and reconstruction capabilities. However, their application in thermoacoustic imaging is still exploratory. In particular, systematic research and engineering implementation are still lacking on how to integrate them with physical priors (such as the sound velocity distribution) to enhance image reconstruction quality.

[0006] Therefore, a new solution to the above problems needs to be proposed. Summary of the Invention

[0007] The purpose of the present invention is to provide an X-ray thermoacoustic imaging method based on diffusion model and sound velocity compensation, so as to solve the problems of image blur, serious reconstruction artifacts and insensitivity to tissue sound velocity changes in existing X-ray thermoacoustic imaging.

[0008] To achieve the above object, 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: Generate X-ray thermoacoustic signal simulation data based on MCX and K-Wave joint modeling;

[0010] S2: Constructing 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 diffusion and reverse diffusion processes;

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

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

[0013] S5: Perform domain transformation on the X-ray thermoacoustic signals collected under sparse viewing angles, use the trained network to complete and enhance the images, and obtain high-quality X-ray thermoacoustic image reconstruction results.

[0014] Furthermore, the S1 at least includes 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 divided;

[0016] According to the tissue information in the CT image, optical parameters of each tissue are set, wherein the optical parameters include absorption coefficient μa, scattering coefficient μs, anisotropy coefficient g, refractive index n, and thermoacoustic parameters, wherein the thermoacoustic parameters include density ρ, Grünitz coefficient Γ, and speed of sound c;

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

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

[0019] p0(x,y,z)=Γ(x,y,z)·H(x,y,z)

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

[0021] According to the distribution of tissue in the CT image, a corresponding sound velocity value is assigned to each voxel, a non-uniform sound velocity field is constructed, and simulation parameters are configured, including spatial resolution, ultrasonic sensor array arrangement, and sensor parameters;

[0022] Preset the spatial resolution of the image and define the number, shape, and spatial distribution of the sensor array in K-Wave; set key parameters of the ultrasonic sensor, including 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 thermal acoustic waves in inhomogeneous media, 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 viewing angles, sensor arrangement density, and bandwidth configuration), the simulation is repeated to generate multiple sets of ultrasonic signal data, and sparse and fully sampled data pairs are constructed to serve as training inputs and targets for training the image reconstruction model.

[0026] Furthermore, the diffusion model is used to: denoise the received signal and learn 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, is the drift coefficient; is the diffusion coefficient; represents Brownian motion; dw represents Brownian incremental motion;

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

[0031]

[0032] where dt is an infinitesimal time step, is the reverse Brownian motion;

[0033] The diffusion model is constructed by training the score network S θ (x t ,t) to estimate the score function of time step t The parameter optimization goal is:

[0034]

[0035] in, represents the square error between the score network output and the true score function; λ(t) is the weighting factor at time t, which is used to control the influence of different time steps; It is the expectation operation, which means that for different random variables x0 and x t Perform expectation calculation; t～N(0,T) Also for expected operation.

[0036] Furthermore, the sound velocity compensation method of the compensation module in S3 is: calculating the propagation path Γ of the sound wave in the inhomogeneous medium, and calculating the time delay t by path integration. m :

[0037]

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

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

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

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

[0042] An end-to-end training method is used, with the difference between the reconstructed image and the original image as the optimization target;

[0043] The perceptual loss, structural similarity loss and mean square error loss are combined to jointly optimize the score network parameters and sound speed compensation parameters.

[0044] Compared with the prior art, the present invention has the following beneficial effects:

[0045] The present invention fully considers the impact of sound velocity non-uniformity on the propagation path and time delay of thermoacoustic waves, and proposes a new method that combines 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; at the same time, the image inversion framework based on the conditional diffusion model has excellent noise suppression and high-frequency detail recovery capabilities, 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; in addition, the method of the present invention can fuse multimodal images (such as CT, MRI, etc.) as reconstruction conditions, improve tissue structure recognition, enhance the adaptability and robustness of the system, and provide technical support and engineering feasibility for the practical application of X-ray thermoacoustic imaging in scenarios such as early lesion detection, tumor localization, and material testing; the present invention will provide technical support for the implementation of X-ray thermoacoustic imaging in clinical disease early screening, tissue damage assessment, and non-destructive material testing. BRIEF DESCRIPTION OF THE DRAWINGS

[0046] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the following briefly introduces the drawings required for describing the embodiments. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without creative work.

[0047] Figure 1 Schematic diagram of the comparative effect experiment of the sound velocity compensation of the present invention;

[0048] Figure 2 The present invention adopts USCT to obtain the target tissue sound velocity distribution map;

[0049] Figure 3 This is an overall flow chart of the present 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 claim 4. DETAILED DESCRIPTION

[0051] The technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, rather than all the embodiments.

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

[0053] S1: Generate X-ray thermoacoustic signal simulation data based on MCX and K-Wave joint modeling:

[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 tissue, soft tissue, and adipose tissue), specific optical parameters are assigned to them, including the 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 light source are set in the tissue model. The propagation path of photons in the tissue is simulated by Monte Carlo simulation 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 by MCX simulation is converted into an initial acoustic pressure distribution map, where the conversion parameters include tissue density and Grünitz coefficient. The acoustic pressure conversion relationship is:

[0056] p0(x,y,z)=Γ(x,y,z)·H(x,y,z)

[0057] Where H is the deposited energy per unit volume, and Γ is the Grünigssen coefficient, which can be obtained by looking up or fitting tissue-specific parameters. Based on the distribution of tissue 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.

[0058] The initial sound pressure distribution is set as the initial condition in the K-Wave simulation to simulate the propagation process of thermal acoustic waves in inhomogeneous media; 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 arrangement density, and bandwidth configuration) to output multiple ultrasonic simulation signal data sets; sparse and non-sparse data pairs are constructed as training input and targets for the subsequent image reconstruction model.

[0059] S2: Constructing 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 diffusion and reverse diffusion processes;

[0060] like Figure 3 As shown, in Figure 3 Specifically, in the upper part of the figure, the initial X-ray thermoacoustic signal sinusoidal graph is first processed, and the sinusoidal graph is reshaped into a symmetrical size that is more suitable for convolution operations using nearest neighbor interpolation. Then, a mask is applied to the reshaped sinusoidal graph 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 its shape is the same as that of the sinusoidal graph. The corresponding mask matrix can be defined according to the required sampling perspective or sensor array setting. 1 indicates that the data at this position is sampled, and 0 indicates that the data at this position is discarded.

[0061] In order to ensure that the reconstructed signal is consistent with the target data sample, the model uses data consistency and regularization constraints to optimize the reconstruction. During the optimization process, the data consistency term is defined as:

[0062]

[0063] Where P(Λ) represents the extraction mask of the sinusoidal graph, which is used to determine and distinguish the sampling strategy, x represents the current reconstructed sinusoidal graph, Represents the sinusoidal graph reshaped by interpolation. The data consistency constraint is used in the optimization reconstruction to measure the L2 norm of x and The error between them is reduced to ensure that during the optimization process, the generated sinusoidal data is consistent with the original sparse data in terms of sampling points.

[0064] Since sparse sampling makes the optimization problem ill-posed, a regularization term needs to be introduced. To alleviate the pathological nature of sparse reconstruction.

[0065]

[0066] Where τ is the regularization coefficient, which is used to balance the consistency of data and the weight of regularization constraints. R(x,y) is the regularization function, which is the high-quality prior knowledge obtained by the noisy gradient score network generated by the diffusion model, and is used to generate more reliable data. Therefore, the final optimization goal is:

[0067]

[0068] The present invention achieves iterative reconstruction of the sinusoidal graph through two alternating update operations: prediction and correction. The prediction goal is to predict the initial estimate of the sample at the current time step, and try to approach the target distribution through the current gradient information and random perturbations. The target signal data is generated from the prior distribution by the following formula:

[0069]

[0070] Among them, σ i is the noise scale, i is the number of model solution iterations, and z is Gaussian white noise that conforms to the standard normal distribution.

[0071] Sample x at the current state i Will be based on the score network S θ (x i ,σ i+1 ) is shifted towards the direction of gradient ascent and adjusted to a more likely position. Determines the amplitude of the adjustment. After the adjustment, a portion of random noise is added Introducing randomness provides diversity while ensuring that the prediction results are more consistent with the real distribution.

[0072] The prediction operation gives the direction of sample adjustment, but it is not accurate. Therefore, a correction operation is needed to further correct the sample so that the sample is more consistent with the target probability distribution. The Langevin Markov correction algorithm is used for further optimization. The following formula is used:

[0073]

[0074] First, the predicted value obtained in the prediction operation Based on this, we use the gradient information provided by the model For further fine-tuning, ε iThe adjustment step size is determined to move closer to the target distribution area with higher probability. In addition, some random perturbations are added during the adjustment process. Ensure the diversity of generated samples.

[0075] In the prediction and correction stages, the replacement fidelity operation is performed alternately according to formula (6). During the iterative reconstruction process, the intermediate results generated each time will be affected by the fidelity term. Part of the data in the result generated by the current iteration will be replaced by the data in the original sinusoidal diagram to ensure data consistency.

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

[0077] like Figure 3 As shown, in Figure 3 In the lower half of the , specifically, after obtaining the iteratively reconstructed sinusoidal graph in S2, the delayed sum (DAS) method is used to convert the sinusoidal graph 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 of the reconstruction results when the sound velocity distribution is complex. Taking into account the practical factors of X-ray thermoacoustic imaging, the sound velocity distribution information c(x,y) is embedded in the DAS method, and the 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 sinusoidal graph m (t m ) time delay t m By the detector position r m (x m ,y m ) and reconstruct the target point position r p The Euclidean distance between (x,y) is determined as:

[0078]

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

[0080]

[0081] Where Γ represents the signal from the sound pressure position r p Propagate to the detector position r m The actual propagation path of , dl represents the small arc length element along the signal propagation path Γ.

[0082] S4: training the joint network of the diffusion model and the sound speed compensation module using the training data set;

[0083] S5: Perform domain transformation on the X-ray thermoacoustic signals collected under sparse viewing angles, use the trained network to complete and enhance the images, and obtain high-quality X-ray thermoacoustic image reconstruction results.

[0084] like Figure 4 Specifically, as the number of iterations increases, the noise in the sinusoidal graph of the X-ray thermoacoustic signal gradually decreases, and the reconstructed image gradually becomes clearer. After the 500th iteration, the method of the present invention has essentially restored all the X-ray thermoacoustic signal information, while also obtaining a good reconstructed image result, further improving image quality.

[0085] It is further clarified that the image reconstruction process includes:

[0086] (1) Mapping sparsely sampled thermoacoustic signals into a latent space;

[0087] (2) Image completion and enhancement through trained back-diffusion network;

[0088] (3) Combined with the guidance of sound velocity compensation, the reconstructed image is made to approach the real distribution in the latent space;

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

[0090] In summary:

[0091] The present invention uses a generative model to learn the prior information in the X-ray thermoacoustic signal. By learning the prior distribution of the X-ray thermoacoustic signal sampled under non-sparse conditions, the X-ray thermoacoustic signal under sparse viewing angle is gradually guided and optimized, effectively overcoming the signal sparsity problem in sparse reconstruction. The delayed sum (DAS) method is used to apply domain transformation to the X-ray thermoacoustic signal to map the X-ray thermoacoustic signal to the image space. When using DAS to perform domain transformation, sound velocity compensation is used to further optimize the reconstruction quality. An image reconstruction network based on a diffusion probability model is constructed. Noise is introduced in the diffusion process to simulate signal degradation, and then the true thermoacoustic image is gradually restored through the reverse diffusion process, thereby improving the details and contrast of the reconstructed image. The diffusion model is trained using a training data set containing X-ray thermoacoustic images under full viewing angle and sparse viewing angle to improve its ability to restore images under different sampling conditions. The present invention can achieve high-quality reconstruction and enhancement of X-ray thermoacoustic images under limited acquisition viewing angle or poor signal quality, improve the robustness and imaging effect of the imaging system, and promote the practical application of X-ray thermoacoustic imaging technology in biomedical diagnosis, tissue functional imaging and other fields.

[0092] 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 embodied in other specific forms without departing from the spirit or essential characteristics of the invention. Therefore, the embodiments should be considered in all respects as illustrative and non-restrictive, and the scope of the invention is defined by the appended claims, not the foregoing description, and all variations within the meaning and range of equivalents of the claims are intended to be included therein. Any reference sign in a claim should not be construed as limiting the claim to which it relates.

Claims

1. An X-ray thermoacoustic imaging method based on a diffusion model and sound velocity compensation, characterized by: At least the following steps are included: S1: Generate X-ray thermoacoustic signal simulation data based on MCX and K-Wave joint modeling; S2: Constructing a diffusion model for learning the distribution of X-ray thermoacoustic signals. The diffusion model simulates the distribution evolution of thermoacoustic signals through forward diffusion and reverse diffusion processes. S3: Construct a sound velocity compensation module to incorporate the sound velocity information of the tissue area into the delay-sum reconstruction process to correct the sound wave propagation time and improve the accuracy of the reconstructed image; S4: training the joint network of the diffusion model and the sound speed compensation module using the training data set; S5: Perform domain transformation on the X-ray thermoacoustic signals collected under sparse viewing angles, use the trained network to complete and enhance the images, and 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: Said S1 at least comprises 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 divided; According to the tissue information in the CT image, optical parameters of each tissue are set, wherein the optical parameters include absorption coefficient μa, scattering coefficient μs, anisotropy coefficient g, refractive index n, and thermoacoustic parameters, wherein the thermoacoustic parameters include density ρ, Grünitz coefficient Γ, and speed of sound c; The incident direction, energy spectrum, and photon number of the X-ray light source were set in the tissue model, and the energy deposition distribution H(x, y, z) in three-dimensional space was calculated using Monte Carlo simulation (MCX). Based on the theory of thermoelastic expansion, the energy deposition H(x, y, z) is converted to the initial sound pressure distribution p0(x, y, z). The conversion formula is: p0(x,y,z)=Γ(x,y,z)·H(x,y,z) Where H is the deposition energy per unit volume, Γ(x,y,z) is the Grünigssen coefficient of the tissue; According to the distribution of tissue in the CT image, a corresponding sound velocity value is assigned to each voxel, a non-uniform sound velocity field is constructed, and simulation parameters are configured, including spatial resolution, ultrasonic sensor array arrangement, and sensor parameters; Preset the spatial resolution of the image and define the number, shape, and spatial distribution of the sensor array in K-Wave; set key parameters of the ultrasonic sensor, including 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 thermal acoustic waves in inhomogeneous media, 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; The simulation is repeated under different sampling conditions to generate multiple sets of ultrasonic signal data, and sparse and fully sampled data pairs are constructed to serve as training inputs and targets 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 diffusion model is used to: denoise the received signal and learn 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, is the drift coefficient; is the diffusion coefficient; represents Brownian motion; dw represents Brownian incremental motion; The reverse diffusion process of the diffusion model satisfies: where dt is an infinitesimal time step, is the reverse Brownian motion; The diffusion model is constructed by training the score network S θ (x t ,t) to estimate the score function of time step t The parameter optimization goal is: in, represents the square error between the score network output and the true score function; λ(t) is the weighting factor at time t, which is used to control the influence of different time steps; It is the expectation operation, which means that for different random variables x0 and x t Perform expectation calculation; t～N(0,T) Also for expected operation.

4. The X-ray thermoacoustic imaging method based on diffusion model and sound velocity compensation according to claim 1, characterized in that: The sound velocity compensation method of the compensation module in S3 is: calculate the propagation path Γ of the sound wave in the inhomogeneous medium, calculate the time delay t by path integration m : Where c(x,y) represents the speed of sound at the position (x,y); dl represents the micro-arc length along the propagation path Γ; The sound velocity compensation method uses a finite difference method or a finite element method to simulate sound wave propagation to calculate the sound wave delay compensation parameters in each tissue layer.

5. The X-ray thermoacoustic imaging method based on diffusion model and sound velocity compensation according to claim 1, characterized in that: The S5 at least includes the following steps: Construct a dataset containing original images and corresponding thermoacoustic signals; An end-to-end training method is used, with the difference between the reconstructed image and the original image as the optimization target; The perceptual loss, structural similarity loss and mean square error loss are combined to jointly optimize the score network parameters and sound speed compensation parameters.

Images