Three-dimensional ultrasonic tomography sound velocity imaging method based on multi-template fast marching method

Abstract

The invention relates to the technical field of ultrasonic tomography, and discloses a three-dimensional ultrasonic tomography sound velocity imaging method based on a multi-template fast marching method, which comprises the following steps: acquiring parameters of a three-dimensional ultrasonic tomography system and initial sound velocity distribution of biological tissues, and establishing a curved ray model based on an eikonal equation; discretizing the propagation time gradient in the eikonal equation by adopting second-order finite difference approximation, and covering 18 neighborhoods of the grid by adopting templates in four directions to carry out combined solution, so as to obtain ultrasonic propagation time distribution of the whole grid; and iteratively reconstructing a three-dimensional sound velocity distribution image through a maximum posterior probability estimation algorithm based on the ultrasonic propagation time distribution. According to the method, the propagation time gradient in the second-order finite difference approximation eikonal equation is adopted to reduce discretization errors, a plurality of direction templates are used, more accurate difference approximation is carried out on gradient operators, and the non-orthogonal propagation path of the ultrasonic waves in the three-dimensional space is completely captured.

Description

Technical Field

[0001] This invention relates to the field of ultrasonic tomography technology, and specifically to a three-dimensional ultrasonic tomography sound velocity imaging method based on a multi-template rapid travel method. Background Technology

[0002] Ultrasound computed tomography (USCT) is a novel non-invasive medical imaging technology that has emerged in recent years. With its core advantages of being radiation-free, having a short examination time, and offering diverse imaging modalities, it has become an important research direction in clinical scenarios such as breast tissue imaging and soft tissue lesion detection. Its technical principle is as follows: A transducer array with a specific spatial geometry (classified by structure as two-dimensional arrays such as linear / ring arrays, and three-dimensional arrays such as hemispherical / ellipsoidal arrays) controls each element to emit ultrasound waves sequentially at a specific frequency. When ultrasound waves propagate in biological tissue, they are refracted and transmitted due to differences in the acoustic properties (such as sound velocity and acoustic impedance) of different tissues, and the transmitted signal waveform is ultimately received by other array elements. Based on the signals acquired across the entire range, the spatial distribution of the acoustic properties of the biological tissue (sound velocity image) is inverted through image reconstruction algorithms, enabling visualization of tissue structure and lesions.

[0003] In the ultrasound tomography technology system, data acquisition and image reconstruction are the two core links. The quality of sound velocity image reconstruction directly depends on the accuracy of the forward modeling of sound wave propagation—that is, whether the propagation process of ultrasound waves in biological tissues can be accurately simulated. This link has become a key technical bottleneck that determines the clinical application value of ultrasound tomography.

[0004] Forward modeling of ultrasound tomography aims to describe the propagation of ultrasound waves in biological tissues. The core is establishing a mapping relationship between "tissue acoustic properties and ultrasound propagation time." Currently, the mainstream technical approaches fall into three categories, each with significant differences and limitations: (1) Full-wave inversion (FWI) model, which has high accuracy but low efficiency: The full-wave inversion model is built upon the complete acoustic wave equations and can accurately simulate finite-frequency phenomena such as diffraction, interference, and dispersion of ultrasound waves, making it the most accurate forward model currently available. However, it requires solving complex three-dimensional wave equations, resulting in an exponential increase in computational load and extremely low imaging efficiency (e.g., reconstruction of a single ultrasound tomography imaging session can take several hours to several days). This makes it difficult to meet the clinical demand for rapid imaging and limits its application to small-scale, high-precision research scenarios, hindering its large-scale deployment.

[0005] (2) Straight-ray model: High efficiency but insufficient accuracy. To balance efficiency and accuracy, early studies proposed a linear model, which assumes that ultrasound waves propagate in a straight line within biological tissues, ignoring refraction effects at interfaces between media with different acoustic properties. While this model significantly simplifies calculations and improves imaging speed, it violates the true physical laws of propagation—significant refraction inevitably occurs at interfaces between different structures in biological tissues (such as glands and fat in the breast, or normal tissue and lesions in the liver) leading to a large deviation between the propagation time in the forward modeling and the actual value. Reconstruction results based on this model are prone to problems such as blurred boundaries and lesion localization errors, failing to meet the clinical needs for detecting subtle lesions.

[0006] (3) Bent-ray model, the mainstream choice that balances efficiency and accuracy: To resolve the aforementioned contradictions, researchers approximate the wave equation at high frequencies into an Eikonalequation, introducing geometric acoustics theory to model ultrasonic wave propagation as a "refracted ray," thus proposing a curved ray model. The core advantage of this model lies in its consideration of refraction effects at interfaces between media with different acoustic properties, while neglecting finite-frequency phenomena such as diffraction and interference that significantly impact computational complexity, achieving a balance between accuracy and efficiency. Compared to full-wave inversion, its numerical solution efficiency is improved by 1-2 orders of magnitude; compared to linear models, its reconstruction accuracy is significantly improved, making it the mainstream technical approach for forward modeling in current ultrasonic tomography, and also the acoustic wave propagation forward model used in this invention.

[0007] Currently, the mainstream numerical method for solving the equation is the Fast Marching Method (FMM), which is based on Dijkstra's shortest path algorithm. It iteratively propagates outward from the source grid, updating the propagation time of adjacent grids according to the rule of "narrow-band concentrated minimum propagation time priority", and finally obtains the propagation time distribution of the entire space.

[0008] The Fast Marching Method (FMM) is widely used in 2D ultrasound tomography reconstruction due to its high computational efficiency and acceptable accuracy in solving 2D scenes. However, when extended to 3D ultrasound tomography (3D-USCT), it reveals a bottleneck in solution accuracy, with the following specific drawbacks: (1) The solution accuracy of the functional equation is insufficient, and the discretization error is significant: The traditional Fast Marching Method (FMM) uses a first-order finite difference approximation to discretize the partial derivatives in the equation. It has low sensitivity to regions in biological tissues where the sound velocity changes rapidly (such as the tumor-normal tissue interface), and cannot accurately capture subtle changes in propagation time. This leads to the accumulation of discretization errors, which directly affects the fidelity of the sound velocity reconstruction.

[0009] (2) Incomplete directional information coverage and missing simulation of non-orthogonal propagation paths: The traditional Fast Marching Method (FMM) solves the discrete-time function equations using only one orthogonal axis template, which covers only the "6-neighborhood" of the mesh (along...). Adjacent grids in orthogonal directions are ignored, while grid information in non-orthogonal (diagonal) directions is disregarded. In 3D ultrasound tomography, ultrasound waves propagate diagonally in many complex biological tissues (e.g., obliquely entering glandular tissue from adipose tissue). Traditional Fast Traveling Method (FMM) lacks simulation of these propagation paths, leading to large deviations in propagation time calculations and causing "strip artifacts" and "blurred boundaries" in the reconstructed image. The template refers to the three directions selected when performing a difference approximation of the gradient operator in a 3D Cartesian coordinate system.

[0010] (3) Iterative reconstruction has poor stability and errors are prone to accumulation and divergence: Due to the aforementioned discretization errors and lack of directional information, the traditional Fast Marching Method (FMM) will continuously accumulate errors during the iterative reconstruction of three-dimensional ultrasound tomography, making it impossible to maintain stable high precision and meet the clinical requirements for the reliability and repeatability of reconstruction results.

[0011] (4) The three-dimensional ultrasound tomography imaging has poor scene adaptability and cannot meet the requirements of high-fidelity imaging: The traditional Fast Marching Method (FMM) was originally designed for two-dimensional scenes. Its technical architecture did not take into account the complexity of ultrasound propagation in three-dimensional space. The aforementioned defects in accuracy, directional coverage, and stability resulted in low image quality (numerous artifacts, high noise, and blurred boundaries) of three-dimensional ultrasound tomography sound velocity reconstruction based on this method. It could not meet the needs of accurate identification of subtle lesions in clinical scenarios such as breast tumor screening and soft tissue lesion detection. Summary of the Invention

[0012] To address the aforementioned shortcomings of the traditional Fast Movement Method (FMM) in three-dimensional ultrasound tomography sound velocity reconstruction, this invention provides a three-dimensional ultrasound tomography sound velocity imaging method based on the Multi-Template Fast Movement Method (MSFM). By adapting to the Multi-Template Fast Movement Method (MSFM) for three-dimensional ultrasound tomography, the following objectives are achieved: (1) Improve the accuracy of solving the equation: enhance the sensitivity to changes in the sound velocity of biological tissues, reduce discretization errors, accurately capture the propagation law of ultrasound in areas where the sound velocity changes rapidly, and lay the foundation for high-quality sound velocity reconstruction.

[0013] (2) Complete the coverage of directional information and accurately simulate non-orthogonal propagation paths: avoid reconstruction artifacts caused by missing directional information and clearly restore tissue interfaces (such as tumor boundaries).

[0014] (3) Ensure the stability of iterative reconstruction and avoid error accumulation and divergence: Through the multi-template joint solution strategy, error compensation is achieved, the error accumulation in the iteration process is suppressed, and the three-dimensional ultrasound tomography reconstruction can still maintain stable high accuracy in multiple iterations.

[0015] (4) Optimize the adaptability of three-dimensional ultrasound tomography to meet the clinical high-fidelity imaging requirements: Through the above optimization of accuracy, directional coverage and stability, the multi-template rapid travel method (MSFM) can be adapted to the complex imaging scenarios of three-dimensional ultrasound tomography, and output a sound velocity distribution map with low artifacts, low noise and high boundary fidelity, providing technical support for soft tissue imaging and breast tumor screening, and promoting three-dimensional ultrasound tomography from scientific research to clinical application.

[0016] To solve the above-mentioned technical problems, the present invention adopts the following technical solution: A three-dimensional ultrasonic tomographic velocity imaging method based on a multi-template fast travel method includes: The array parameters of the three-dimensional ultrasound tomography system and the sound velocity distribution of the medium in which the biological tissue is located are obtained. The computational domain is discretized into a grid, and a curved ray model based on the equation of process function is established. The multi-template fast travel method is used to solve the equation: the propagation time gradient in the equation is discretized by the second-order finite difference approximation, and the 18 neighborhoods of the grid are covered by four directional templates for joint solution to obtain the ultrasonic propagation time distribution of the whole grid. Based on the ultrasonic wave propagation time distribution, a three-dimensional sound velocity distribution image is reconstructed iteratively using the maximum a posteriori probability estimation algorithm.

[0017] In one embodiment, the functional equation is: ; Indicates coordinates as The propagation time of the grid gradient, Coordinates are The sound velocity value of the grid.

[0018] In one embodiment, solving the equation using the multi-template fast travel method further includes a process of classifying the mesh based on whether the propagation time has been determined: The grid within the computational domain is divided into a frozen set, a narrowband set, and an unknown set; The propagation time of the grid within the frozen set has been determined and will not be updated. It serves as the benchmark for solving the propagation time from the sound source to other grids. Initially, only the grid containing the sound source belongs to the frozen set. The grid within the narrowband set is the boundary grid between the frozen set and the unknown set. The propagation time from the sound source point to the current grid has been calculated using adjacent grids belonging to the frozen set. The propagation time from the sound source point to the current grid has not yet been calculated for the grid in the unknown set. After the grid in the narrowband set is processed, the grid in the unknown set will be gradually included in the solution range of the equation. In solving the equations, the grid with the shortest propagation time is always selected from the narrowband set for priority processing.

[0019] In one embodiment, the discretization of the propagation time gradient in the equation using a second-order finite difference approximation specifically includes: ; Indicates the time of transmission. Indicates coordinates as The propagation time of the grid, Indicates coordinates as The propagation time of the grid, Indicates coordinates as The propagation time of the grid, Indicates coordinates as The propagation time of the grid, Indicates coordinates as The propagation time of the grid, Indicates coordinates as The propagation time of the grid, Indicates coordinates as The propagation time of the grid, These are the x-axis, y-axis, and z-axis coordinates of the grid, respectively. These represent the spatial grid step size along the x-axis, y-axis, and z-axis, respectively, which is the distance between two adjacent discrete grid points on the orthogonal axes.

[0020] In one embodiment, the 18 neighborhoods of the grid include 6 orthogonal neighborhoods and 12 diagonal neighborhoods of the grid.

[0021] In one embodiment, the method of using four-directional templates to cover 18 neighborhoods of the mesh for joint solution to obtain the ultrasonic wave propagation time distribution of the entire mesh specifically includes: Adjust the finite difference approximation of the derivative of the fast travel method according to a specific direction template: (9) in, This represents the approximate step size of the difference in the current direction template. Indicates the propagation time of the current grid. This represents the propagation time sign substitution formula for adjacent grids in the three differential directions corresponding to the current direction template. This represents the three difference directions corresponding to the current direction template. This represents the coordinate vector of the current grid. This represents the coordinate vectors of the adjacent grids in the three differential directions corresponding to the current direction template. This represents the current sound speed value in the grid. Transform equation (9) into a quadratic equation and calculate the numerical solution; Numerical solutions are calculated for each direction template, and the minimum value among these numerical solutions is selected as the final solution of the equation for the current grid. .

[0022] In one embodiment, the quadratic equation is in the form of: ; in, These are the coefficients corresponding to the current direction template when solving in the difference direction u. The spatial grid step size, coefficient ,coefficient ,coefficient , The specific form is determined by the difference approximation order and the direction template.

[0023] In one embodiment, the step of iteratively reconstructing the three-dimensional sound velocity distribution image based on the ultrasonic wave propagation time distribution using a maximum a posteriori probability estimation algorithm specifically includes: The imaging problem is transformed into an optimization problem that maximizes the posterior probability. The ultrasonic propagation time distribution and the prior probability of the estimated independent variable are combined using Bayes' theorem to iteratively solve the optimization problem. The sound velocity distribution is updated using the conjugate gradient method, thereby suppressing measurement noise and stabilizing the reconstruction process using prior constraints, ultimately obtaining a three-dimensional sound velocity distribution image.

[0024] In one embodiment, transforming the imaging problem into an optimization problem of maximizing the posterior probability specifically includes: In the curved ray model, Indicates the first The slowness value of each grid cell, where slowness is the reciprocal of the speed of sound, and the grid index is... Arranged linearly in rows; Indicates the index of the ray; Indicates the first The ray passed through the first The intercept size of each grid; Indicates the first The propagation time of a ray; the sound wave propagation process in ultrasound tomography is represented as: (12) Transform equation (12) into matrix form: ; This is the ray path matrix, used to store the intercept information of the rays on each grid. Let be the time vector of ray propagation. This is the distribution vector of slowness values. The total number of grid cells. The total number of rays; The problem with three-dimensional sound velocity distribution images is that they are observed... In the case of iterative solution Minimize the value of the objective function: ; Let be the objective function. It is the F2 norm.

[0025] In one embodiment, the step of combining the ultrasonic wave propagation time distribution with the prior probability of the estimated independent variable using Bayes' theorem to iteratively solve the optimization problem, and updating the sound velocity distribution using the conjugate gradient method, thereby suppressing measurement noise while stabilizing the reconstruction process using prior constraints, and finally obtaining a three-dimensional sound velocity distribution image, specifically includes: set up For the first The unknown noise during the ray measurement process follows a Gaussian distribution. And they are independent of each other. For the first The vector formed by the intercepts of the rays on all grid lines. for The standard deviation of the Gaussian distribution it follows If it is a Gaussian distribution, then probability density function for: ; Similarly, conditional probability density function for: ; Then the joint probability density function for: ; in, noise vector The covariance matrix, Represents a determinant. Represents the identity matrix. Represents matrix transpose and noise vector ; Let the slowness value distribution vector be... Follows Gaussian distribution ,in, Indicates the first The distribution vector of the slowness value in the next iteration, then Prior probability density function for: ; in, The distribution vector of slowness values The covariance matrix, for The standard deviation of the Gaussian distribution it follows; According to Bayes' theorem, solving the problem of finding a three-dimensional sound velocity distribution image is equivalent to maximizing... posterior probability : ; After derivation, we obtain Approximate iterative formula: ; The current ray path matrix is ​​obtained by solving the conjugate gradient descent method. Below, the distribution vector of the slowness value in the (k+1)th iteration. ; After iteratively solving for the slowness distribution vector that satisfies the set convergence threshold, taking the reciprocal yields the current ray path matrix. A three-dimensional sound velocity distribution image.

[0026] For the ray path matrix With a fixed value, the slowness value The iterative solution process is called the inner loop. When the solution result of the inner loop makes the objective function reach the preset convergence threshold... At this point, the inner loop ends, the current sound velocity reconstruction result is input into the multi-template fast traversal method to obtain the updated propagation time distribution, and the ray path matrix is ​​updated accordingly. Update the ray path matrix. The process of iteratively solving the problem within the inner loop is called the outer loop. When the solution obtained from the outer loop causes the objective function to reach a preset convergence threshold... At this point, the outer loop ends, and the final three-dimensional sound velocity distribution image is obtained.

[0027] Compared with the prior art, the beneficial technical effects of the present invention are: This invention addresses the shortcomings of the traditional Fast Marching Method (FMM) for solving the equations in 3D ultrasound tomography, with its core innovation being "second-order finite difference + multi-template joint solution." This invention reduces discretization errors by approximating the propagation time gradient in the equations using second-order finite difference; and by using four directional templates to cover 18 neighborhoods of the grid (6 orthogonal neighborhoods + 12 diagonal neighborhoods), it performs a more accurate difference approximation of the gradient operator, fully capturing the non-orthogonal propagation path of ultrasound waves in 3D space. This invention overcomes the shortcomings of existing methods that only use first-order finite difference and single-directional template information to solve the equations, leading to inaccuracies in the ultrasound propagation time distribution. The low accuracy of previous methods affected the quality of subsequent 3D ultrasound tomography sound velocity image reconstruction. The 3D sound velocity reconstruction method in this invention offers better image quality (reduced artifacts and noise, improved boundary fidelity) and higher stability in the iterative reconstruction process. Attached Figure Description

[0028] Figure 1 This is a schematic diagram of the ultrasonic transducer array acquiring signals in an embodiment of the present invention, and a curved ray model taking a two-dimensional scene as an example.

[0029] Figure 2 This is a schematic diagram of the three-dimensional spatial distribution of 1024 transducer array elements in an embodiment of the present invention.

[0030] Figure 3 This is a schematic diagram of the calculation process of the fast traversal algorithm.

[0031] Figure 4 This is a schematic diagram of a three-dimensional Cartesian domain solution template in an embodiment of the present invention.

[0032] Figure 5 This is a schematic diagram of a three-dimensional electronic phantom in an embodiment of the present invention.

[0033] Figure 6 This is a comparison of the true velocity of sound (GT) image, fast moving method (FMM) reconstruction results, and multi-template fast moving method (MSFM) reconstruction results of the Shepp-Logan phantom in the embodiments of the present invention.

[0034] Figure 7 This is a comparison of the true velocity of sound (GT) image, fast-moving method (FMM) reconstruction results, and multi-template fast-moving method (MSFM) reconstruction results of the bronchial phantom in the embodiments of the present invention.

[0035] Figure 8 This is a curve showing the sound velocity values ​​distributed along the dashed line in the reconstructed image of this invention.

[0036] Figure 9 This is a flowchart of the method in an embodiment of the present invention. Detailed Implementation

[0037] A preferred embodiment of the present invention will now be described in detail with reference to the accompanying drawings.

[0038] like Figure 9 As shown, a three-dimensional ultrasonic tomographic velocity imaging method based on a multi-template rapid travel method according to the present invention includes the following steps: S1: Obtain the array parameters of the three-dimensional ultrasound tomography system and the sound velocity distribution of the medium in which the biological tissue is located; discretize the computational domain into a grid; and establish a curved ray model based on the process function equation. S2, use the multi-template fast travel method to solve the equation: the propagation time gradient in the equation is discretized by the second-order finite difference approximation, and the 18 neighborhoods of the grid are covered by four directional templates for joint solution to obtain the ultrasonic propagation time distribution of the whole grid. S3, Based on the ultrasonic wave propagation time distribution, the three-dimensional sound velocity distribution image is reconstructed iteratively using the maximum a posteriori probability estimation algorithm.

[0039] The technical problem this invention aims to solve is the insufficient accuracy of the traditional Fast Marching Method (FMM) in solving the equations in three-dimensional ultrasound computed tomography (3D-USCT) scenarios. This results in artifacts and noise in the reconstructed sound velocity images, and the iterative process is prone to divergence, making it difficult to meet the requirements for high-fidelity simulation of ultrasound propagation in complex acoustic media and high-quality reconstruction of three-dimensional sound velocity images. This invention improves the accuracy of solving the equations in three-dimensional scenarios based on the Multi-stencil Fast Marching Method (MSFM), thereby improving the image quality of reconstructed sound velocity images from three-dimensional ultrasound computed tomography and ensuring iterative stability.

[0040] I. Overall framework of the technical solution.

[0041] The core of this invention is "improving the accuracy of forward modeling of ultrasonic propagation." It optimizes the solution process of the equation using the Multi-Template Fast Traveling (MSFM) method, providing high-fidelity propagation time data for subsequent sound velocity reconstruction. The overall process is divided into four stages: "data input and model establishment - equation solving - 3D sound velocity image reconstruction - result output and quality verification." Each stage works collaboratively to generate high-quality sound velocity images for 3D ultrasonic tomography. The framework logic is as follows: (1) Data Input and Model Establishment: Obtain the parameters of the three-dimensional ultrasound tomography system (transducer array, computational domain grid) and the initial sound velocity distribution of biological tissue, and establish a curved ray model (forward model), such as Figure 1 As shown. Figure 1 (a) in the figure is a schematic diagram of the ultrasonic transducer array acquiring signals in an embodiment of the present invention; Figure 1 (b) is a schematic diagram of the curved ray model in a two-dimensional scene in an embodiment of the present invention.

[0042] (2) Solving the equation: Based on the "second-order finite difference + multi-template joint solution" strategy of the multi-template fast travel method (MSFM), the ultrasonic propagation time distribution of the whole grid is calculated. This addresses the issue of insufficient accuracy when using the traditional Fast Marching Method (FMM) to solve equations in 3D scenes.

[0043] (3) Three-dimensional sound velocity image reconstruction: Ultrasonic wave propagation time distribution output by the multi-template fast travel method (MSFM). Based on the forward modeling of the curved ray model, the sound velocity distribution image in the 3D scene is reconstructed iteratively using the Maximum a posteriori (MAP) estimation algorithm. .

[0044] (4) Results output and quality verification: A three-dimensional sound velocity distribution map is generated, and the improvement effect of reconstruction quality compared with the traditional reconstruction method (FMM) is verified by quantitative evaluation indicators such as root mean square error (RMSE) and peak signal-to-noise ratio (PSNR).

[0045] 2. Simulation parameters for 2D and 3D ultrasound tomography.

[0046] The simulation parameters are set as follows: the transducer array consists of 1024 elements, uniformly distributed on a sphere with a radius of 130 mm (e.g., ...). Figure 2 (As shown in the figure); the computational domain is a cube of 300mm×300mm×300mm, which is discretized into 150×150×150 grids, with each grid being 2mm in size; the forward propagation model adopts a curved ray model based on the process function equation, which is solved by the traditional fast marching method (FMM) and the multi-template fast marching method (MSFM), respectively, both of which use the second-order finite difference approximation; the reconstruction process uses the maximum a posteriori probability (MAP) estimation algorithm, which is iterated 4 times to optimize the reconstruction results.

[0047] III. Technical details of solving equations using the Multi-Template Fast Marching Method (MSFM).

[0048] The Multi-Template Fast Marching Method (MSFM) is the core innovation of this invention. It achieves high-precision solution of the equation (formula (1)) in a three-dimensional scene through "mesh classification - second-order difference discretization - multi-template joint solution", providing a foundation for subsequent three-dimensional sound velocity reconstruction.

[0049] (1) The specific technical details are as follows.

[0050] 3.1 Grid classification mechanism: Based on the traditional Fast Marching Method (FMM) mesh management logic, the mesh in the 3D computational domain is divided into three categories to ensure an orderly and efficient solution process: Frozen set: Propagation time has been determined. The grid, the points of this type Fixed and no longer updated, serving as the benchmark for solving the propagation time from the sound source point to other grids; initially, only the sound source point (transducer element position) belongs to the frozen set, and its propagation time... .

[0051] Narrow set: The boundary grid between the frozen set and the unknown set, which has been preliminarily calculated using adjacent frozen points. However, it requires subsequent verification and updates and is the core transition region for iterative solution.

[0052] Unknown set: Not yet calculated The grid will be gradually incorporated into the solution range after the narrow-band concentrated points have been processed.

[0053] Throughout the solution process, narrowband sets are always selected. Minimum grid size priority processing (analogous to the "shortest path first" logic of Dijkstra's algorithm) ensures the causality of ultrasonic wave propagation (it reaches the near-field grid first, then propagates to the far-field), such as... Figure 3 As shown (for ease of understanding, the solution process is visualized using a two-dimensional scene).

[0054] 3.2 Second-order finite difference discretization (improving discretization accuracy): Traditional Fast Marching Method (FMM) addresses the propagation time gradient in the equation function. The first-order difference approximation scheme is adopted, as shown in formula (2): (2) To improve the propagation time gradient in the three-dimensional Cartesian coordinate domain To improve the accuracy of the difference approximation and reduce discretization errors, a second-order finite difference approximation is used to enhance sensitivity to regions of rapid sound velocity changes in biological tissues. The specific discretization formula is as follows: (3) 3.3, Multi-template joint solution: To overcome the limitation of the traditional Fast Marching Method (FMM) in 3D scenes where "a single template only covers 6 neighborhoods," this invention designs templates with four directions to cover 18 neighborhoods of the mesh (6 orthogonal neighborhoods + 12 diagonal neighborhoods). A more accurate differential approximation is applied to the gradient operator, resulting in a more complete capture of the non-orthogonal propagation paths of ultrasonic waves in 3D space. The template structure and coverage area are as follows: Figure 4 As shown, the traditional Fast Movement Method (FMM) only considers template a, while the Multi-Template Fast Movement Method (MSFM) comprehensively considers the direction template ad. Figure 4 middle, , , These represent the three difference directions of the template in that direction. Figure 4 In this context, (a) is the template for the traditional rapid movement method, denoted as direction template a. Figure 4 In the above, (b), (c), and (d) are templates introduced by the multi-template rapid movement method, and are denoted as direction template b, direction template c, and direction template d, respectively.

[0055] The specific technical differences are as follows: (1) Traditional Fast Marching Method (FMM) The propagation time gradient in the equation (Equation (1)) is approximated using the finite difference method. Discretize the differential operator by replacing it with a discrete form (taking the first-order difference approximation as an example): (4) To simplify formula (4), the following symbol substitutions are introduced: (5) Substitute formula (5) into formula (4), and assume Formula (4) can be simplified to formula (6): (6) Observation shows that formula (6) can be transformed into... A quadratic equation with as the independent variable: (7) in, , The specific form is determined by the order of the difference approximation (see direction template a in Table 1).

[0056] Table 1. Coefficients corresponding to different numerical solution schemes of the Multi-Template Fast Marching Method (MSFM):

[0057] Solving this quadratic equation yields two numerical solutions to the functional equation. Selecting the larger solution satisfies the causality requirement. The result obtained using the traditional Fast Marching Method (FMM) is then obtained. The solution result for a single-direction template (direction template a) : (8) (2) Multi-template Fast Marching Method (MSFM) The Multi-Template Fast Marching (MSFM) method solves the equations using multiple templates in different directions. This method employs four templates to cover 18 neighborhoods of the central grid (e.g., ...). Figure 4 The use of a medium template (ad) significantly improves the numerical solution accuracy in non-orthogonal directions (diagonal directions), thus enabling the capture of more complex sound wave propagation paths. Correspondingly, the solution algorithm also needs to be improved to adapt to this design.

[0058] First, the finite difference approximation of the derivative needs to be adjusted according to a specific directional template. Therefore, equation (6) is replaced with: (9) in, This represents the approximate step size of the differential template in the current direction.

[0059] Similarly, formula (9) can also be transformed into a quadratic equation similar to that shown in formula (7) for numerical solution: (10) in, The specific form is determined by the difference approximation order and the direction template (see the direction template ad in Table 1).

[0060] Finally, for each template, numerical solutions that satisfy the causality requirement are calculated, and the minimum value among these solutions is selected as the final solution of the grid equation. : (11) These are the numerical solutions calculated for direction templates a, b, c, and d, respectively.

[0061] IV. MAP-based 3D sound velocity image reconstruction algorithm.

[0062] Output using Multi-Template Fast Marching (MSFM) Based on the forward modeling of the curved ray model, the sound velocity distribution image in the 3D scene is iteratively reconstructed using the maximum a posteriori probability (MAP) estimation algorithm. .

[0063] The principle of the curved ray model is as follows: Figure 1As shown in (b) (taking a two-dimensional scene as an example), Indicates the first The slowness value (the reciprocal of the speed of sound) of each grid, with the grid number index arranged linearly by row; Indicates the sequence number of the ray; Indicates the first The ray passed through the first The intercept size of each grid; Indicates the first The propagation time of the rays. Therefore, the propagation process of the sound waves in ultrasound tomography can be represented by equation (12), or it can be transformed into the matrix form of equation (13).

[0064] (12) (13) Rewriting equation (13) in least squares form, the problem of three-dimensional sound velocity image reconstruction is essentially about observing... In the case of iterative solution Minimize the value of the objective function: (14) This invention uses the maximum a posteriori probability estimation (MAP) method based on Bayes' theorem to iteratively solve the least squares problem in order to reconstruct the optimal three-dimensional sound speed distribution image.

[0065] The main idea of ​​the maximum a posteriori probability estimation algorithm is to use the prior probability of the pre-estimated independent variable and several experimental results to infer and maximize the posterior probability of the pre-estimated independent variable according to Bayes' theorem, and finally use this value as the true value of the estimated independent variable. The pre-estimated independent variable is the finally reconstructed three-dimensional slowness distribution vector.

[0066] Assumption For the first The unknown noise during the ray measurement process follows a Gaussian distribution. And they are independent of each other, then The probability density function is: (15) Similarly, in the slowness distribution... Arrival time The conditional probability density function is: (16) Then the joint probability density function for: (17) in, noise vector The covariance matrix.

[0067] Assume the distribution vector of slowness values Follows Gaussian distribution ,in Let represent the distribution vector of the slowness values ​​in the k-th iteration. Similarly, its prior probability density function can be obtained as follows: (18) in, Slowness value vector The covariance matrix.

[0068] According to Bayes' theorem, solving the problem of a three-dimensional sound velocity distribution image is equivalent to maximizing... Posterior probability: (19) Through a series of derivations, we can obtain... The approximate iterative formula is: (20) Because calculating the inverse matrix of a large-dimensional matrix is ​​difficult and time-consuming, formula (20) can be decomposed as follows: ;(twenty one) .(twenty two) First, use the conjugate gradient descent method to find an approximate solution for equation (21). Substituting this into formula (22) yields the distribution vector of the slowness value in this iteration. .

[0069] After iteratively solving for the slowness distribution vector that satisfies the set convergence threshold, taking the reciprocal yields the current ray path matrix. A three-dimensional sound velocity distribution image.

[0070] For the ray path matrix With a fixed value, the slowness value The iterative solution process is called the inner loop. When the solution result of the inner loop makes the objective function reach the preset convergence threshold... At this point, the inner loop ends, the current sound velocity reconstruction result is input into the multi-template fast traversal method to obtain the updated propagation time distribution, and the ray path matrix is ​​updated accordingly. Update the ray path matrix. The process of iteratively solving the problem within the inner loop is called the outer loop. When the solution obtained from the outer loop causes the objective function to reach a preset convergence threshold... At this point, the outer loop ends, and the final three-dimensional sound velocity distribution image is obtained.

[0071] V. Results Output and Quality Verification.

[0072] Two complex three-dimensional numerical phantoms (Shepp-Logan phantom and bronchial phantom) are set as the three-dimensional ground truth (GT) distribution of sound velocity to be reconstructed, such as... Figure 5 As shown in the figure, the bronchial phantom is mainly used to test the ability of the two methods to handle situations with drastic changes in sound velocity in space. Iterative reconstruction of three-dimensional sound velocity images was performed using curved ray models based on the traditional fast-moving model (FMM) and the multi-template fast-moving model (MSFM), respectively. Finally, the results of the two methods were quantitatively compared, demonstrating that the multi-template fast-moving model (MSFM) can improve the quality of reconstructed sound velocity images from three-dimensional ultrasound tomography. Figure 5 Image (a) shows a three-dimensional schematic diagram of the Shepp-Logan phantom and three representative two-dimensional cross-sections. Figure 5 (b) in the figure is a three-dimensional schematic diagram of the bronchial phantom and three representative two-dimensional cross-sections.

[0073] Reconstructed images of the Shepp-Logan phantom and bronchial phantom are as follows: Figure 6 , Figure 7 As shown (a representative two-dimensional cross-section within a three-dimensional volume), the sound velocity reconstruction results demonstrate that, compared to the traditional Fast Marching Method (FMM), the Multi-Template Fast Marching Method (MSFM) significantly suppresses reconstruction artifacts and noise in three-dimensional ultrasound tomography sound velocity images.

[0074] Figure 8 It shows the transverse central axis ( Figure 6 and Figure 7 The figure shows the sound velocity distribution curve (indicated by the dashed line). It can be observed that in regions where the sound velocity changes drastically, the reconstruction results based on the Multi-Template Fast Traveling Method (MSFM) are significantly more consistent with the true values ​​than the traditional Fast Traveling Method (FMM). This phenomenon is due to the fact that when ultrasound waves pass through interfaces with large acoustic impedance differences, a large refraction angle is generated, leading to more significant and complex changes in the propagation trajectory. The core advantage of the Multi-Template Fast Traveling Method (MSFM) lies in using multi-directional templates to solve the equations, which improves the computational accuracy in multiple directions in three-dimensional space, thus more accurately modeling such complex wave propagation behavior. Therefore, the Multi-Template Fast Traveling Method (MSFM) significantly improves the fidelity of sound velocity reconstruction, especially at tissue boundaries with severe acoustic impedance mismatch. Figure 8 (a) in the image is a Shepp-Logan phantom reconstruction image along... Figure 6 The sound velocity value curve is distributed by dashed lines. Figure 8 Image (b) is a reconstructed image of the bronchus phantom along... Figure 7 The sound velocity curve is distributed by the dashed line.

[0075] The root mean square error (RMSE) and peak signal-to-noise ratio (PSNR) were used as two evaluation metrics to quantitatively evaluate the three-dimensional sound velocity images reconstructed by the two methods. The results are shown in Table 2, with the optimal values ​​in Table 2 marked by underline.

[0076] Table 2. Comparison of quantitative evaluation metrics for reconstructing sound velocity images based on the Fast Moving Model (FMM) and Multi-Template Fast Moving Model (MSFM) methods:

[0077] It can be observed that as the reconstruction proceeds iteratively, the RMSE of the Multi-Template Fast Marching Method (MSFM) is significantly lower than that of the Traditional Fast Marching Method (FMM), while the PSNR is significantly higher. For the reconstruction of the Shepp-Logan phantom, MSFM reduced the RMSE from 37.62 to 29.86, while increasing the PSNR from 16.62 to 18.63; for the reconstruction of the bronchial phantom, MSFM reduced the RMSE from 57.61 to 30.87, and increased the PSNR from 12.92 to 18.34.

[0078] It is worth noting that in some cases, the reconstruction results of the traditional Fast Marching Method (FMM) gradually diverge with iterations, while the Multi-Template Fast Marching Method (MSFM) maintains stable convergence throughout the entire reconstruction process. This robust convergence characteristic ensures that 3D ultrasound tomography can consistently output high-quality imaging results, especially exhibiting superior stability when dealing with complex tissue boundaries and high-contrast interfaces. The multi-directional template strategy employed by the Multi-Template Fast Marching Method (MSFM) effectively suppresses the error accumulation problem commonly found in the traditional Fast Marching Method (FMM), thus providing more reliable quantitative results during continuous iterations. In clinical applications, the stability of reconstruction results is crucial for diagnostic reliability; even small fluctuations in results can affect treatment decisions. Therefore, this advantage of the Multi-Template Fast Marching Method (MSFM) has significant clinical implications.

[0079] This invention proposes a scheme that combines "using second-order finite difference to approximate the partial derivatives of the discrete process function equation" with "multi-template design that combines 4 directional templates to cover 18 neighborhoods" to improve the solution accuracy of the process function equation in the curved ray model. This can improve the reconstruction quality of the three-dimensional ultrasonic tomography sound velocity image (reduce artifacts and noise, improve boundary fidelity and image accuracy) and improve the stability of the iterative reconstruction process.

[0080] The terminology used herein is for the purpose of describing particular embodiments only and is not intended to limit the invention. The terms “comprising,” “including,” etc., as used herein indicate the presence of the stated features, steps, operations, and / or components, but do not exclude the presence or addition of one or more other features, steps, operations, or components.

[0081] It should be understood that although the steps in the flowcharts of the accompanying drawings are shown sequentially according to the arrows, these steps are not necessarily executed in the order indicated by the arrows. Unless explicitly stated in this invention, there is no strict order restriction on the execution of these steps, and they can be executed in other orders. Moreover, at least some of the steps in the flowcharts of the accompanying drawings may include multiple steps or multiple stages, which are not necessarily completed at the same time, but may be executed at different times. The execution order of these steps or stages is not necessarily sequential, but may be performed alternately or in turn with other steps or at least a portion of the steps or stages of other steps.

[0082] The technical features of the above embodiments can be combined in any way. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

[0083] 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, and no reference numerals in the claims should be construed as limiting the scope of the claims.

[0084] Furthermore, it should be understood that although this specification describes embodiments, not every embodiment contains only one independent technical solution. This narrative style is merely for clarity. Those skilled in the art should consider the specification as a whole, and the technical solutions in each embodiment can also be appropriately combined to form other embodiments that can be understood by those skilled in the art.

Claims

1. A three-dimensional ultrasonic tomographic velocity imaging method based on a multi-template rapid travel method, characterized in that, include: The array parameters of the three-dimensional ultrasound tomography system and the sound velocity distribution of the medium in which the biological tissue is located are obtained. The computational domain is discretized into a grid, and a curved ray model based on the equation of process function is established. The multi-template fast travel method is used to solve the equation: the propagation time gradient in the equation is discretized by the second-order finite difference approximation, and the 18 neighborhoods of the grid are covered by four directional templates for joint solution to obtain the ultrasonic propagation time distribution of the whole grid. Based on the ultrasonic wave propagation time distribution, a three-dimensional sound velocity distribution image is reconstructed iteratively using the maximum a posteriori probability estimation algorithm.

2. The three-dimensional ultrasonic tomographic velocity imaging method based on multi-template rapid travel method according to claim 1, characterized in that, The equation is: ； Indicates coordinates as The propagation time of the grid gradient, Coordinates are The sound velocity value of the grid.

3. The three-dimensional ultrasonic tomographic velocity imaging method based on multi-template rapid travel method according to claim 1, characterized in that, The method of solving the equation using the multi-template fast travel method also includes a process of classifying the grid based on whether the propagation time has been determined: The grid within the computational domain is divided into a frozen set, a narrowband set, and an unknown set; The propagation time of the grid within the frozen set has been determined and will not be updated. It serves as the benchmark for solving the propagation time from the sound source to other grids. Initially, only the grid containing the sound source belongs to the frozen set. The grid within the narrowband set is the boundary grid between the frozen set and the unknown set. The propagation time from the sound source point to the current grid has been calculated using adjacent grids belonging to the frozen set. The propagation time from the sound source point to the current grid has not yet been calculated for the grid in the unknown set. After the grid in the narrowband set is processed, the grid in the unknown set will be gradually included in the solution range of the equation. In solving the equations, the grid with the shortest propagation time is always selected from the narrowband set for priority processing.

4. The three-dimensional ultrasonic tomographic velocity imaging method based on multi-template rapid travel method according to claim 1, characterized in that, The discretization of the propagation time gradient in the equation using the second-order finite difference approximation specifically includes: ； Indicates the time of transmission. Indicates coordinates as The propagation time of the grid, Indicates coordinates as The propagation time of the grid, Indicates coordinates as The propagation time of the grid, Indicates coordinates as The propagation time of the grid, Indicates coordinates as The propagation time of the grid, Indicates coordinates as The propagation time of the grid, Indicates coordinates as The propagation time of the grid, These are the x-axis, y-axis, and z-axis coordinates of the grid, respectively. These represent the spatial grid step size along the x-axis, y-axis, and z-axis, respectively.

5. The three-dimensional ultrasonic tomographic velocity imaging method based on multi-template rapid travel method according to claim 1, characterized in that, The 18 neighborhoods of the grid include 6 orthogonal neighborhoods and 12 diagonal neighborhoods.

6. The three-dimensional ultrasonic tomographic velocity imaging method based on multi-template rapid travel method according to claim 1, characterized in that, The method employs a four-directional template covering 18 neighborhoods of the mesh for joint solution to obtain the ultrasonic wave propagation time distribution of the entire mesh, specifically including: Adjust the finite difference approximation of the derivative of the fast travel method according to a specific direction template: ；（9） in, This represents the approximate step size of the difference in the template for the current direction. Indicates the propagation time of the current grid. This represents the propagation time sign substitution formula for adjacent grids in the three differential directions corresponding to the current direction template. This represents the three difference directions corresponding to the current direction template. This represents the coordinate vector of the current grid. This represents the coordinate vectors of adjacent grids in the three differential directions corresponding to the current direction template. This represents the current sound speed value in the grid. Transform equation (9) into a quadratic equation and calculate the numerical solution; Numerical solutions are calculated for each direction template, and the minimum value among these numerical solutions is selected as the final solution of the equation for the current grid. .

7. The three-dimensional ultrasonic tomographic velocity imaging method based on multi-template rapid travel method according to claim 6, characterized in that, The quadratic equation is in the form of: ； in, These are the coefficients corresponding to the current direction template when solving in the difference direction u. The spatial grid step size, coefficient ,coefficient ,coefficient , The specific form is determined by the difference approximation order and the direction template.

8. The three-dimensional ultrasonic tomographic velocity imaging method based on multi-template rapid travel method according to claim 1, characterized in that, The process of iteratively reconstructing the three-dimensional sound velocity distribution image based on the ultrasonic wave propagation time distribution using a maximum a posteriori probability estimation algorithm specifically includes: The imaging problem is transformed into an optimization problem that maximizes the posterior probability. The ultrasonic propagation time distribution and the prior probability of the estimated independent variable are combined using Bayes' theorem to iteratively solve the optimization problem. The sound velocity distribution is updated using the conjugate gradient method, thereby suppressing measurement noise and stabilizing the reconstruction process using prior constraints, ultimately obtaining a three-dimensional sound velocity distribution image.

9. A three-dimensional ultrasonic tomographic velocity imaging method based on a multi-template rapid travel method according to claim 8, characterized in that, The transformation of the imaging problem into an optimization problem of maximizing the posterior probability specifically includes: In the curved ray model, Indicates the first The slowness value of each grid cell, where slowness is the reciprocal of the speed of sound, and the grid index is... Arranged linearly in rows; Indicates the index of the ray; Indicates the first The ray passed through the first The intercept size of each grid; Indicates the first The propagation time of a ray; the sound wave propagation process in ultrasound tomography is represented as: ；（12） Transform equation (12) into matrix form: ； The ray path matrix, Let be the time vector of ray propagation. This is the distribution vector of slowness values. The total number of grid cells. The total number of rays; The problem with three-dimensional sound velocity distribution images is that they are observed... In the case of iterative solution Minimize the value of the objective function: ； Let be the objective function. It is the F2 norm.

10. A three-dimensional ultrasonic tomographic velocity imaging method based on a multi-template rapid travel method according to claim 9, characterized in that, The process involves combining the ultrasonic wave propagation time distribution with the prior probability of the estimated independent variable using Bayes' theorem, iteratively solving the optimization problem, and updating the sound velocity distribution using the conjugate gradient method. This suppresses measurement noise while stabilizing the reconstruction process using prior constraints, ultimately obtaining a three-dimensional sound velocity distribution image. Specifically, this includes: set up For the first The unknown noise during the ray measurement process follows a Gaussian distribution. And they are independent of each other. For the first The vector formed by the intercepts of the rays on all grid lines. for The standard deviation of the Gaussian distribution it follows If it is a Gaussian distribution, then probability density function for: ； Similarly, conditional probability density function for: ； Then the joint probability density function for: ； in, Noise vector The covariance matrix, Represents a determinant. Represents the identity matrix. Represents matrix transpose and noise vector ; Let the slowness value distribution vector be... Follows Gaussian distribution ,in, Indicates the first The distribution vector of the slowness value in the next iteration, then Prior probability density function for: ； in, The distribution vector of slowness values The covariance matrix, for The standard deviation of the Gaussian distribution it follows; According to Bayes' theorem, solving the problem of finding a three-dimensional sound velocity distribution image is equivalent to maximizing... posterior probability : ； After derivation, we obtain Approximate iterative formula: ； The current ray path matrix is ​​obtained by solving the conjugate gradient descent method. Below, the distribution vector of the slowness value in the (k+1)th iteration. ; After iteratively solving for the slowness distribution vector that satisfies the set convergence threshold, taking the reciprocal yields the current ray path matrix. A three-dimensional sound velocity distribution image.

Images