ULTRASOUND CT IMAGE RECONSTRUCTION METHOD AND SYSTEM BASED ON RAY THEORY

Abstract

The disclosure is related to the technical field of functional imaging, and discloses an ultrasound CT image reconstruction method and system based on ray theory, wherein the method includes an ultrasound CT sound speed reconstruction method and an ultrasound CT attenuation coefficient reconstruction method based on ray theory; the ultrasound CT sound speed reconstruction method based on ray theory includes: (1) extraction of the difference in travel time; (2) calculate the ray path that the acoustic wave passes from the transmitting array element to the receiving array element; (3) solution for inverse problem: by using the quasi-Newton method to solve the path-slowness-time equation system, the speed reconstruction value vector of the object to be measured can be obtained.

Claims

1. An ultrasound CT image reconstruction method using an ultrasound CT sound speed reconstruction method based on ray theory, comprising the following steps: (1) extracting a travel time difference: first, based on the same transmitting array element and receiving array element, ultrasound transmission wave data from pure water and an object to be measured are collected respectively after an ultrasound wave is transmitted, and respective corresponding pure water data and the data of the object to be measured are obtained respectively according to channels; then use an AIC method to extract a travel time of the pure water data and the data is recorded as tof.sub.water; then, determine a matching window, take tof.sub.water as a starting point, and take a time t.sub.water_max at the maximum amplitude of the pure water data as the end of window; a window length is recorded as w, and a window data in the matching window is recorded as W.sub.water; then look for a sliding window whose window length is maintained at w on the data of the object to be measured in the corresponding channel, the window data in the sliding window is recorded as W.sub.object; then, W.sub.object and W.sub.water are correlated to each other so as to calculate a cross-correlation coefficient; adjust a starting point of the sliding window, thereby sliding the sliding window to obtain a series of cross-correlation coefficients, choose the sliding window corresponding to the cross-correlation coefficient with the largest value among the cross-correlation coefficients as a search result of the sliding window, and record the starting point of the search result of the sliding window as tof.sub.object, then the difference in travel time is Δtof=tof.sub.object−tof.sub.water; next, adjust the channel and repeat the extraction process, and finally obtain a series of travel time difference Δtof corresponding to a series of channels; (2) calculate a ray path that acoustic wave passes from the transmitting array element to the receiving array element: record the total number of a series of effective travel time differences Δtof obtained in step (1) as n.sub.t, and an imaging area is divided so that a grid number of the imaging area satisfies Σ×Σ, wherein Σ is a positive integer, when √{square root over (n.sub.t)} is an integer, is equal to √{square root over (n.sub.t)}, when √{square root over (n.sub.t)} is a non-integer, Σ is an integer obtained by ceiling, flooring, or rounding off √{square root over (n.sub.t)}; the imaging area is established in a two-dimensional plane rectangular array element coordinate system corresponding to the transmitting array element and the receiving array element, then, for each group of transmitting array element-receiving array element, a straight line is connected between a coordinate position of the transmitting array element in an array element coordinate system and a coordinate position of the receiving array element in the array element coordinate system, thereby obtaining a path length of the connection line in each grid in the imaging area, and then obtaining a two-dimensional Σ×Σ matrix associated with the path length of each grid, thereafter, the matrix is vectorized to obtain a vector about the path length in each grid, finally, vectors regarding the path length in each grid obtained by the each group of transmitting array element-receiving array element are arranged to form a two-dimensional Σ.sup.2×Σ.sup.2 matrix path matrix L regarding all of the transmitting array element-receiving array element groups; (3) solution for inverse problem: vectorize the effective travel time difference Δtof obtained in step (1) to obtain ΔT; then construct a path-slowness-time equation system as shown in equation (3):
LΔS=ΔT   (3) wherein, ΔS is an amount of change in slowness to be solved; then, a quasi-Newton method is utilized to solve the equation system, and a one-dimensional vector ΔS containing Σ.sup.2 elements is obtained, thereafter, slowness of ultrasound wave in water is added to the amount of change in slowness ΔS, take a reciprocal, and a speed reconstruction value vector of the object to be measured can be obtained.

2. An ultrasound CT image reconstruction method using an ultrasound CT attenuation coefficient reconstruction method based on ray theory, comprising the following steps: (1) first, based on the same transmitting array element and receiving array element, energy parameters from pure water and the object to be measured are collected respectively after the ultrasound wave is transmitted, and the respective corresponding pure water energy parameters and the energy parameter of the object to be measured are obtained respectively according to the channels, then a ratio of the energy parameter of the object to be measured to the energy parameter of pure water is calculated; next, adjust the channel, repeat the extraction and calculation process, and finally obtain an energy parameter ratio corresponding to a series of channels; (2) calculate the ray path that the acoustic wave passes from the transmitting array element to the receiving array element; record the total number of the series of energy parameter ratios obtained in step (1) as n.sub.t, and the imaging area is divided so that the grid number of the imaging area satisfies Σ×Σ, wherein Σ is a positive integer, when √{square root over (.sub.t)} is an integer, Σ is equal to √{square root over (n.sub.t)}, when √{square root over (n.sub.t)} is a non-integer, Σ is an integer obtained by ceiling, flooring, or rounding off √{square root over (n.sub.t)}; the imaging area is established in the two-dimensional plane rectangular array element coordinate system corresponding to the transmitting array element and the receiving array element, then, for each group of transmitting array element-receiving array element, the straight line is connected between the coordinate position of the transmitting array element in the array element coordinate system and the coordinate position of the receiving array element in the array element coordinate system, thereby obtaining the path length of the connection line in each grid in the imaging area, and then obtaining the two-dimensional Σ×Σ matrix associated with the path length of each grid, thereafter, the matrix is vectorized to obtain the vector about the path length in each grid, finally, the vectors regarding the path length in each grid obtained by the each group of transmitting array element-receiving array element are arranged to form the two-dimensional Σ.sup.2×Σ.sup.2 matrix path matrix L regarding all of the transmitting array element-receiving array element groups; (3) solution for inverse problem: vectorize the series of energy parameter ratios obtained in step (1) to obtain ΔP; then construct a path-attenuation-energy parameter equation system as shown in equation (4):
LΔA=ΔP   (4) wherein, ΔA is an amount of change in attenuation coefficient to be solved; then, the quasi-Newton method is used to solve the equation system, and a one-dimensional vector ΔA containing Σ.sup.2 elements is obtained; thereafter, an attenuation coefficient of ultrasound wave in water is added to the amount of change in attenuation coefficient ΔA, and the attenuation coefficient reconstruction value vector of the object to be measured can be obtained.

3. The ultrasound CT image reconstruction method using the ultrasound CT sound speed reconstruction method based on ray theory according to claim 1, wherein the method utilizes the ultrasound CT sound speed reconstruction method based on ray theory according to claim 1 and further comprises the following steps: (4) imaging: two-dimensionalize the obtained speed reconstruction value vector of the object to be measured to form a Σ×Σ matrix; then obtain a two-dimensional pixel image based on the Σ×Σ matrix, and each pixel in the two-dimensional pixel image corresponds to a sound speed value.

4. The ultrasound CT image reconstruction method using the ultrasound CT sound speed reconstruction method based on ray theory according to claim 3, wherein in the step (4), the two-dimensional pixel image is obtained by performing logarithmic compression, grayscale mapping on the sound speed value, and finally displayed.

5. The ultrasound CT image reconstruction method using the ultrasound CT attenuation coefficient reconstruction method based on ray theory according to claim 2, wherein the method utilizes the ultrasound CT attenuation coefficient reconstruction method based on ray theory according to claim 2 and further comprises the following steps: (4) imaging: two-dimensionalize the obtained attenuation coefficient reconstruction value vector of the object to be measured to form a Σ×Σ matrix; then obtain a two-dimensional pixel image based on the Σ×Σ matrix, and each pixel in the two-dimensional pixel image corresponds to an attenuation coefficient value.

6. The ultrasound CT image reconstruction method using the ultrasound CT attenuation coefficient reconstruction method based on ray theory according to claim 5, characterized in that, in the step (4), the two-dimensional pixel image is obtained by performing logarithmic compression, grayscale mapping on the attenuation coefficient value, and finally displayed.

7. An ultrasound CT sound speed reconstruction system based on ray theory, wherein the system comprising: a travel time difference extraction module, configured to: based on the same transmitting array element and receiving array element, collect ultrasound transmission wave data from pure water and an object to be measured respectively after transmitting an ultrasound wave, and obtain respective corresponding pure water data and data of objects to be measured according to channels; utilize an AIC method to extract a travel time of the pure water data and then record the travel time as tof.sub.water; determine a matching window and record a starting point as tof.sub.water, record a time t.sub.water_max at the maximum amplitude of the pure water data as the end of window, then a window length is recorded as w, and a window data in the matching window is recorded as W.sub.water; find a sliding window whose window length is maintained at w on the data of the object to be measured in the corresponding channel, the window data in the sliding window is recorded as W.sub.object; W.sub.object and W.sub.water are correlated to each other to calculate a cross-correlation coefficient; a starting point of the sliding window is adjusted to slide the sliding window to obtain a series of cross-correlation coefficients, choose the sliding window corresponding to the cross-correlation coefficient with the largest value among the cross-correlation coefficients as a search result of the sliding window, and record the starting point of the search result of the sliding window as tof.sub.object, then the difference in travel time is Δtof=tof.sub.object−tof.sub.water; adjust the channel, repeat the extraction process, and finally obtain a series of travel time differences Δtof corresponding to a series of channels; a calculation module for calculating a ray path that acoustic waves pass from the transmitting array element to the receiving array element, configured to: record the total number of the series of effective travel time differences Δtof as n.sub.t, and an imaging area is divided so that the grid number of the imaging area satisfies Σ×Σ, wherein Σ is a positive integer, when √{square root over (n.sub.t)} is an integer, Σ is equal to √{square root over (n.sub.t)}; when √{square root over (n.sub.t)} is a non-integer, Σ is an integer obtained by ceiling, flooring, or rounding off √{square root over (n.sub.t)}, the imaging area is established in a two-dimensional plane rectangular array element coordinate system corresponding to the transmitting array element and the receiving array element, for each group of transmitting array element-receiving array element, a straight line is connected between a coordinate position of the transmitting array element in a array element coordinate system and a coordinate position of the receiving array element in the array element coordinate system, thereby obtaining a path length of the connection line in each grid in the imaging area, and then obtaining a two-dimensional Σ×Σ matrix associated with the path length of each grid, the matrix is vectorized to obtain a vector about the path length in each grid, vectors regarding the path length in each grid obtained by the each group of transmitting array element-receiving array element are arranged to form a two-dimensional Σ.sup.2×Σ.sup.2 matrix path matrix L regarding all of the transmitting array element-receiving array element groups; an inverse problem solution module, configured to vectorize the obtained effective travel time difference Δof to obtain ΔT; then construct a path-slowness-time equation system as shown in equation (3):
LΔS=ΔT   (3) wherein, ΔS is an amount of change in slowness to be solved; a quasi-Newton method is used to solve the equation system, and a one-dimensional vector ΔS containing Σ2 elements is obtained; thereafter, the slowness of ultrasound wave in water is added to an amount of change in slowness ΔS, take a reciprocal, and a speed reconstruction value vector of the object to be measured can be obtained.

8. An ultrasound CT attenuation coefficient reconstruction system based on ray theory, wherein the system comprising: an energy parameter extraction module, configured to, based on the same transmitting array element and receiving array element, energy parameters from pure water and the object to be measured are collected respectively after an ultrasound wave is transmitted, and the respective corresponding pure water energy parameters and the energy parameter of the object to be measured are obtained respectively according to channels, then a ratio of the energy parameter of the object to be measured to the energy parameter of pure water is calculated, adjust the channel, repeat the extraction and calculation process, and finally obtain an energy parameter ratio corresponding to a series of channels; a calculation module for calculating a ray path that the acoustic waves pass from the transmitting array element to the receiving array element, configured to: record the total number of a series of energy parameter ratios as n.sub.t, and an imaging area is divided so that a grid number of the imaging area satisfies Σ×Σ, wherein Σ is a positive integer, when √{square root over (n.sub.t)} is an integer, Σ is equal to √{square root over (n.sub.t)}, when √{square root over (n.sub.t)} is a non-integer, Σ is an integer obtained by ceiling, flooring, or rounding off √{square root over (n.sub.t)}, the imaging area is established in a two-dimensional plane rectangular array element coordinate system corresponding to the transmitting array element and the receiving array element, for each group of transmitting array element-receiving array element, a straight line is connected between a coordinate position of the transmitting array element in an array element coordinate system and a coordinate position of the receiving array element in the array element coordinate system, thereby obtaining a path length of the connection line in each grid in the imaging area, and then obtaining a two-dimensional Σ×Σ matrix associated with the path length of each grid, the matrix is vectorized to obtain a vector about the path length in each grid, vectors regarding the path length in each grid obtained by the each group of transmitting array element-receiving array element are arranged to form a two-dimensional Σ.sup.2×Σ.sup.2 matrix path matrix L regarding all of the transmitting array element-receiving array element groups; an inverse problem solution module, configured to vectorize the obtained series of energy parameter ratios to obtain ΔP; then construct a path-attenuation-energy parameter equation system as shown in equation (4):
LΔA=ΔP   (4) wherein, ΔA is an amount of change in attenuation coefficient to be solved; then, a quasi-Newton method is used to solve the equation system, and a one-dimensional vector ΔA containing Σ.sup.2 elements is obtained, thereafter, the attenuation coefficient of ultrasound wave in water is added to the amount of change in attenuation coefficient ΔA, and an attenuation coefficient reconstruction value vector of the object to be measured can be obtained.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0056] FIG. 1 is a schematic diagram of AIC method.

[0057] FIG. 2 is a schematic diagram of AIC-MAX-CC method, wherein the upper diagram corresponds to pure water data and the lower diagram corresponds to phantom data.

[0058] FIG. 3 is a schematic diagram of a straight path.

[0059] FIG. 4 is an ultrasound CT ring array (i.e., UCT ring array) and a breast custom phantom 1768-00 (CIRSINC, USA).

[0060] FIG. 5 is a reflection image reconstructed from a certain layer of data obtained by ultrasound CT scanning of the phantom.

[0061] FIG. 6 is the result of reconstructing the layer of data by the sound speed reconstruction algorithm provided by the disclosure (i.e., the sound speed image reconstructed by the algorithm of the disclosure).

[0062] FIG. 7 is an attenuation coefficient image obtained through reconstruction by a sound speed reconstruction algorithm in Embodiment 2 of the disclosure.

DESCRIPTION OF THE EMBODIMENTS

[0063] In order to make the purpose, technical solutions and advantages of the disclosure clearer, the disclosure will be further described in detail in combination with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the disclosure, and are not intended to limit the disclosure. In addition, the technical features involved in the various embodiments of the disclosure described below can be combined with each other as long as there is no conflict with each other.

Embodiment 1

[0064] The sound speed reconstruction method based on the ray theory in the disclosure is to properly process the transmission data collected by the ultrasound CT system and finally make a reconstruction to obtain the ultrasound sound speed image. The steps include extracting the travel time, calculating the ray path, resolving the inverse problem and displaying the sound speed images.

[0065] The collected data is derived from the ultrasound CT hardware system. The hardware system mainly includes an ultrasound probe, a transmitting circuit module, a receiving circuit module, a data acquisition module, and a computer system. The connection between them and the way of they work with each other can be directly derived from existing technologies, such as literature (Junjie Song, S. W. L. Z., A Prototype System for Ultrasound Computer Tomography with Ring Array, in International conference Biomedical Image and signal processing. 2017). The hardware modules such as the transmitting circuit module, the receiving circuit module, and the data acquisition module can be constructed by referring to the existing technology. The ultrasound probe may be a ring array probe, a linear array probe, an area array probe, etc. Taking the ring array probe as an example, the transmission data is collected by using the ring transmission mode. The transmission wave is mainly received by the array element opposite to the position of the transmitting array element. Therefore, the disclosure uses the data received by the array element opposite to the transmitting array element to reconstruct the sound speed. If other linear array probes and area array probes are used, they can be processed in a similar manner.

[0066] First, the AIC method is adopted to extract the travel time of pure water data, which is recorded as tof.sub.water.

[0067] The algorithm of the AIC method operates as follows: first, select a suitable window near the estimated travel time on the data, and the window length is N. Take the current retrieval point as the dividing point, divide the window into two segments, calculate the AIC value of the current retrieval point, and search the points in the window in sequence, and record the point with the smallest AIC value as the travel time point.

[0068] Specifically, reference for the initially estimated travel time point can be derived from the conventional operation in the related art. For example, the approximate travel time point can be calculated based on the theoretical sound speed of water, or the approximate location of travel time point can be observed by naked eyes on the waveform. The window length N is also selected based on the actual shape of the waveform and experiences, it can be increased or decreased as needed if the selection result does not meet the need. Reference for the calculation method of the AIC value can be directly derived from the related art.

[0069] Then determine the matching window, take tof.sub.water as the starting point, and take the time t.sub.water_max at the maximum amplitude of pure water data as the end of window; the window length is recorded as w, and the window data in the matching window is recorded as W.sub.water.

[0070] Then look for the sliding window with window length w on the phantom data of the corresponding channel, and the data of the sliding window is recorded as W.sub.object. W.sub.object and W.sub.water are correlated to calculate the cross-correlation coefficient, and the sliding point with the largest cross-correlation coefficient is selected and recorded as p (p is negative when the starting point of sliding window is before tof.sub.water, p is positive when the starting point of sliding window is after tof.sub.water), then Δtof=p, as shown in FIG. 2. That is to say, the window data length of the sliding window is maintained as w, and then the sliding is performed. When sliding one point, one correlation coefficient is calculated, and the sliding point corresponding to the largest correlation coefficient is selected and recorded as p.

[0071] Since there are multiple channels, the travel time difference of each channel needs to be solved separately.

[0072] The imaging area is then divided. The size of the grid needs to be determined according to the number of effective travel time differences, that is, the number of effective equations of the path-slowness-time equation system. Invalid travel time difference such as some missing channel signals is an invalid one.

[0073] In order to maintain the positive definiteness of the equation, the number of unknown numbers and the number of equations need to maintain consistent. After extracting the travel time difference, it is necessary to remove the wrong channel. Assuming that the number of valid travel time differences after removal is n.sub.t, then the grid number is √{square root over (n.sub.t)}*√{square root over (n.sub.t)}. If √{square root over (n.sub.t)} is a decimal, it can be rounded (any rounding rule such as ceiling, flooring or rounding off can be adopted).

[0074] Then calculate the ray path that the acoustic wave passes from the transmitting array element to the receiving array element. In the disclosure, on the premise that the propagation medium is relatively uniform, so the propagation path of sound wave is approximated by a straight line. Then, the focus of the problem is shifted to calculating an intersection point of two connection lines passing through the grid between two array element positions. After acquiring the intersection point, the distance between every two intersection points is calculated in sequence, and the distance is the length of the path in a single grid (if there is only one intersection point in a grid, the path length in the grid is 0). As shown in FIG. 3, it is a schematic diagram of a straight path, S is the transmitting array element, and R is the receiving array element. The √{square root over (n.sub.t)}*√{square root over (n.sub.t)} is vectorized (that is, the two-dimensional matrix √{square root over (n.sub.t)}*√{square root over (n.sub.t)} is turned into a one-dimensional column vector). After vectorization, the blank grid is zero, representing the grid is not crossed by the grid, the gray grid is a non-zero value. The non-zero value is the path length in the network, which means that the grid is crossed by the path. Then the size of the path matrix L of n.sub.t paths is n.sub.t×n.sub.t.

[0075] For the above grid, according to the dimensional parameters given by the probe manufacturer, the processing methods in the related art can be utilized to correspond each array element to a coordinate system, and the corresponding array element coordinates can be obtained.

[0076] After the path calculation is completed, the path-slowness-time equation system is constructed as shown in the following equation (3), wherein ΔS is the amount of change in slowness, ΔT is the travel time difference, and L is the n.sub.t×n.sub.t path matrix; the dimension of the equation (3) is the number of effective travel time difference.

LΔS=ΔT   (3)

[0077] Use the quasi-Newton method to solve the equation system, and output ΔS (the obtained ΔS is a vector, which can be a positive value or negative value). The slowness of water plus this amount of change in slowness, take the reciprocal, then the speed reconstruction value of the phantom can be obtained. The obtained speed reconstruction value of the phantom is also in the form of vector.

[0078] Imaging step: Two-dimensionalize the speed reconstruction value of the phantom in the form of a vector (that is, the inverse process of the above √{square root over (n.sub.t)}*√{square root over (n.sub.t)} matrix vectorization) and then obtain a two-dimensional pixel image, each pixel in the two-dimensional pixel image is a sound speed value.

Embodiment 2

[0079] Similar to Embodiment 1, the attenuation coefficient reconstruction method based on the ray theory in the disclosure includes:

[0080] The attenuation coefficient reconstruction method based on the ray theory in the disclosure is to properly process the transmission data collected by the ultrasound CT system, and finally make a reconstruction to obtain the ultrasound attenuation image. The steps include extracting signal energy parameters, calculating the ray path, solving the inverse problem, and displaying attenuation coefficient image.

[0081] The energy parameters may be the amplitude, intensity, and energy of the signal.

[0082] The calculation of the ray path is the same as in Embodiment 1.

[0083] The solution to the inverse problem is the same as in Embodiment 1.

[0084] After the path calculation is completed, the path-attenuation-energy parameter equation system is constructed, as shown in the following equation (4), wherein ΔA is the amount of change in the attenuation coefficient and ΔP is the energy parameter ratio, that is, the energy parameter ratio of the energy parameter of phantom data to the energy parameter of the water data.

LΔA=ΔP   (4)

[0085] Adopt the quasi-Newton method to solve (4), then ΔA is obtained. The attenuation coefficient of water plus the amount of change in this attenuation coefficient, then the reconstruction value of the attenuation coefficient of the phantom can be obtained.

[0086] Finally, an image display of the attenuation coefficient is shown in FIG. 7.

[0087] In the above embodiment, for the obtained ultrasound image, logarithmic compression and grayscale mapping may be sequentially performed according to the method in the related art to obtain an ultrasound image with a grayscale value within a specified range.

[0088] The disclosure is applicable to various commercial ultrasound CT system probes, such as ring probes, linear array probes, area array probes, concave arrays, etc.

[0089] Those skilled in the art can easily understand that the above are only preferred embodiments of the disclosure and are not intended to limit the disclosure. Any modification, equivalent replacement and improvement made within the spirit and principle of the disclosure should fall within the scope of the disclosure.