A sparse reconstruction method for photoacoustic tomography based on diffusion model

Abstract

The application discloses a photoacoustic tomography sparse reconstruction method based on a diffusion model, and in the photoacoustic tomography reconstruction process, a sparse reconstruction strategy combining a fractional-based diffusion model and a model-based iterative reconstruction method is proposed, a fractional network is used to learn the data distribution of a target image, and the output of the final network is used as prior information of an optimization problem in model iteration to obtain an optimal solution. The method can make the photoacoustic image reconstructed under sparse detection view angles have fewer artifacts than the reconstructed image of the traditional method, and can more effectively and accurately reflect the real information of the target object.

Description

Technical Field

[0001] This invention relates to the field of biomedical imaging technology, specifically to a sparse reconstruction method for photoacoustic tomography based on a diffusion model. Background Technology

[0002] Photoacoustic imaging is a non-invasive imaging technology based on photoacoustic effect that has emerged in the field of biomedical imaging in recent years. It combines the advantages of high contrast in optical imaging and high penetration in acoustic imaging, and has great potential for clinical translation and application prospects.

[0003] Photoacoustic tomography (PAT), an important branch of photoacoustic imaging, can obtain high-resolution images of deep tissues. Currently, PAT has been widely used in tissue cell imaging, cancer detection, cardiovascular disease detection, and image-guided surgery, and is gradually becoming an important tool in preclinical and clinical settings. However, under ill-posed conditions (such as limited or sparse viewpoints), PAT-reconstructed images suffer from low resolution and may contain severe artifacts. Therefore, achieving high-quality PAT reconstruction under sparse sampling remains a pressing problem to be solved. Summary of the Invention

[0004] The purpose of this invention is to address the shortcomings of existing technologies by providing a sparse reconstruction method for photoacoustic tomography based on a diffusion model, in order to solve the problems of low image resolution and severe artifacts in photoacoustic tomography reconstruction under ill-posed conditions (such as finite viewing angle, sparse viewing angle, etc.).

[0005] To achieve the above-mentioned technical objectives, the present invention adopts the following technical solution.

[0006] A sparse reconstruction method for photoacoustic tomography based on a diffusion model includes the following steps:

[0007] Step S1: Based on the photoacoustic signal detected by the ultrasonic transducer, the reconstruction process of the photoacoustic tomography image is transformed into a solution process of a constraint problem based on the least squares method.

[0008] Step S2: Find the optimal solution to the constraint problem. Introduce a prior data distribution output by the diffusion model as a regularization term for constraint. Decouple the prior data term and data fidelity term output by the diffusion model. Use an alternating solution method to alternately update and obtain the optimal solution between the two sub-problems shown in the following formula, thereby obtaining the photoacoustic tomography reconstruction image:

[0009]

[0010] In this equation, the upper part of equation (1) is generated by the diffusion model and serves as the regularization part in the iterative generation process; the lower part of equation (1) is the data fidelity term, which is generated using the gradient descent iterative formula x.k-1 =x k -aA * (Ax k -y) solves the minimum optimization problem.

[0011] Step S1 involves converting the reconstruction process of the photoacoustic tomography image into a solution process for a constraint problem based on the least squares method. The conversion process is as follows:

[0012] In photoacoustic tomography, the object being imaged is irradiated by a short-time pulsed laser. The object absorbs the laser energy and converts some of it into heat. This heat is further converted into a momentarily rising local pressure through thermoelastic expansion and propagates outwards in the form of photoacoustic waves. These photoacoustic waves can be detected by an ultrasonic transducer and converted into photoacoustic signals, which are then used to reconstruct the photoacoustic image. Assuming that the laser excitation simultaneously satisfies both thermal and pressure constraints, the relative changes in absorbed heat diffusion and volume expansion can be neglected. The initial sound pressure generated by the imaged object propagates outwards through the medium. The mathematical model of the propagation equation is expressed as:

[0013]

[0014] In the above formula, ▽ represents the Hamiltonian operator, c represents the speed of sound, β represents the thermoelastic expansion coefficient, and c p It is the specific heat capacity, and p(r,t) is the sound pressure level at position r at time t;

[0015] Furthermore, p(r,t) can be obtained using Green's formula:

[0016]

[0017] The above equation describes the propagation process of photoacoustic imaging, which can be simplified to a linear process as shown in the following equation:

[0018] y = Ax (4)

[0019] In the above formula, y represents the photoacoustic signal p(r,t) detected by the ultrasonic transducer, the linear operator A represents the photoacoustic forward process, and x represents the initial sound pressure p0.

[0020] As shown in equation (3), the reconstruction problem of photoacoustic tomography image is to obtain the initial sound pressure x using y through an optimization algorithm. This optimization problem can be solved by minimizing the least squares error:

[0021]

[0022] In the solution process, a regularization term R(x) is often introduced to constrain the optimization problem to obtain a more accurate solution. Then, equation (4) can be further expressed as:

[0023]

[0024] In the above formula, R(x) is a data consistency term, which is a regularization term that contains prior information about the target image, and λ is a regularization parameter.

[0025] The optimal solution of the objective function shown in equation (5) is obtained by iteratively applying the gradient descent method:

[0026]

[0027] In the above formula, x k This represents the result of the k-th iteration.

[0028] In step S2, a prior data distribution output by the diffusion model is introduced as a regularization term for constraint. The process of the diffusion model outputting the prior data distribution is as follows:

[0029] Let R(x) be the probability density function logp t (x), through training the scoring network s θ (O) Predicting the gradient of the photoacoustic tomography image dataset distribution ▽ x logp t (x), where O is the dataset of the target image, and the preliminary predicted reconstructed image is obtained by numerically solving the inverse SDE based on the score matching generation model:

[0030]

[0031] In the above formula, σ i is the noise scale, i is the iteration number, and z is zero-mean Gaussian white noise; after correcting the initially predicted reconstructed image using the prediction-correction sampler PC, the output prior data distribution is used as a regularization term to constrain the objective function.

[0032] Specifically, the training process of the scoring network uses a denoised score matching method. During training, the parameters θ of the scoring network are optimized using the following formula:

[0033]

[0034] After the network training is complete, use approximate conditions to let

[0035] Specifically, the method for correcting the initially predicted reconstructed image using a prediction-correction sampler (PC) is as follows:

[0036] An Euler discretization method is used to introduce a predictor-corrector sampler (PC) to correct errors in the evolution of the discretization direction SDE, resulting in a corrected reconstructed image.

[0037]

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

[0039] This invention proposes a sparse reconstruction method for photoacoustic tomography based on a diffusion model, which significantly improves reconstruction quality compared to traditional methods under extremely sparse probe angles (64 and 32 probe angles). Under a 64 projection, the PSNR and SSIM of the reconstructed mouse abdominal image reach 27.52 dB and 0.93, respectively, representing improvements of 5.94 dB and 0.25 compared to traditional methods. Under a 32 extremely sparse projection, the PSNR and SSIM reach 22.35 dB and 0.90, respectively, representing improvements of 5.09 dB and 0.32 compared to traditional methods. This effectively solves the problems of low image resolution and severe artifacts inherent in traditional photoacoustic tomography reconstruction methods under ill-posed conditions (such as limited or sparse probe angles). Attached Figure Description

[0040] To more clearly illustrate the technical solutions of the embodiments of this disclosure, the accompanying drawings used in the embodiments will be briefly described below. It should be understood that the following drawings only show some embodiments of this disclosure and should not be regarded as a limitation of the scope. For those skilled in the art, other related drawings can be obtained based on these drawings without creative effort.

[0041] Figure 1 This is a flowchart of the network model algorithm of the present invention;

[0042] Figure 2 This is a schematic diagram of forward and backward propagation in a diffusion model;

[0043] Figure 3 This is a simulation of vascular reconstruction results;

[0044] Figure 4 The result of the circular phantom reconstruction;

[0045] Figure 5 The results show the sparse reconstruction of the mouse abdomen. Detailed Implementation

[0046] To facilitate understanding and implementation of the present invention by those skilled in the art, the various steps of the method proposed in this invention are described in detail below. It should be understood that these embodiments are for illustrative purposes only and are not intended to limit the scope of the invention. Furthermore, it should be understood that after reading the teachings of this invention, those skilled in the art can make various modifications or alterations to the invention, and these equivalent forms also fall within the scope defined by the appended claims.

[0047] Example

[0048] like Figure 1 As shown, this invention provides a sparse reconstruction method for photoacoustic tomography based on a diffusion model, comprising the following steps:

[0049] Step S1: Based on the photoacoustic signal detected by the ultrasonic transducer, the reconstruction process of the photoacoustic tomography image is transformed into a solution process of a constraint problem based on the least squares method.

[0050] Step S2: Find the optimal solution to the constraint problem. Introduce a prior data distribution output by the diffusion model as a regularization term for constraint. Decouple the prior data term and data fidelity term output by the diffusion model. Use an alternating solution method to alternately update and obtain the optimal solution between the two sub-problems shown in the following formula, thereby obtaining the photoacoustic tomography reconstruction image:

[0051]

[0052] In this equation, the upper part of equation (1) is generated by the diffusion model and serves as the regularization part in the iterative generation process; the lower part of equation (1) is the data fidelity term, which is generated using the gradient descent iterative formula x. k-1 =x k -aA * (Ax k -y) solves the minimum optimization problem.

[0053] Step S1 involves converting the reconstruction process of the photoacoustic tomography image into a solution process for a constraint problem based on the least squares method. The conversion process is as follows:

[0054] In photoacoustic tomography, the object being imaged is irradiated by a short-time pulsed laser. The object absorbs the laser energy and converts some of it into heat. This heat is further converted into a momentarily rising local pressure through thermoelastic expansion and propagates outwards in the form of photoacoustic waves. These photoacoustic waves can be detected by an ultrasonic transducer and converted into photoacoustic signals, which are then used to reconstruct the photoacoustic image. Assuming that the laser excitation simultaneously satisfies both thermal and pressure constraints, the relative changes in absorbed heat diffusion and volume expansion can be neglected. The initial sound pressure generated by the imaged object propagates outwards through the medium. The mathematical model of the propagation equation is expressed as:

[0055]

[0056] In the above formula, ▽ represents the Hamiltonian operator, c represents the speed of sound, β represents the thermoelastic expansion coefficient, and c p It is the specific heat capacity, and p(r,t) is the sound pressure level at position r at time t;

[0057] Furthermore, p(r,t) can be obtained using Green's formula:

[0058]

[0059] The above equation describes the propagation process of photoacoustic imaging, which can be simplified to a linear process as shown in the following equation:

[0060] y = Ax (4)

[0061] In the above formula, y represents the photoacoustic signal p(r,t) detected by the ultrasonic transducer, the linear operator A represents the photoacoustic forward process, and x represents the initial sound pressure p0.

[0062] As shown in equation (3), the reconstruction problem of photoacoustic tomography image is to obtain the initial sound pressure x using y through an optimization algorithm. This optimization problem can be solved by minimizing the least squares error:

[0063]

[0064] In the solution process, a regularization term R(x) is often introduced to constrain the optimization problem to obtain a more accurate solution. Then, equation (4) can be further expressed as:

[0065]

[0066] In the above formula, R(x) is a data consistency term, which is a regularization term that contains prior information about the target image, and λ is a regularization parameter.

[0067] The optimal solution of the objective function shown in equation (5) is obtained by iteratively applying the gradient descent method:

[0068]

[0069] In the above formula, x k This represents the result of the k-th iteration.

[0070] In the above process, choosing an appropriate regularization term helps in reconstructing a complete image. Common regularization methods in photoacoustic imaging include Tikhonov regularization and TV regularization, but these conventional regularization methods each have their limitations. This embodiment employs a method combining a diffusion model and gradient descent.

[0071] Specifically, in the solution process of step S2, such as Figure 2 As shown, a prior data distribution output by the diffusion model is introduced as a regularization term for constraint. The process of the diffusion model outputting the prior data distribution is as follows:

[0072] Let R(x) be the probability density function logp t (x), through training the scoring network s θ (O) Predicting the gradient of the photoacoustic tomography image dataset distribution ▽ x logpt (x), where O is the dataset of the target image, and the preliminary predicted reconstructed image is obtained by numerically solving the inverse SDE based on the score matching generation model:

[0073]

[0074] In the above formula, σ i is the noise scale, i is the iteration number, and z is zero-mean Gaussian white noise; after correcting the initially predicted reconstructed image using the prediction-correction sampler PC, the output prior data distribution is used as a regularization term to constrain the objective function.

[0075] Specifically, to improve sampling quality, the scoring network needs to be trained. The training process uses a denoising score matching method, and during training, the parameters θ of the scoring network are optimized using the following formula:

[0076]

[0077] After the network training is complete, use approximate conditions to let

[0078] Specifically, the method for correcting the initially predicted reconstructed image using a prediction-correction sampler (PC) is as follows:

[0079] An Euler discretization method is used to introduce a predictor-corrector sampler (PC) to correct errors in the evolution of the discretization direction SDE, resulting in a corrected reconstructed image.

[0080]

[0081] To verify the effectiveness of this invention, this embodiment simulates the generation of a large number of segmented blood vessels, circular phantoms, and raw photoacoustic signals and initial pressure distributions of the mouse abdomen. The experiments were conducted using a ring-shaped array of sensors, with 32, 64, and 128 sensors respectively, a ring radius of 21.6 mm, a sound velocity of 1500 m / s, and all images were 256×256 pixels in size. The center frequency of the sensors was set to 2.5 MHz. The entire dataset consisted of 1200 training samples and 300 test samples. All experimental procedures were implemented on the deep learning open-source framework PyTorch. The sparse reconstruction effect was observed by comparing the labeled results with the results reconstructed using the method of this invention, and the effectiveness was judged by comparing quantitative indicators.

[0082] The experimental platform was configured with a GPU graphics card (GeForce RTX 2080Ti).

[0083] The data verified in this experiment are simulated blood vessels, circular phantoms, and mouse abdomen reconstructed through 1000 iterations under 32, 64, and 128 probes.

[0084] like Figure 3 The image shows the results of reconstructing simulated blood vessels using different methods under 32, 64, and 128 probes; Figure 4 The image shows the results of reconstructing a circular phantom using different methods under 32, 64, and 128 probes; Figure 5 The figure shows the results of reconstructing the mouse abdomen using different methods under 32, 64, and 128 probes. As can be seen from the results, the method of this invention, compared to DAS and U-net methods, yields the closest reconstruction result to the labeled image using the diffusion model.

[0085] The results of these reconstruction methods were further quantitatively compared by comparing peak signal-to-noise ratio (PSNR) and structural similarity (SSIM).

[0086] The results shown in Tables 1, 3, and 5 below are the SSIM quantitative results of reconstructed images of simulated blood vessels, circular phantoms, and mouse abdomens using traditional reconstruction methods and the method of this invention under 32, 64, and 128 probes. The results shown in Tables 2, 4, and 6 below are the PSNR quantitative results of reconstructed images of simulated blood vessels, circular phantoms, and mouse abdomens using traditional reconstruction methods and the method of this invention under 32, 64, and 128 probes.

[0087] Table 1 Figure 3 PSNR values ​​of each image

[0088] DAS 7.62 8.77 9.06 U-Net 24.39 32.37 36.08 Diffusion Model 38.50 41.10 45.09

[0089] Table 2 Figure 3 SSIM values ​​of each image

[0090] DAS 0.15 0.19 0.21 U-Net 0.88 0.94 0.95 Diffusion Model 0.97 0.98 0.99

[0091] Table 3 Figure 4 PSNR values ​​of each image

[0092] DAS 17.26 21.58 23.90 U-Net 21.19 26.07 27.33 Diffusion Model 22.35 27.52 31.34

[0093] Table 4 Figure 4 SSIM values ​​of each image

[0094] DAS 0.25 0.41 0.53 U-Net 0.69 0.73 0.78 Diffusion Model 0.90 0.93 0.95

[0095] Table 5 Figure 5 PSNR values ​​of each image

[0096]

[0097]

[0098] Table 6 Figure 5 PSNR values ​​of each image

[0099] DAS 0.53 0.73 0.80 U-Net 0.91 0.93 0.89 Diffusion Model 0.97 0.98 0.98

[0100] The data in the table above clearly demonstrate the superior quantitative performance of the method of this invention compared to traditional reconstruction methods.

[0101] It should be understood that the above description of the preferred embodiments is quite detailed, but it should not be considered as a limitation on the scope of protection of this invention. Those skilled in the art, under the guidance of this invention, can make substitutions or modifications without departing from the scope of protection of the claims of this invention, and all such substitutions or modifications fall within the scope of protection of this invention. The scope of protection of this invention should be determined by the appended claims.

Claims

1. A sparse reconstruction method for photoacoustic tomography based on a diffusion model, characterized in that, Includes the following steps: Step S1: Based on the photoacoustic signal detected by the ultrasonic transducer, the reconstruction process of the photoacoustic tomography image is transformed into a solution process of a constraint problem based on the least squares method. Step S2: Find the optimal solution to the constraint problem. Introduce a prior data distribution output by the diffusion model as a regularization term for constraint. Optimize and decouple the prior data term and data fidelity term output by the diffusion model. Use an alternating solution method to alternately update and obtain the optimal solution between the two sub-problems shown in the following equation, thereby obtaining the photoacoustic tomography reconstruction image: （1）； In this equation, the first formula is generated by the diffusion model and serves as the regularization part in the iterative generation process; the second formula in equation (1) is the data fidelity term, which is generated using the gradient descent iterative formula. Solve the minimum optimization problem; The process by which the diffusion model outputs the prior data distribution is as follows: make The probability density function By training the scoring network Predicting the gradient of the photoacoustic tomography image dataset data distribution ,in, This is a dataset of target images. Based on a score-matching generation model, preliminary predicted reconstructed images are obtained by numerically solving the inverse SDE. （8）； In the above formula, It is a noise scale. It is the number of iterations. It is zero-mean Gaussian white noise; after correcting the initially predicted reconstructed image using a prediction-correction sampler (PC), the output prior data distribution is used as a regularization term to constrain the objective function.

2. The sparse reconstruction method for photoacoustic tomography based on a diffusion model according to claim 1, characterized in that, Step S1 involves converting the reconstruction process of the photoacoustic tomography image into a solution process for a constraint problem based on the least squares method. The conversion process is as follows: In photoacoustic tomography, the object being imaged is irradiated by a short-time pulsed laser. The object absorbs the laser energy and converts some of it into heat. This heat is further converted into a momentarily rising local pressure through thermoelastic expansion and propagates outwards in the form of photoacoustic waves. These photoacoustic waves can be detected by an ultrasonic transducer and converted into photoacoustic signals, which are then used to reconstruct the photoacoustic image. Assuming that the laser excitation simultaneously satisfies both thermal and pressure constraints, the relative changes in absorbed heat diffusion and volume expansion can be neglected. The initial sound pressure generated by the imaged object propagates outwards through the medium. The mathematical model of the propagation equation is expressed as: （2）； In the above formula, Represents the Hamiltonian operator. Indicates the speed of sound. Indicates the coefficient of thermoelastic expansion. It is specific heat capacity. Is Time and location The magnitude of the sound pressure level; Furthermore, We can continue to use Green's theorem to solve for the following: （3）； The above equation describes the propagation process of photoacoustic imaging, which can be simplified to a linear process as shown in the following equation: （4）； In the above formula, Represents the photoacoustic signal detected by the ultrasonic transducer Linear operators Represents the positive process of photoacoustic transmission. Represents the initial sound pressure level ; As can be seen from equation (3), the problem of reconstructing photoacoustic tomography images is to utilize... The initial sound pressure was obtained through an optimization algorithm. This optimization problem can be solved by minimizing the least squares error: （5）； In the solution process, a regularization term is often introduced. By imposing constraints on the optimization problem to obtain a more accurate solution, equation (4) can be further expressed as: （6）； In the above formula, It is a data consistency item. It is a regularization term that includes prior information about the target image. It is a regularization parameter; The optimal solution of the objective function shown in equation (5) is obtained by iteratively applying gradient descent: （7）； 3. The sparse reconstruction method for photoacoustic tomography based on a diffusion model according to claim 1, characterized in that, The training process of the scoring network uses a denoised score matching method. During training, the parameters of the scoring network are adjusted. Optimize using the following formula: （9）； After the network training is complete, use approximate conditions to let .

4. The sparse reconstruction method for photoacoustic tomography based on a diffusion model according to claim 1, characterized in that, The method for correcting the initially predicted reconstructed image using a prediction-correction sampler (PC) is as follows: An Euler discretization method is used to introduce a predictor-corrector sampler (PC) to correct errors in the evolution of the discretization direction SDE, resulting in a corrected reconstructed image. （10）。

Images