阅读说明：本技术 一种基于费马原理的超声ct图像重建方法及系统 (Ultrasonic CT image reconstruction method and system based on Fermat principle ) 是由 尉迟明 丁明跃 方小悦 武云 宋俊杰 周亮 张求德 周权 于 2019-08-27 设计创作，主要内容包括：本发明属于功能成像技术领域,公开了一种基于费马原理的超声CT图像重建方法及系统,其中基于费马原理的超声CT声速重建方法包括步骤：(1)渡越时间的提取；(2)准备进入迭代；(3)根据有限差分方法计算从每个发射阵元出发到成像区域每个像素点的时间τ,并根据费马原理计算胖射线路径,其中胖射线路径的宽度随迭代次数的增加而收窄；(4)反问题的求解,并更新声速值；(5)判断迭代是否终止。本发明通过对方法整体流程进行改进,通过路径的优化,基于费马原理,尤其是变参数的、路径由宽至窄的变路径宽度的胖射线路径,使得本发明中超声CT声速及衰减系数重建方法及系统具有快速、稳定、成像效果更好的特点。(The invention belongs to the technical field of functional imaging and discloses an ultrasonic CT image reconstruction method and system based on the Fermat principle, wherein the ultrasonic CT sound velocity reconstruction method based on the Fermat principle comprises the following steps: (1) extracting the transit time; (2) preparing to enter iteration; (3) calculating the time tau from each emission array element to each pixel point of an imaging region according to a finite difference method, and calculating a fat ray path according to the Fermat principle, wherein the width of the fat ray path is narrowed along with the increase of the iteration times; (4) solving an inverse problem and updating a sound velocity value; (5) it is determined whether the iteration is terminated. The method and the system for reconstructing the ultrasonic CT sound velocity and the attenuation coefficient have the characteristics of rapidness, stability and better imaging effect by improving the overall process of the method, optimizing the path and based on the Fermat principle, particularly the fat ray path with the variable parameter and the variable path width from wide to narrow.)

1. an ultrasonic CT sound velocity reconstruction method based on the Fermat principle is characterized by comprising the following steps:

(1) Extracting the transit time:

Aiming at an object to be detected, acquiring transmission wave data received by all receiving array elements within a preset receiving array element angle requirement range by adopting an annular array probe according to a preset sequence of the transmitting array elements and based on one transmitting array element; then, according to the preset sequence of the receiving array elements, for each piece of transmitted wave data, selecting a matching window near the estimated transition time point on the transmitted wave data, presetting the window length as N, then, taking the current retrieval point as a demarcation point, dividing the window into two sections, calculating the AIC value of the current retrieval point by using a formula (1), then, sequentially retrieving the points in the window, and recording the point with the minimum AIC value as the transition time point;

AIC(k)＝k*log(var(d(1，k)))+(N-k-1)*log(var(d(k+1，N))) (1)

In the formula (1), var represents a variance, and d (1, k) represents data from the first point to the kth point;

thus obtaining a series of transit times corresponding to the same transmitting array element and according to the sequence of the receiving array element, finally obtaining a total set of the transit times corresponding to all the transmitting array elements and according to the sequence of the transmitting array element and the receiving array element, vectorizing the total set of the transit times to obtain a transit time column vector [ T_tof]；

(2) giving a preset sound velocity value as an initial sound velocity value, recording the iteration number as 1, and entering first iteration;

(3) Calculating the time tau from each emission array element to each pixel point of an imaging area according to a finite difference method, and calculating a fat ray path according to the Fermat principle, wherein the method specifically comprises the following steps:

dividing an imaging region into m₁×m₂A grid of m₁、m₂Are all preset positive integers;

Then, according to the preset sequence of the transmitting array elements, for one transmitting array element, selecting one receiving array element according to the preset sequence of the receiving array elements from all the receiving array elements which correspond to the transmitting array element and meet the angle requirement range of the preset receiving array element to obtain a pair of transmitting array elements-receiving array elements; then, calculating the fat ray path of the sound wave from the transmitting array element to the receiving array element by the pair of transmitting array element-receiving array element; specifically, a wired difference method is adopted to calculate the time tau of each transmitting array element transmitting sound wave and spreading to each pixel point in an imaging area, each pixel point P in the imaging area is taken as an object, and if the time tau meets the requirement that each transmitting array element transmits the sound wave, each pixel point P in the imaging area is taken as an object

τ_SP+τ_RP-τ_SR≤Δt (2)

then, the fat ray path value of the grid is a valid value and is recorded as a preset value; otherwise, the value is an invalid value and is marked as 0; thus obtaining a series of fat ray path values corresponding to each grid of the imaging area under the pair of transmitting array elements and receiving array elements, and then vectorizing the series of fat ray path values to obtain fat ray path row vectors;

In the formula (2), S represents a transmitting array element, R represents a receiving array element, S and R are mutually transmitting and receiving, and tau_SPRepresents the time of departure of the sound wave from S to P, tau_RPrepresents the time of departure of the sound wave from R to P, tau_SRrepresenting the time of the sound wave from S to R; delta t is the ith time threshold selected from a preset time threshold sequence { delta t } according to the iteration times, i is equal to the iteration times in value, the time thresholds in the time threshold sequence { delta t } are arranged in a descending order, and any time threshold in the time threshold sequence { delta t } is related to the central frequency f of the probe;

thus, a series of fat ray path row vectors corresponding to the same transmitting array element and according to the sequence of the receiving array element can be obtained, and finally, a fat ray path row vector total set corresponding to all transmitting array elements and according to the sequence of the transmitting array element and the receiving array element is obtained; regarding each row vector as an element whole, arranging all elements in the fat ray path row vector total set into a column, and obtaining a fat ray path total matrix [ L ] with the row position corresponding to the sequence of the transmitting array element and the receiving array element;

(4) Solving an inverse problem, and updating a sound speed value:

Based on the transition time column vector [ T ] obtained in the step (1)_tof]And the fat ray path total matrix [ L ] obtained in the step (3)]Constructing a path-slowness-time equation system shown in the formula (3):

[L][S]＝[T_tof](3)

In the formula (3), [ S ]]Is the total number of rows m₁×m₂The slowness column vector to be solved with the total column number of 1;

then solving the equation set by adopting a random gradient descent method to obtain [ S ], then taking the reciprocal of each element in the [ S ] to obtain a speed reconstruction value aiming at the object to be detected, and then updating the sound velocity value by the speed reconstruction value; updating the iteration times at the same time, and adding 1 to the iteration times;

(5) judging whether the iteration times are smaller than a preset iteration time threshold value or not, and if so, repeating the step (3) and the step (4); if not, terminating iteration and outputting a final sound velocity value;

in addition, the total number of time thresholds in the preset time threshold number series { Δ t } is equal to or greater than the preset iteration number threshold.

2. the Fermat principle-based ultrasonic CT sound velocity reconstruction method according to claim 1, wherein in step (1), the predetermined receiving array element angle requirement range is based on the right opposite side of the transmitting array element, the annular array, the center of the probe as the center of the circle, and the central angle of the probe from- α to + α to form a receiving array element within an angle range of 2 α, wherein α is a predetermined positive angle, and 2 α satisfies 90 ° to 270 °.

3. The Fermat principle-based ultrasonic CT sound velocity reconstruction method as claimed in claim 1, wherein in the step (2), the value of sound velocity in water is used as the initial sound velocity value;

In the step (3), the effective value of the fat transmission path value of the mesh is recorded as 1.

4. An ultrasonic CT attenuation coefficient reconstruction method based on the Fermat principle is characterized by comprising the following steps:

(1) Extracting energy parameters:

Aiming at an object to be detected, acquiring energy parameters of transmitted waves received by all receiving array elements within a preset receiving array element angle requirement range by adopting an annular array probe according to a preset sequence of the transmitting array elements and based on one transmitting array element, and arranging the energy parameters according to the preset sequence of the receiving array elements, thereby obtaining a series of energy parameters corresponding to the same transmitting array element and according to the sequence of the receiving array elements; finally, obtaining an energy parameter total set aiming at all transmitting array elements and corresponding to the transmitting array elements and the receiving array elements in sequence, and vectorizing the parameter total set to obtain an energy parameter column vector [ P ];

(2) giving a preset attenuation coefficient value as an initial attenuation coefficient value, recording the iteration times as 1, and entering first iteration;

(3) calculating the time tau from each emission array element to each pixel point of an imaging area according to a finite difference method, and calculating a fat ray path according to the Fermat principle, wherein the method specifically comprises the following steps:

dividing an imaging region into m₁×m₂A grid of m₁、m₂Are all preset positive integers;

Then, according to the preset sequence of the transmitting array elements, for one transmitting array element, selecting one receiving array element according to the preset sequence of the receiving array elements from all the receiving array elements which correspond to the transmitting array element and meet the angle requirement range of the preset receiving array element to obtain a pair of transmitting array elements-receiving array elements; then, calculating the fat ray path of the sound wave from the transmitting array element to the receiving array element by the pair of transmitting array element-receiving array element; specifically, each grid center point P of the imaging area is taken as an object, and if the requirement is met, the point P is used as the object

τ_SP+τ_RP-τ_SR≤Δt (2)

then, the fat ray path value of the grid is a valid value and is recorded as a preset value; otherwise, the value is an invalid value and is marked as 0; thus obtaining a series of fat ray path values corresponding to each grid of the imaging area under the pair of transmitting array elements and receiving array elements, and then vectorizing the series of fat ray path values to obtain fat ray path row vectors;

in the formula (2), S represents a transmitting array element, R represents a receiving array element, S and R are mutually transmitting and receiving, and tau_SPRepresents the time of departure of the sound wave from S to P, tau_RPRepresents the time of departure of the sound wave from R to P, tau_SRrepresenting the time of the sound wave from S to R; delta t is the ith time threshold selected from a preset time threshold sequence { delta t } according to the iteration times, i is equal to the iteration times in value, the time thresholds in the time threshold sequence { delta t } are arranged in a descending order, and any time threshold in the time threshold sequence { delta t } is related to the central frequency f of the probe;

Thus, a series of fat ray path row vectors corresponding to the same transmitting array element and according to the sequence of the receiving array element can be obtained, and finally, a fat ray path row vector total set corresponding to all transmitting array elements and according to the sequence of the transmitting array element and the receiving array element is obtained; regarding each row vector as an element whole, arranging all elements in the fat ray path row vector total set into a column, and obtaining a fat ray path total matrix [ L ] with the row position corresponding to the sequence of the transmitting array element and the receiving array element;

(4) Solving the inverse problem and updating the attenuation coefficient value:

Constructing a path-attenuation-energy parameter equation set shown in a formula (4) based on the energy parameter column vector [ P ] obtained in the step (1) and the fat ray path total matrix [ L ] obtained in the step (3):

[L][A]＝[P] (4)

in the formula (4), [ A ]]Is the total number of rows m₁×m₂The total column number is 1, and the column vector of the attenuation coefficient to be solved;

then, solving the equation set by adopting a random gradient descent method to obtain [ A ], thus obtaining an attenuation coefficient reconstruction value aiming at the object to be detected, and then updating the attenuation coefficient value by the attenuation coefficient reconstruction value; updating the iteration times at the same time, and adding 1 to the iteration times;

(5) judging whether the iteration times are smaller than a preset iteration time threshold value or not, and if so, repeating the step (3) and the step (4); if not, terminating iteration and outputting a final attenuation coefficient value;

in addition, the total number of time thresholds in the preset time threshold number series { Δ t } is equal to or greater than the preset iteration number threshold.

5. the Fermat principle-based ultrasonic CT attenuation coefficient reconstruction method as claimed in claim 4, wherein in step (1), the predetermined receiving array element angle requirement range is based on the receiving array element which is directly opposite to the transmitting array element, is centered at the center of the ring array probe, and has a central angle ranging from- α to + α to 2 α, wherein α is a predetermined positive angle, and 2 α satisfies 90 ° to 270 °; preferably, 2 α is 180 °.

6. the Fermat principle-based ultrasonic CT attenuation coefficient reconstruction method as claimed in claim 4, wherein in the step (2), the initial attenuation coefficient value is the attenuation coefficient value in water;

in the step (3), the effective value of the fat transmission path value of the mesh is recorded as 1.

7. An ultrasonic CT image reconstruction method using the ultrasonic CT sound velocity reconstruction method based on the Fermat principle according to any one of claims 1 to 3, characterized in that the method using the ultrasonic CT sound velocity reconstruction method based on the Fermat principle further comprises the following steps:

(6) Imaging: arranging the final sound velocity values obtained in the step (5) into a two-dimensional matrix to form m₁×m₂then, a two-dimensional pixel map is obtained based on the two-dimensional matrix of the sound velocity values, wherein each pixel in the two-dimensional pixel map corresponds to the sound velocity value;

the two-dimensional pixel map is obtained by carrying out logarithmic compression, gray mapping and final display on the sound velocity value.

8. an ultrasound CT image reconstruction method using the ultrasound CT attenuation coefficient reconstruction method based on Fermat principle as claimed in any one of claims 4-6, characterized in that the method using the ultrasound CT attenuation coefficient reconstruction method based on Fermat principle further comprises the following steps:

(6) Imaging: arranging the final attenuation coefficient values obtained in the step (5) into a two-dimensional matrix to form m₁×m₂Then obtaining a two-dimensional pixel map based on the two-dimensional matrix, wherein each pixel in the two-dimensional pixel map corresponds to the attenuation coefficient value;

The two-dimensional pixel map is obtained by carrying out logarithmic compression, gray scale mapping and final display on the attenuation coefficient value.

9. An ultrasonic CT sound velocity reconstruction system based on the Fermat principle, which is characterized by comprising:

The transit time extraction module is to: aiming at an object to be detected, acquiring transmission wave data received by all receiving array elements within a preset receiving array element angle requirement range by adopting an annular array probe according to a preset sequence of the transmitting array elements and based on one transmitting array element; according to the preset sequence of the receiving array elements, for each piece of transmitted wave data, selecting a matching window near a predicted transition time point on the transmitted wave data, presetting the window length as N, taking a current retrieval point as a demarcation point, dividing the window into two sections, calculating the AIC value of the current retrieval point by using a formula (1), sequentially retrieving the points in the window, and recording the point with the minimum AIC value as the transition time point;

AIC(k)＝k*log(var(d(1，k)))+(N-k-1)*log(var(d(k+1，N))) (1)

in the formula (1), var represents a variance, and d (1, k) represents data from the first point to the kth point;

Thereby obtaining a series of transit times corresponding to the same transmitting array element and according to the sequence of the receiving array element; finally, a total transition time set corresponding to all the transmitting array elements according to the sequence of the transmitting array elements and the receiving array elements is obtained, and the total transition time set is vectorized to obtain a transition time column vector [ T_tof]；

The iterative computation function module comprises a fat ray path computation module and an inverse problem solving module, wherein,

The fat ray path computation module is configured to: calculating the time tau from each emission array element to each pixel point of an imaging area according to a finite difference method, and calculating a fat ray path according to the Fermat principle, wherein the method specifically comprises the following steps:

dividing an imaging region into m₁×m₂a grid of m₁、m₂are all preset positive integers;

according to the preset sequence of the transmitting array elements, for one transmitting array element, selecting one receiving array element according to the preset sequence of the receiving array elements from all the receiving array elements which correspond to the transmitting array element and meet the angle requirement range of the preset receiving array element to obtain a pair of transmitting array elements-receiving array elements; calculating the fat ray path of the sound wave from the transmitting array element to the receiving array element by the pair of transmitting array element-receiving array element; specifically, each grid center point P of the imaging area is taken as an object, and if the requirement is met, the point P is used as the object

τ_SP+τ_RP-τ_SR≤Δt (2)

then, the fat ray path value of the grid is a valid value and is recorded as a preset value; otherwise, the value is an invalid value and is marked as 0; obtaining a series of fat ray path values corresponding to each grid of the imaging region under the pair of transmitting array elements and receiving array elements, and performing row vectorization on the series of fat ray path values to obtain fat ray path row vectors;

In the formula (2), S represents a transmitting array element, R represents a receiving array element, S and R are mutually transmitting and receiving, and tau_SPRepresents the time of departure of the sound wave from S to P, tau_RPRepresents the time of departure of the sound wave from R to P, tau_SRrepresenting the time of the sound wave from S to R; delta t is the ith time threshold selected from a preset time threshold sequence { delta t } according to the iteration times, i is equal to the iteration times in value, the time thresholds in the time threshold sequence { delta t } are arranged in a descending order, and any time threshold in the time threshold sequence { delta t } is related to the central frequency f of the probe;

Thereby obtaining a series of fat ray path row vectors corresponding to the same transmitting array element and the receiving array element sequence, and finally obtaining a fat ray path row vector total set corresponding to all transmitting array elements and the transmitting array element and the receiving array element sequence; regarding each row vector as an element whole, arranging all elements in the fat ray path row vector total set into a column to obtain a fat ray path total matrix [ L ] with the row position corresponding to the sequence of the transmitting array element and the receiving array element;

The inverse problem solving module is configured to: based on the transit time column vector [ T_tof]and the fat ray path total matrix [ L]Constructing a path-slowness-time equation system shown in the formula (3):

[L][S]＝[T_tof](3)

In the formula (3), [ S ]]is the total number of rows m₁×m₂the slowness column vector to be solved with the total column number of 1;

solving the equation set by adopting a random gradient descent method to obtain [ S ], and taking the reciprocal of each element in the [ S ] to obtain a speed reconstruction value aiming at the object to be detected, wherein the speed reconstruction value is the updated sound velocity value;

The function module for judging whether to terminate is used for: judging whether the iteration times are smaller than a preset iteration time threshold value or not, and repeating the iteration if the iteration times are smaller than the preset iteration time threshold value; if not, terminating iteration and outputting a final sound velocity value;

In addition, the total number of time thresholds in the preset time threshold number series { Δ t } is equal to or greater than the preset iteration number threshold.

10. an ultrasonic CT attenuation coefficient reconstruction system based on the fermat principle, the system comprising:

The energy parameter extraction module is used for: aiming at an object to be detected, acquiring energy parameters of transmitted waves received by all receiving array elements within a preset receiving array element angle requirement range by adopting an annular array probe according to a preset sequence of the transmitting array elements and based on one transmitting array element, and arranging the energy parameters according to the preset sequence of the receiving array elements to obtain a series of energy parameters corresponding to the same transmitting array element and according to the sequence of the receiving array elements; finally, obtaining an energy parameter total set aiming at all transmitting array elements and corresponding to the transmitting array elements and the receiving array elements in sequence, and vectorizing the parameter total set to obtain an energy parameter column vector [ P ];

the iterative computation function module comprises a fat ray path computation module and an inverse problem solving module, wherein,

the fat ray path computation module is configured to: calculating the time tau from each emission array element to each pixel point of an imaging area according to a finite difference method, and calculating a fat ray path according to the Fermat principle, wherein the method specifically comprises the following steps:

dividing an imaging region into m₁×m₂a grid of m₁、m₂Are all preset positive integers;

According to the preset sequence of the transmitting array elements, for one transmitting array element, selecting one receiving array element according to the preset sequence of the receiving array elements from all the receiving array elements which correspond to the transmitting array element and meet the angle requirement range of the preset receiving array element to obtain a pair of transmitting array elements-receiving array elements; calculating the fat ray path of the sound wave from the transmitting array element to the receiving array element by the pair of transmitting array element-receiving array element; specifically, each grid center point P of the imaging area is taken as an object, and if the requirement is met, the point P is used as the object

τ_SP+τ_RP-τ_SR≤Δt (2)

then, the fat ray path value of the grid is a valid value and is recorded as a preset value; otherwise, the value is an invalid value and is marked as 0; thus obtaining a series of fat ray path values corresponding to each grid of the imaging region under the pair of transmitting array elements and receiving array elements, and carrying out row vectorization on the series of fat ray path values to obtain fat ray path row vectors;

In the formula (2), S represents a transmitting array element, R represents a receiving array element, S and R are mutually transmitting and receiving, and tau_SPRepresents the time of departure of the sound wave from S to P, tau_RPrepresents the time of departure of the sound wave from R to P, tau_SRrepresenting the time of the sound wave from S to R; delta t is the ith time threshold selected from a preset time threshold sequence { delta t } according to the iteration times, i is equal to the iteration times in value, the time thresholds in the time threshold sequence { delta t } are arranged in a descending order, and any time threshold in the time threshold sequence { delta t } is related to the central frequency f of the probe;

thereby obtaining a series of fat ray path row vectors corresponding to the same transmitting array element and the receiving array element sequence, and finally obtaining a fat ray path row vector total set corresponding to all transmitting array elements and the transmitting array element and the receiving array element sequence; regarding each row vector as an element whole, arranging all elements in the fat ray path row vector total set into a column, and obtaining a fat ray path total matrix [ L ] with the row position corresponding to the sequence of the transmitting array element and the receiving array element;

the inverse problem solving module is used for: constructing a path-attenuation-energy parameter equation set shown in a formula (4) based on the energy parameter column vector [ P ] and the fat ray path total matrix [ L ]:

[L][A]＝[P](4)

In the formula (4), [ A ]]is the total number of rows m₁×m₂the total column number is 1, and the column vector of the attenuation coefficient to be solved;

solving the equation set by adopting a random gradient descent method to obtain [ A ], so as to obtain an attenuation coefficient reconstruction value aiming at the object to be detected, wherein the attenuation coefficient reconstruction value is the updated attenuation coefficient value;

The function module for judging whether to terminate is used for: judging whether the iteration times are smaller than a preset iteration time threshold value or not, and repeating the iteration if the iteration times are smaller than the preset iteration time threshold value; if not, terminating iteration and outputting a final attenuation coefficient value;

in addition, the total number of time thresholds in the preset time threshold number series { Δ t } is equal to or greater than the preset iteration number threshold.

Technical Field

The invention belongs to the technical field of functional imaging, and particularly relates to an ultrasonic CT image reconstruction method and system based on the Fermat principle, belonging to a transmission type ultrasonic imaging mode in ultrasonic tomography and being capable of realizing reconstruction of ultrasonic CT sound velocity and attenuation coefficient.

Background

Ultrasonic CT is a tomographic imaging mode, which uses an ultrasonic probe to transmit and receive signals, generates reflection data and transmission data, and uses the data to reconstruct ultrasonic tomographic images of different modalities so as to better observe parameter information inside an object. The ultrasonic CT detection has the advantages of low price, no radiation to human bodies and the like, and along with the rapid development of a probe processing technology and high-performance operation of a computer, the ultrasonic tomography technology becomes a research hotspot in the industrial field and the medical field again in recent years.

The data acquired by the ultrasonic CT system can be used for reconstructing various parameter information such as sound velocity, attenuation, density and the like, and belongs to the field of functional images. Here we use transmission data. It has been shown that in the early stage of the disease, the functional parameters of the diseased tissue change earlier than the structural changes. Therefore, the ultrasonic CT functional imaging has important auxiliary significance for early diagnosis of the lesion.

Taking the sound velocity as an example, the ultrasonic CT sound velocity reconstruction method includes a reconstruction algorithm based on a ray theory and a reconstruction algorithm based on a fluctuation theory. The method based on the fluctuation theory has better imaging resolution capability, but is easily disturbed by tiny errors, so the robustness is not high, the calculation amount is large, the requirements on the system precision and the signal-to-noise ratio of data are higher, and the method is not suitable for practical application at present. The reconstruction algorithm model based on the ray theory is simpler, has higher stability and smaller computation amount, and is a high-efficiency and stable sound velocity reconstruction method suitable for clinic at present.

In recent years, research on ultrasonic CT sound velocity reconstruction methods is gradually active at home and abroad. In the research field based on the ray theory, the ray path tracking methods are numerous, the applicability is different for different application scenes, and the imaging effect is different. There are those using a straight path and those using a broken path. The straight-line path is adopted to simulate an X-ray CT reconstruction method, and theoretically, the propagation of the ultrasound is more similar to a broken line. However, due to the limited bandwidth of the clinical ultrasound probe, the propagation path of the ultrasound is not a polyline in practice. The tracing of the path is thus a difficulty of this type of method.

Disclosure of Invention

In view of the above defects or improvement requirements of the prior art, the present invention aims to provide an ultrasound CT image reconstruction method and system based on the fermat principle, wherein the overall process of the method is improved, and by optimizing the path, based on the fermat principle, especially the fat ray path with variable parameters and a path width varying from wide to narrow, on one hand, the method and system for reconstructing ultrasound CT sound velocity and attenuation coefficient in the present invention have the characteristics of rapidness, stability and better imaging effect by considering not only the influence of points on the ray path on receiving information, but also the influence of other points on receiving information in the ray field, and using a specific path tracking process, on the other hand, more importantly, the method for reconstructing from low resolution to deepening high resolution is implemented, so that the reconstruction process is more stable and the reconstruction resolution is higher, reconstruction artifacts are reduced.

to achieve the above object, according to one aspect of the present invention, there is provided an ultrasonic CT sound velocity reconstruction method based on the fermat principle, comprising the steps of:

(1) extracting the transit time:

Aiming at an object to be detected, acquiring transmission wave data received by all receiving array elements within a preset receiving array element angle requirement range by adopting an annular array probe according to a preset sequence of the transmitting array elements and based on one transmitting array element; then, according to the preset sequence of the receiving array elements, for each piece of transmitted wave data, selecting a matching window near the estimated transition time point on the transmitted wave data, presetting the window length as N, then, taking the current retrieval point as a demarcation point, dividing the window into two sections, calculating the AIC value of the current retrieval point by using a formula (1), then, sequentially retrieving the points in the window, and recording the point with the minimum AIC value as the transition time point;

AIC(k)--k*log(var(d(1,k)))+(N-k-1)*log(var(d(k+1，N))) (1)

In the formula (1), var represents a variance, and d (1, k) represents data from the first point to the kth point;

Thus obtaining a series of transit times corresponding to the same transmitting array element and according to the sequence of the receiving array element, finally obtaining a total set of the transit times corresponding to all the transmitting array elements and according to the sequence of the transmitting array element and the receiving array element, vectorizing the total set of the transit times to obtain a transit time column vector [ T_tof]；

(2) Giving a preset sound velocity value as an initial sound velocity value, recording the iteration number as 1, and entering first iteration;

(3) Calculating the time tau from each emission array element to each pixel point of an imaging area according to a finite difference method, and calculating a fat ray path according to the Fermat principle, wherein the method specifically comprises the following steps:

Dividing an imaging region into m₁×m₂a grid of m₁、m₂are all preset positive integers;

Then, according to the preset sequence of the transmitting array elements, for one transmitting array element, selecting one receiving array element according to the preset sequence of the receiving array elements from all the receiving array elements which correspond to the transmitting array element and meet the angle requirement range of the preset receiving array element to obtain a pair of transmitting array elements-receiving array elements; then, calculating the fat ray path of the sound wave from the transmitting array element to the receiving array element by the pair of transmitting array element-receiving array element; specifically, a wired difference method is adopted to calculate the time tau of each transmitting array element transmitting sound wave and spreading to each pixel point in an imaging area, each pixel point P in the imaging area is taken as an object, and if the time tau meets the requirement that each transmitting array element transmits the sound wave, each pixel point P in the imaging area is taken as an object

τ_SP+τ_RP-τ_SR≤Δt (2)

Then, the fat ray path value of the grid is a valid value and is recorded as a preset value; otherwise, the value is an invalid value and is marked as 0; thus obtaining a series of fat ray path values corresponding to each grid of the imaging area under the pair of transmitting array elements and receiving array elements, and then vectorizing the series of fat ray path values to obtain fat ray path row vectors;

In the formula (2), S represents a transmitting array element, R represents a receiving array element, S and R are mutually transmitting and receiving, and tau_SPRepresents the time of departure of the sound wave from S to P, tau_RPRepresents the time of departure of the sound wave from R to P, tau_SRRepresenting the time of the sound wave from S to R; delta t is the ith time threshold selected from a preset time threshold sequence { delta t } according to the iteration times, i is equal to the iteration times in value, the time thresholds in the time threshold sequence { delta t } are arranged in a descending order, and any time threshold in the time threshold sequence { delta t } is related to the central frequency f of the probe;

Thus, a series of fat ray path row vectors corresponding to the same transmitting array element and according to the sequence of the receiving array element can be obtained, and finally, a fat ray path row vector total set corresponding to all transmitting array elements and according to the sequence of the transmitting array element and the receiving array element is obtained; regarding each row vector as an element whole, arranging all elements in the fat ray path row vector total set into a column, and obtaining a fat ray path total matrix [ L ] with the row position corresponding to the sequence of the transmitting array element and the receiving array element;

(4) Solving an inverse problem, and updating a sound speed value:

Based on the transition time column vector [ T ] obtained in the step (1)_tof]And the fat ray path total matrix [ L ] obtained in the step (3)]constructing a path-slowness-time equation system shown in the formula (3):

[L][S]＝[T_tof] (3)

In the formula (3), [ S ]]Is the total number of rows m₁×m₂The slowness column vector to be solved with the total column number of 1;

Then solving the equation set by adopting a random gradient descent method to obtain [ S ], then taking the reciprocal of each element in the [ S ] to obtain a speed reconstruction value aiming at the object to be detected, and then updating the sound velocity value by the speed reconstruction value; updating the iteration times at the same time, and adding 1 to the iteration times;

(5) Judging whether the iteration times are smaller than a preset iteration time threshold value or not, and if so, repeating the step (3) and the step (4); if not, terminating iteration and outputting a final sound velocity value;

in addition, the total number of time thresholds in the preset time threshold number series { Δ t } is equal to or greater than the preset iteration number threshold.

as a further preferred aspect of the present invention, in the step (1), the predetermined receiving array element angle requirement range is a receiving array element in an angle range of 2 α in total, where α is a predetermined positive angle and 2 α satisfies 90 ° to 270 °, based on the opposite face of the transmitting array element, the annular array, the center of the probe, and the central angle of the probe from- α to + α.

as a further preferred aspect of the present invention, in the step (2), a sound velocity value in water is used as the initial sound velocity value;

In the step (3), the effective value of the fat transmission path value of the mesh is recorded as 1.

According to another aspect of the present invention, the present invention provides an ultrasonic CT attenuation coefficient reconstruction method based on the fermat principle, which is characterized by comprising the following steps:

(1) Extracting energy parameters:

aiming at an object to be detected, acquiring energy parameters of transmitted waves received by all receiving array elements within a preset receiving array element angle requirement range by adopting an annular array probe according to a preset sequence of the transmitting array elements and based on one transmitting array element, and arranging the energy parameters according to the preset sequence of the receiving array elements, thereby obtaining a series of energy parameters corresponding to the same transmitting array element and according to the sequence of the receiving array elements; finally, obtaining an energy parameter total set aiming at all transmitting array elements and corresponding to the transmitting array elements and the receiving array elements in sequence, and vectorizing the parameter total set to obtain an energy parameter column vector [ P ];

(2) Giving a preset attenuation coefficient value as an initial attenuation coefficient value, recording the iteration times as 1, and entering first iteration;

(3) Calculating the time tau from each emission array element to each pixel point of an imaging area according to a finite difference method, and calculating a fat ray path according to the Fermat principle, wherein the method specifically comprises the following steps:

Dividing an imaging region into m₁×m₂a grid of m₁、m₂are all preset positive integers;

Then, according to the preset sequence of the transmitting array elements, for one transmitting array element, selecting one receiving array element according to the preset sequence of the receiving array elements from all the receiving array elements which correspond to the transmitting array element and meet the angle requirement range of the preset receiving array element to obtain a pair of transmitting array elements-receiving array elements; then, calculating the fat ray path of the sound wave from the transmitting array element to the receiving array element by the pair of transmitting array element-receiving array element; specifically, each grid center point P of the imaging area is taken as an object, and if the requirement is met, the point P is used as the object

τ_SP+τ_RP-τ_SR≤Δt (2)

Then, the fat ray path value of the grid is a valid value and is recorded as a preset value; otherwise, the value is an invalid value and is marked as 0; thus obtaining a series of fat ray path values corresponding to each grid of the imaging area under the pair of transmitting array elements and receiving array elements, and then vectorizing the series of fat ray path values to obtain fat ray path row vectors;

In the formula (2), S represents a transmitting array element, R represents a receiving array element, S and R are mutually transmitting and receiving, and tau_SPRepresenting the time when the sound wave starts from S and reaches Pm, t_RPRepresents the time of departure of the sound wave from R to P, tau_SRRepresenting the time of the sound wave from S to R; delta t is the ith time threshold selected from a preset time threshold sequence { delta t } according to the iteration times, i is equal to the iteration times in value, the time thresholds in the time threshold sequence { delta t } are arranged in a descending order, and any time threshold in the time threshold sequence { delta t } is related to the central frequency f of the probe;

thus, a series of fat ray path row vectors corresponding to the same transmitting array element and according to the sequence of the receiving array element can be obtained, and finally, a fat ray path row vector total set corresponding to all transmitting array elements and according to the sequence of the transmitting array element and the receiving array element is obtained; regarding each row vector as an element whole, arranging all elements in the fat ray path row vector total set into a column, and obtaining a fat ray path total matrix [ L ] with the row position corresponding to the sequence of the transmitting array element and the receiving array element;

(4) Solving the inverse problem and updating the attenuation coefficient value:

Constructing a path-attenuation-energy parameter equation set shown in a formula (4) based on the energy parameter column vector [ P ] obtained in the step (1) and the fat ray path total matrix [ L ] obtained in the step (3):

[L][A]＝[P] (4)

in the formula (4), [ A ]]Is the total number of rows m₁×m₂The total column number is 1, and the column vector of the attenuation coefficient to be solved;

Then, solving the equation set by adopting a random gradient descent method to obtain [ A ], thus obtaining an attenuation coefficient reconstruction value aiming at the object to be detected, and then updating the attenuation coefficient value by the attenuation coefficient reconstruction value; updating the iteration times at the same time, and adding 1 to the iteration times;

(5) Judging whether the iteration times are smaller than a preset iteration time threshold value or not, and if so, repeating the step (3) and the step (4); if not, terminating iteration and outputting a final attenuation coefficient value;

In addition, the total number of time thresholds in the preset time threshold number series { Δ t } is equal to or greater than the preset iteration number threshold.

As a further preferable aspect of the present invention, in the step (1), the preset receiving array element angle requirement range is a receiving array element with the opposite face of the transmitting array element as a reference, the center of the circular array probe as a center of a circle, and the central angle of the circular array probe from- α to + α to 2 α, where α is a preset positive angle, and 2 α satisfies 90 ° to 270 °; preferably, 2 α is 180 °.

As a further preferred aspect of the present invention, in the step (2), the initial attenuation coefficient value is an attenuation coefficient value in water;

in the step (3), the effective value of the fat transmission path value of the mesh is recorded as 1.

according to another aspect of the present invention, the present invention provides an ultrasound CT image reconstruction method using the above mentioned ultrasound CT sound velocity reconstruction method based on the fermat principle, which is characterized in that the method uses the above mentioned ultrasound CT sound velocity reconstruction method based on the fermat principle, and further includes the following steps:

(6) Imaging: arranging the final sound velocity values obtained in the step (5) into a two-dimensional matrix to form m₁×m₂then, a two-dimensional pixel map is obtained based on the two-dimensional matrix of the sound velocity values, wherein each pixel in the two-dimensional pixel map corresponds to the sound velocity value;

the two-dimensional pixel map is obtained by carrying out logarithmic compression, gray mapping and final display on the sound velocity value.

According to another aspect of the present invention, the present invention provides an ultrasound CT image reconstruction method using the ultrasound CT attenuation coefficient reconstruction method based on the fermat principle, which is characterized in that the method uses the ultrasound CT attenuation coefficient reconstruction method based on the fermat principle, and further includes the following steps:

(6) Imaging: arranging the final attenuation coefficient values obtained in the step (5) into a two-dimensional matrix to form m₁×m₂then obtaining a two-dimensional pixel map based on the two-dimensional matrix, wherein each pixel in the two-dimensional pixel map corresponds to the attenuation coefficient value;

the two-dimensional pixel map is obtained by carrying out logarithmic compression, gray scale mapping and final display on the attenuation coefficient value.

according to another aspect of the present invention, the present invention provides an ultrasonic CT sound velocity reconstruction system based on the fermat principle, which is characterized in that the system comprises:

the transit time extraction module is to: aiming at an object to be detected, acquiring transmission wave data received by all receiving array elements within a preset receiving array element angle requirement range by adopting an annular array probe according to a preset sequence of the transmitting array elements and based on one transmitting array element; according to the preset sequence of the receiving array elements, for each piece of transmitted wave data, selecting a matching window near a predicted transition time point on the transmitted wave data, presetting the window length as N, taking a current retrieval point as a demarcation point, dividing the window into two sections, calculating the AIC value of the current retrieval point by using a formula (1), sequentially retrieving the points in the window, and recording the point with the minimum AIC value as the transition time point;

AIC(k)＝k*log(var(d(1,k)))+(N-k-1)*log(var(d(k+1，N))) (1)

In the formula (1), var represents a variance, and d (1, k) represents data from the first point to the kth point;

thereby obtaining a series of transit times corresponding to the same transmitting array element and according to the sequence of the receiving array element; finally, a total transition time set corresponding to all the transmitting array elements according to the sequence of the transmitting array elements and the receiving array elements is obtained, and the total transition time set is vectorized to obtain a transition time column vector [ T_tof]；

The iterative computation function module comprises a fat ray path computation module and an inverse problem solving module, wherein,

The fat ray path computation module is configured to: calculating the time tau from each emission array element to each pixel point of an imaging area according to a finite difference method, and calculating a fat ray path according to the Fermat principle, wherein the method specifically comprises the following steps:

dividing an imaging region into m₁×m₂A grid of m₁、m₂are all preset positive integers;

According to the preset sequence of the transmitting array elements, for one transmitting array element, selecting one receiving array element according to the preset sequence of the receiving array elements from all the receiving array elements which correspond to the transmitting array element and meet the angle requirement range of the preset receiving array element to obtain a pair of transmitting array elements-receiving array elements; calculating the fat ray path of the sound wave from the transmitting array element to the receiving array element by the pair of transmitting array element-receiving array element; specifically, each grid center point P of the imaging area is taken as an object, and if the requirement is met, the point P is used as the object

τ_SP+τ_RP-τ_SR≤Δt (2)

then, the fat ray path value of the grid is a valid value and is recorded as a preset value; otherwise, the value is an invalid value and is marked as 0; obtaining a series of fat ray path values corresponding to each grid of the imaging region under the pair of transmitting array elements and receiving array elements, and performing row vectorization on the series of fat ray path values to obtain fat ray path row vectors;

in the formula (2), S represents a transmitting array element, R represents a receiving array element, S and R are mutually transmitting and receiving, and tau_SPrepresents the time of departure of the sound wave from S to P, tau_RPRepresents the time of departure of the sound wave from R to P, tau_SRRepresenting the time of the sound wave from S to R; delta t is the ith time threshold selected from a preset time threshold sequence { delta t } according to the iteration times, i is equal to the iteration times in value, the time thresholds in the time threshold sequence { delta t } are arranged in a descending order, and any time threshold in the time threshold sequence { delta t } is related to the central frequency f of the probe;

thereby obtaining a series of fat ray path row vectors corresponding to the same transmitting array element and the receiving array element sequence, and finally obtaining a fat ray path row vector total set corresponding to all transmitting array elements and the transmitting array element and the receiving array element sequence; regarding each row vector as an element whole, arranging all elements in the fat ray path row vector total set into a column to obtain a fat ray path total matrix [ L ] with the row position corresponding to the sequence of the transmitting array element and the receiving array element;

The inverse problem solving module is configured to: based on the transit time column vector [ T_tof]and the fat ray path total matrix [ L]Constructing a path-slowness-time equation system shown in the formula (3):

[L][S]＝[T_tof] (3)

In the formula (3), [ S ]]is the total number of rows m₁×m₂The slowness column vector to be solved with the total column number of 1;

Solving the equation set by adopting a random gradient descent method to obtain [ S ], and taking the reciprocal of each element in the [ S ] to obtain a speed reconstruction value aiming at the object to be detected, wherein the speed reconstruction value is the updated sound velocity value;

the function module for judging whether to terminate is used for: judging whether the iteration times are smaller than a preset iteration time threshold value or not, and repeating the iteration if the iteration times are smaller than the preset iteration time threshold value; if not, terminating iteration and outputting a final sound velocity value;

in addition, the total number of time thresholds in the preset time threshold number series { Δ t } is equal to or greater than the preset iteration number threshold.

According to a final aspect of the present invention, the present invention provides an ultrasound CT attenuation coefficient reconstruction system based on the fermat principle, characterized in that the system comprises:

The energy parameter extraction module is used for: aiming at an object to be detected, acquiring energy parameters of transmitted waves received by all receiving array elements within a preset receiving array element angle requirement range by adopting an annular array probe according to a preset sequence of the transmitting array elements and based on one transmitting array element, and arranging the energy parameters according to the preset sequence of the receiving array elements to obtain a series of energy parameters corresponding to the same transmitting array element and according to the sequence of the receiving array elements; finally, obtaining an energy parameter total set aiming at all transmitting array elements and corresponding to the transmitting array elements and the receiving array elements in sequence, and vectorizing the parameter total set to obtain an energy parameter column vector [ P ];

the iterative computation functional module comprises a fat ray path computation module and an inverse problem solving module, wherein the fat ray path computation module is used for: calculating the time tau from each emission array element to each pixel point of an imaging area according to a finite difference method, and calculating a fat ray path according to the Fermat principle, wherein the method specifically comprises the following steps:

dividing an imaging region into m₁×m₂A grid of m₁、m₂are all preset positive integers;

According to the preset sequence of the transmitting array elements, for one transmitting array element, selecting one receiving array element according to the preset sequence of the receiving array elements from all the receiving array elements which correspond to the transmitting array element and meet the angle requirement range of the preset receiving array element to obtain a pair of transmitting array elements-receiving array elements; calculating the fat ray path of the sound wave from the transmitting array element to the receiving array element by the pair of transmitting array element-receiving array element; specifically, each grid center point P of the imaging area is taken as an object, and if the requirement is met, the point P is used as the object

τ_SP+τ_RP-τ_SR≤Δt (2)

Then, the fat ray path value of the grid is a valid value and is recorded as a preset value; otherwise, the value is an invalid value and is marked as 0; thus obtaining a series of fat ray path values corresponding to each grid of the imaging region under the pair of transmitting array elements and receiving array elements, and carrying out row vectorization on the series of fat ray path values to obtain fat ray path row vectors;

In the formula (2), S represents a transmitting array element, R represents a receiving array element, S and R are mutually transmitting and receiving, and tau_SPrepresents the time of departure of the sound wave from S to P, tau_RPRepresents the time of departure of the sound wave from R to P, tau_SRrepresenting the time of the sound wave from S to R; delta t is the ith time threshold selected from a preset time threshold sequence { delta t } according to the iteration times, i is equal to the iteration times in value, the time thresholds in the time threshold sequence { delta t } are arranged in a descending order, and any time threshold in the time threshold sequence { delta t } is related to the central frequency f of the probe;

thereby obtaining a series of fat ray path row vectors corresponding to the same transmitting array element and the receiving array element sequence, and finally obtaining a fat ray path row vector total set corresponding to all transmitting array elements and the transmitting array element and the receiving array element sequence; regarding each row vector as an element whole, arranging all elements in the fat ray path row vector total set into a column, and obtaining a fat ray path total matrix [ L ] with the row position corresponding to the sequence of the transmitting array element and the receiving array element;

the inverse problem solving module is used for: constructing a path-attenuation-energy parameter equation set shown in a formula (4) based on the energy parameter column vector [ P ] and the fat ray path total matrix [ L ]:

[L][A]＝[P] (4)

in the formula (4), [ A ]]is the total number of rows m₁×m₂The total column number is 1, and the column vector of the attenuation coefficient to be solved;

solving the equation set by adopting a random gradient descent method to obtain [ A ], so as to obtain an attenuation coefficient reconstruction value aiming at the object to be detected, wherein the attenuation coefficient reconstruction value is the updated attenuation coefficient value;

the function module for judging whether to terminate is used for: judging whether the iteration times are smaller than a preset iteration time threshold value or not, and repeating the iteration if the iteration times are smaller than the preset iteration time threshold value; if not, terminating iteration and outputting a final attenuation coefficient value;

In addition, the total number of time thresholds in the preset time threshold number series { Δ t } is equal to or greater than the preset iteration number threshold.

Compared with the prior art, the technical scheme of the invention has the advantages that the influence of points on the ray path on the received information is considered based on the Fermat principle, the influence of other points in the ray field on the received information is also considered, the fat ray path total matrix is obtained, the path tracking process is simple and rapid, and the imaging resolution capability is higher. On the other hand, the invention realizes the deepening reconstruction method from low resolution to high resolution based on the fat ray path with variable parameters and the path width from wide to narrow, thereby ensuring that the reconstruction process is more stable, the reconstruction resolution is higher and the reconstruction artifacts are reduced.

For a certain pair of transmitting array elements-receiving array elements, when a fat ray path through which the acoustic wave passes from the transmitting array element to the receiving array element is judged, specifically based on formula (2), when an imaging point of an imaging region satisfies formula (2), it is judged that the imaging point is located on the fat ray path, and a fat ray path value of the imaging point is marked as an effective value, otherwise, it is marked as an invalid value 0. Moreover, the effective value can be uniformly set to 1 in advance, so that a fat ray path matrix which is not 0 or 1 can be obtained for all the imaging points, and the subsequent calculation processing is convenient.

τ_SP+τ_RP-τ_SR≤Δt (2)

during each iteration, the value of delta t in the formula (2) is changed, the larger the iteration times, the smaller the delta t is, namely, the delta t is a larger value in the initial stage of the iteration, and the delta t is gradually reduced along with the deepening of the iteration; for example, Δ t may be decreased in order from (1/f), …, (1/3f), …, (1/5f), …, (1/10f), …. The larger Δ t, the wider the fat ray path, and conversely, the smaller Δ t, the narrower the fat ray path. As for different imaging objects and different imaging systems, the width of the path directly influences the imaging resolution, the change process of the delta t which is gradually reduced along with the deepening of iteration substantially corresponds to the image reconstruction method of the fat ray path with variable parameters (variable path width), and the method is a method for deepening reconstruction from low resolution to high resolution, so that the reconstruction process is more stable, the reconstruction resolution is higher, and the reconstruction artifacts are reduced. Moreover, the iteration from the wide path to the narrow path adopted by the method is beneficial to avoiding falling into the local optimal solution; the preset iteration number threshold value can be a positive integer (such as 50 times, 100 times and the like) which is more than or equal to 10; if one iteration is adopted and the minimum delta t is directly taken, the obtained solution is the local optimal solution, and the method can avoid falling into the local optimal solution through slow change.

Furthermore, the invention can further control the accuracy of the transmitted wave data by selecting the transmitted wave data received by all the receiving array elements within the range meeting the preset angle requirement of the receiving array elements. Allowing for reception in the vicinity of the transmitting array elementThe array element mainly receives a reflection signal, aiming at any one transmitting array element, the invention preferably collects transmission wave data received by the receiving array element which takes the opposite surface of the transmitting array element as a reference, takes the center of a circle of the annular array probe as the center of the circle and has the central angle ranging from-alpha to + alpha to sum up to 2 alpha, 2 alpha preferably satisfies 90-270 degrees (namely the angle range is at least 90 degrees and at most not more than 270 degrees), thereby obtaining a transition time sequence vector [ T ] of the element row position corresponding to the transmitting array element and the receiving array element sequence by utilizing the preset sequence (such as the anticlockwise sequence or the clockwise sequence) of the transmitting array element and the receiving array element_tof](or energy parameter column vector [ P ]]) And fat ray path total matrix [ L]。

In addition, the equation set is solved by adopting a random gradient descent method, and the random gradient descent method is used as a method for iteratively solving the equation set and is between a gradient descent method and a Newton method, so that the method is a good iterative optimization method. In the invention, because the ultrasonic CT reconstruction data volume is extremely large and the calculation amount is extremely large, the equation set is solved by adopting a random gradient descent method, and in the iterative optimization solving process, one sample is randomly selected for each iteration to update the parameters to be reconstructed, so that the optimization process is accelerated, and the operation burden can be greatly reduced. One sample is an equation in the equation (3) or the equation (4) in the inverse problem solving step.

Generally, the method and the system for reconstructing the sound velocity and the attenuation coefficient based on the Fermat principle have the following advantages: 1. based on the Fermat principle, the path tracking process is simple and quick; 2. the path tracking based on the Fermat principle enables the imaging resolution capability to be higher; 3. the image reconstruction based on the fat ray path with the variable path can realize deepening reconstruction from low resolution to high resolution, so that the reconstruction process is more stable, the reconstruction resolution is higher, and the reconstruction artifacts are reduced; 4. when judging whether the path is effective, preferably adopting non-0, namely 1, so as to bring convenience to subsequent calculation; 5. the inverse problem is solved by adopting a random gradient descent method, the calculation of a Hessian matrix is avoided, the calculation amount is small, and the solving precision is high.

drawings

FIG. 1 is a schematic diagram of the AIC extraction transit time.

Fig. 2 is a schematic diagram of fat ray path based on the fermat principle in the present invention.

FIG. 3 is a schematic diagram of a simulation data real model used in the present invention.

FIG. 4 is a graph of sound velocity reconstructed by the method of the present invention.

Detailed Description

in order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention is described in further detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention. In addition, the technical features involved in the embodiments of the present invention described below may be combined with each other as long as they do not conflict with each other.

according to the ultrasonic CT image reconstruction method and system based on the Fermat principle, the influence of points on a ray path on the information is considered, and other points in the ray field also have the influence on the received information. Meanwhile, the fat ray path adopts a parameter-variable (path width-variable) idea, the path width is continuously narrowed along with the deepening of iteration (namely along with the increasing of iteration times), and the deepening of reconstruction from low resolution to high resolution can be realized.

taking sound velocity reconstruction as an example, the sound velocity reconstruction method based on the Fermat principle in the invention generally comprises the following steps:

(1) the transit time is extracted.

(2) An initial sound speed value, such as the sound speed value of water, is given and a first iteration is entered.

(3) And calculating the time tau from each emission array element to each pixel point of the imaging region according to a finite difference method, and calculating a fat ray path according to the Fermat principle.

(4) And solving the inverse problem and updating the sound velocity value.

(5) And (4) judging whether the iteration times are smaller than a given iteration time threshold, if so, returning to the step (3), otherwise, terminating the iteration and outputting a sound velocity value.

specifically, the method comprises the following steps:

the ultrasonic CT system adopts an annular array probe, and transducer array elements are uniformly distributed in an annular shape. The collected data utilizes a single-shot full-connection mode, namely one array element is transmitted, all the array elements are received, all the array elements are transmitted in sequence, and if T array elements exist, T can be collected²and (6) recording the data.

the data is first subjected to an extraction of the time of flight (tof). The transit time is the time at which the waveform arrives at the receive location, which is the time at which the first-arriving waveform in the received waveform takes off its jump point. The method adopts Akaike Information Criterion (AIC) extraction.

the ray path that the acoustic wave traverses from the transmitting array element to the receiving array element is then calculated. The invention is based on the fermat principle, so that the path is not a ray, but an area. Firstly, whether an imaging point is in a path area is judged according to the Fermat principle and is stored in a path matrix. And finally, establishing and solving a path-slowness-time equation set, namely establishing and solving an inverse problem. The slowness is the reciprocal of the speed of sound, the slowness in the present invention refers to the slowness of the object to be reconstructed, the path matrix may be a matrix of non-zero or one, if the point is on the path, the matrix is 1, otherwise the matrix is 0(1 means that the imaging point is in the path area, and 0 means that the imaging point is not in the path area). Time is the time of flight of the object to be reconstructed.

In addition, in the path processing process, the fat ray path width is different along with the difference of the iteration times during each iteration; the larger the iteration times, the narrower the path, thereby realizing a fat ray path with a variable path (variable path width), and corresponding to a method for deepening the reconstruction from low resolution to high resolution, the reconstruction process is more stable, the reconstruction resolution is higher, and the reconstruction artifacts are reduced.

The following are specific examples: