WO2020206755A1 - 一种基于射线理论的超声ct图像重建方法及系统 - Google Patents

一种基于射线理论的超声ct图像重建方法及系统 Download PDF

Info

Publication number
WO2020206755A1
WO2020206755A1 PCT/CN2019/084712 CN2019084712W WO2020206755A1 WO 2020206755 A1 WO2020206755 A1 WO 2020206755A1 CN 2019084712 W CN2019084712 W CN 2019084712W WO 2020206755 A1 WO2020206755 A1 WO 2020206755A1
Authority
WO
WIPO (PCT)
Prior art keywords
array element
transmitting
window
path
grid
Prior art date
Application number
PCT/CN2019/084712
Other languages
English (en)
French (fr)
Inventor
尉迟明
丁明跃
方小悦
武云
宋俊杰
周亮
张求德
刘阔林
Original Assignee
华中科技大学
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 华中科技大学 filed Critical 华中科技大学
Priority to US16/968,151 priority Critical patent/US20210236095A1/en
Publication of WO2020206755A1 publication Critical patent/WO2020206755A1/zh

Links

Images

Classifications

    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/13—Tomography
    • A61B8/15—Transmission-tomography
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S15/00—Systems using the reflection or reradiation of acoustic waves, e.g. sonar systems
    • G01S15/88—Sonar systems specially adapted for specific applications
    • G01S15/89—Sonar systems specially adapted for specific applications for mapping or imaging
    • G01S15/8906—Short-range imaging systems; Acoustic microscope systems using pulse-echo techniques
    • G01S15/8977—Short-range imaging systems; Acoustic microscope systems using pulse-echo techniques using special techniques for image reconstruction, e.g. FFT, geometrical transformations, spatial deconvolution, time deconvolution
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/13—Tomography
    • A61B8/14—Echo-tomography
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/52—Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/5207—Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves involving processing of raw data to produce diagnostic data, e.g. for generating an image
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/52—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S15/00
    • G01S7/52017—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S15/00 particularly adapted to short-range imaging
    • G01S7/52023—Details of receivers
    • G01S7/52036—Details of receivers using analysis of echo signal for target characterisation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/00—2D [Two Dimensional] image generation
    • G06T11/003—Reconstruction from projections, e.g. tomography
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/00—2D [Two Dimensional] image generation
    • G06T11/003—Reconstruction from projections, e.g. tomography
    • G06T11/006—Inverse problem, transformation from projection-space into object-space, e.g. transform methods, back-projection, algebraic methods
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/08—Detecting organic movements or changes, e.g. tumours, cysts, swellings
    • A61B8/0825—Detecting organic movements or changes, e.g. tumours, cysts, swellings for diagnosis of the breast, e.g. mammography

Definitions

  • the invention relates to a transmissive ultrasound imaging mode in ultrasound tomography, belonging to the technical field of functional imaging, and more specifically, to a method and system for reconstructing ultrasound CT images based on radiation theory.
  • Ultrasonic CT refers to transmitting ultrasound waves to an object through an ultrasound probe and receiving reflection data or transmission data, and using these data to reconstruct ultrasound tomographic images in order to observe the three-dimensional information inside the object.
  • Ultrasonic testing has the advantages of low price and harmless to human body. With the rapid development of probe processing technology and computer high-performance computing, ultrasonic tomography technology has once again become a research hotspot in recent years.
  • the reflection image of ultrasound CT has a higher image resolution, so as to assist the doctor to see the smaller diseased tissue.
  • the transmission data can be used to reconstruct sound speed, attenuation coefficient and other functional parameters, which belong to the field of functional imaging. Studies have pointed out that in the early stage of the disease, the diseased tissue first changes in functional parameters and then has structural changes. Therefore, transmission imaging is of great significance to the early imaging and diagnosis of lesions, and it can diagnose lesions earlier.
  • the ultrasonic CT sound velocity reconstruction method includes the reconstruction method based on the ray theory and the reconstruction method based on the wave theory.
  • the reconstruction method based on the wave theory has higher imaging resolution, but is greatly affected by the small error disturbance, so it is not stable enough, and the amount of calculation is huge, and it is not suitable for practical applications at present.
  • the reconstruction method based on ray theory has a simpler model, higher robustness, and a smaller amount of calculation. At present, it is a kind of efficient, stable and practical sound velocity reconstruction method.
  • the purpose of the present invention is to provide a method and system for reconstructing ultrasound CT images based on radiation theory, in which the overall process of the method is improved, especially by optimizing the radiation theory, Using specific ray calculation processing methods, fast and stable ultrasonic CT sound velocity reconstruction and ultrasonic CT attenuation coefficient reconstruction can be realized, and then ultrasonic CT image reconstruction based on ray theory can be realized.
  • a method for reconstruction of ultrasonic CT sound velocity based on ray theory which is characterized in that it comprises the following steps:
  • the data of the ultrasonic transmission wave from the pure water and the object to be measured are collected respectively, and the pure water data and the data of the object to be measured corresponding to each channel are obtained respectively;
  • search for a sliding window with a window length of w on the target data of the corresponding channel and record the window data in the sliding window as W object ; then, correlate W object and W water to calculate the correlation coefficient; adjust The window starting point of the sliding window, thereby sliding the sliding window to obtain a series of cross-correlation coefficients, selecting the sliding window corresponding to the cross-correlation coefficient with the largest value among these cross-correlation coefficients as the sliding window search result, and recording the value of the sliding window search result
  • the total number of the difference ⁇ tof between a series of effective transit times obtained in the step (1) is n t , and the imaging area is divided so that the grid number of the imaging area meets ⁇ ; where ⁇ is a positive integer, when When it is an integer, then ⁇ is equal to when When it is not an integer, ⁇ is Press round up, round down or round to the whole number;
  • the imaging area is established in the two-dimensional planar rectangular array element coordinate system corresponding to the transmitting array element and the receiving array element, and then for each group of transmitting array element-receiving array element group, from the transmitting array element in the array element coordinate system Draw a straight line from the coordinate position of the receiving array element to the coordinate position of the array element coordinate system to obtain the path length of the line in each grid in the imaging area, and then obtain the path length of each grid , ⁇ two-dimensional matrix, and then vectorize the matrix to obtain a vector about the path length in each grid; finally, combine each group of transmitting array element-receiving array element group to obtain information about the path in each grid
  • the length vector arrangement forms a ⁇ 2 ⁇ ⁇ 2 two-dimensional matrix path matrix L of all transmitting array elements-receiving array element groups;
  • ⁇ S is the slowness change to be solved; then, the quasi-Newton method is used to solve the equations to obtain a one-dimensional vector ⁇ S containing ⁇ 2 elements; then, the slowness of the ultrasonic wave in the water plus the slowness change ⁇ S And take the reciprocal to get the velocity reconstruction value vector of the object to be measured.
  • the present invention provides a method for reconstruction of ultrasonic CT attenuation coefficients based on ray theory, which is characterized in that it comprises the following steps:
  • the energy parameters of the ultrasonic transmitted waves from pure water and the object to be measured are collected after the ultrasonic wave is transmitted, and the pure water energy parameters corresponding to each channel and the object to be measured are obtained respectively.
  • Object energy parameters and then calculate the ratio between the energy parameters of the object to be measured and the pure water energy parameters;
  • the total number of ratios of a series of energy parameters obtained in step (1) is n t , and divide the imaging area so that the grid number of the imaging area satisfies ⁇ ; where ⁇ is a positive integer, when When it is an integer, then ⁇ is equal to when When it is not an integer, ⁇ is Press round up, round down or round to the whole number;
  • the imaging area is established in the two-dimensional planar rectangular array element coordinate system corresponding to the transmitting array element and the receiving array element, and then for each group of transmitting array element-receiving array element group, from the transmitting array element in the array element coordinate system Draw a straight line from the coordinate position of the receiving array element to the coordinate position of the array element coordinate system to obtain the path length of the line in each grid in the imaging area, and then obtain the path length of each grid , ⁇ two-dimensional matrix, and then vectorize the matrix to obtain a vector about the path length in each grid; finally, combine each group of transmitting array element-receiving array element group to obtain information about the path in each grid
  • the length vector arrangement forms a ⁇ 2 ⁇ ⁇ 2 two-dimensional matrix path matrix L of all transmitting array elements-receiving array element groups;
  • ⁇ A is the change of the attenuation coefficient to be solved; then, the quasi-Newton method is used to solve the equations to obtain a one-dimensional vector ⁇ A containing ⁇ 2 elements; then, the attenuation coefficient of the ultrasonic wave in the water plus the change of the attenuation coefficient ⁇ A can be used to obtain the reconstructed value vector of the attenuation coefficient of the object under test.
  • the present invention provides an ultrasound CT image reconstruction method using the above-mentioned ray theory-based ultrasound CT sound velocity reconstruction method, characterized in that the method uses the ray theory-based ultrasound as claimed in claim 1.
  • the CT sound velocity reconstruction method also includes the following steps:
  • the obtained velocity reconstruction value vector of the object to be measured is two-dimensionalized to form a ⁇ matrix; then a two-dimensional pixel image is obtained based on the ⁇ matrix, and each pixel in the two-dimensional pixel image is The sound velocity value corresponds to it.
  • the two-dimensional pixel map is obtained by performing logarithmic compression on the sound velocity value, gray-scale mapping, and finally displaying.
  • the present invention provides an ultrasound CT image reconstruction method using the ray theory-based ultrasound CT attenuation coefficient reconstruction method as claimed in claim 2, characterized in that the method uses the method as claimed in claim 2.
  • the attenuation coefficient reconstruction method of ultrasonic CT based on ray theory also includes the following steps:
  • the obtained attenuation coefficient reconstruction value vector of the object to be tested is two-dimensionalized to form a ⁇ matrix; then a two-dimensional pixel image is obtained based on the ⁇ matrix, and each pixel in the two-dimensional pixel image Corresponds to the attenuation coefficient value.
  • the two-dimensional pixel map is obtained by performing logarithmic compression on the attenuation coefficient value, grayscale mapping, and finally displaying.
  • the present invention provides an ultrasonic CT sound velocity reconstruction system based on ray theory, which is characterized in that the system includes:
  • Transit time difference extraction module used for: Based on the same transmitting and receiving array element, after transmitting the ultrasonic wave, the data of the ultrasonic transmission wave from the pure water and the object to be measured are collected respectively, and the pure water corresponding to each channel is obtained.
  • Data, and the data of the object to be measured use the AIC method to extract the transit time of the pure water data, and record it as tof water ; determine the matching window, take tof water as the starting point of the window, and use the time t water_max at the maximum amplitude of the pure water data as the window
  • the window length is denoted as w
  • the window data in the matching window is denoted as W water
  • the sliding window with the window length maintained as w is found on the object data of the corresponding channel, and the window data in the sliding window is denoted as W object ; cross-correlate W object and W water to calculate the cross-correlation coefficient; adjust the window starting point of the sliding window, thereby sliding the sliding window to obtain a series of cross-correlation coefficients, and select the cross-correlation coefficient with the largest value among these cross-correlation coefficients
  • the sliding window corresponding to the number is used as the sliding window search result.
  • the ray path calculation module of the sound wave from the transmitting element to the receiving element is used to: remember the total number of the difference ⁇ tof of a series of effective transit time is n t , and divide the imaging area to make the network of the imaging area
  • the grid number satisfies ⁇ ; where ⁇ is a positive integer, When it is an integer, then ⁇ is equal to when When it is not an integer, ⁇ is An integer obtained by rounding up, rounding down or rounding up; 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.
  • each group of transmitting array elements -Receiving array element group For each group of transmitting array elements -Receiving array element group, from the coordinate position of the transmitting array element in the array element coordinate system to the coordinate position of the receiving array element in the array element coordinate system to make a straight line to obtain each grid of the line in the imaging area The length of the path in each grid is obtained, and a ⁇ two-dimensional matrix of the path length in each grid is obtained, and the matrix is vectorized to obtain a vector of the path length in each grid; each group of transmitting elements-receiving The vector arrangement of the path length in each grid obtained by the array element group forms a ⁇ 2 ⁇ ⁇ 2 two-dimensional matrix path matrix L of ⁇ 2 ⁇ ⁇ 2 for all the transmitting array element-receiving array element group;
  • the inverse problem solving module is used to: vectorize the obtained effective transit time difference ⁇ tof to obtain ⁇ T; construct a path-slowness-time equation set as shown in formula (3):
  • ⁇ S is the slowness change to be solved
  • the quasi-Newton method is used to solve the equations to obtain a one-dimensional vector ⁇ S containing ⁇ 2 elements; the slowness of the ultrasonic wave in the water plus the slowness change ⁇ S, take the reciprocal, and obtain the velocity reconstruction value vector of the object to be measured.
  • the present invention provides an ultrasonic CT attenuation coefficient reconstruction system based on ray theory, which is characterized in that the system includes:
  • the energy parameter extraction module is used to collect the energy parameters of the ultrasonic transmitted waves from pure water and the object to be measured after transmitting the ultrasound based on the same transmitting and receiving array element, and obtain the pure water energy parameters corresponding to each channel one by one. , And the energy parameters of the object to be measured, calculate the ratio between the energy parameter of the object to be measured and the energy parameter of pure water; adjust the channel, repeat the extraction and calculation processing, and finally obtain the energy parameter ratio corresponding to a series of channels;
  • the ray path calculation module of the acoustic wave from the transmitting element to the receiving element is used to: the total number of the memorized series of energy parameter ratios is n t , and the imaging area is divided so that the grid number of the imaging area meets ⁇ ⁇ ; where ⁇ is a positive integer, when When it is an integer, then ⁇ is equal to when When it is not an integer, ⁇ is An integer obtained by rounding up, rounding down or rounding up; 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.
  • each group of transmitting array elements -Receiving array element group For each group of transmitting array elements -Receiving array element group, from the coordinate position of the transmitting array element in the array element coordinate system to the coordinate position of the receiving array element in the array element coordinate system to make a straight line to obtain each grid of the line in the imaging area The length of the path in each grid is obtained, and a ⁇ two-dimensional matrix of the path length in each grid is obtained, and the matrix is vectorized to obtain a vector of the path length in each grid; each group of transmitting elements-receiving The vector arrangement of the path length in each grid obtained by the array element group forms a ⁇ 2 ⁇ ⁇ 2 two-dimensional matrix path matrix L of ⁇ 2 ⁇ ⁇ 2 for all the transmitting array element-receiving array element group;
  • the inverse problem solving module is used to: vectorize the obtained series of energy parameter ratios to obtain ⁇ P; construct the path-attenuation-energy parameter equation set as shown in formula (4):
  • ⁇ A is the change of attenuation coefficient to be solved
  • the quasi-Newton method is used to solve the equations to obtain a one-dimensional vector ⁇ A containing ⁇ 2 elements; the attenuation coefficient of the ultrasonic wave in the water is added to the attenuation coefficient change ⁇ A to obtain the attenuation coefficient reconstruction value vector of the object to be measured.
  • the above technical solution conceived by the present invention is based on ray theory and uses parameters such as the difference in transit time ( ⁇ tof) or the ratio of energy parameters in the reconstruction process, which can achieve fast and stable Ultrasonic CT sound velocity reconstruction and ultrasonic CT attenuation coefficient reconstruction.
  • ⁇ tof difference in transit time
  • the so-called ray-based theory assumes that the propagation medium is relatively uniform, and the propagation path of ultrasound in a homogeneous medium can be approximated as a ray.
  • the main steps include the extraction of transit time, the calculation of ray paths, and the solution of inverse problems:
  • the transit time is the time when the waveform arrives at the receiving position. This time is the time at which the waveform first arrives in the received waveform.
  • the difference in transit time refers to the difference in transit time between phantom data and pure water data.
  • the present invention adopts the method of combining the cross-correlation method (CC) and the mutual information method (AIC), and introduces the maximum amplitude position information (MAX) to extract the difference ⁇ tof between the phantom data and the pure water data, referred to as AIC- MAX-CC method.
  • the present invention is based on the premise that the propagation medium is relatively uniform, so the path is approximated by a straight path.
  • the problem is converted to the problem of calculating the intersection of two points through the grid. After the intersection is obtained, the distance between the two intersections is sequentially calculated, which is the length of the path in a single grid.
  • the last is to establish and solve the path-slowness-time equations, that is, the establishment and solution of the inverse problem.
  • Slowness is the reciprocal of the speed of sound.
  • the slowness in the present invention refers to the increase in slowness of the phantom relative to pure water.
  • the path is the length of the propagation path across the grid, and the time is the transit time of the phantom and pure water data. Difference.
  • the present invention adopts the quasi-Newton method to solve the equations.
  • the quasi-Newton method is between the gradient descent method and the Newton method, and is a good iterative optimization method.
  • the attenuation coefficient reconstruction method based on ray theory, in addition to replacing the difference in transit time ( ⁇ tof) with the ratio of energy parameters, and replacing the path-slowness-time equation set with the path-attenuation-energy parameter equation set, is different from the ray-based
  • the theoretical sound velocity reconstruction method is roughly similar.
  • the advantages of the sound velocity reconstruction method and system based on the ray theory proposed by the present invention are as follows: 1.
  • the construction of the equations adopts the difference in transit time instead of the transit time itself, which reduces the system error Influence; 2.
  • the extraction of the difference in transit time combines the cross-correlation method and the mutual information method, and introduces the maximum amplitude position information, referred to as the AIC-MAX-CC method, which has strong anti-noise ability; 3.
  • Simple path calculation Based on the assumption that the sound velocity of the propagation medium is uniform, the path calculation is transformed into the problem of calculating the intersection of two points through the grid, which is simple and efficient; 4.
  • the inverse problem is solved by the quasi-Newton method, which avoids the calculation of the Hessian matrix, and the calculation is small. High precision.
  • the advantages of the attenuation coefficient reconstruction method and system based on the ray theory are as follows: 1. The construction of the equation system uses the ratio of energy parameters instead of the energy itself, which reduces the influence of system errors; 2. The path calculation is simple, Based on the assumption that the sound velocity of the propagation medium is uniform, the path calculation is transformed into the problem of calculating the intersection of two points through the grid, which is simple and efficient; 3. The inverse problem is solved by the quasi-Newton method, which avoids the calculation of the Hessian matrix, and the calculation is small. High precision.
  • Figure 1 is a schematic diagram of the AIC method.
  • Figure 2 is a schematic diagram of the AIC-MAX-CC method, where the upper figure corresponds to pure water data, and the lower figure corresponds to phantom data.
  • Figure 3 is a schematic diagram of a straight path.
  • Figure 4 shows the ultrasound CT ring array (ie UCT ring array) and breast customized phantom 1768-00 (CIRSINC, USA).
  • Fig. 5 is a reflection image reconstructed from a certain layer of data obtained by ultrasound CT scanning of the phantom.
  • Fig. 6 is the result of reconstruction of the layer data by the sound velocity reconstruction algorithm proposed by the present invention (ie, the sound velocity image reconstructed by the algorithm of the present invention).
  • Fig. 7 is an attenuation coefficient image reconstructed by the sound velocity reconstruction algorithm in Embodiment 2 of the present invention.
  • the sound velocity reconstruction method based on the ray theory in the present invention is to appropriately process the transmission data collected by the ultrasonic CT system, and finally reconstruct the ultrasonic sound velocity image.
  • the steps include the extraction of the transit time, the calculation of the ray path, and the solution of the inverse problem. , The display of sound speed images.
  • the collected data comes 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 relationship between them and the working method of cooperation can directly refer to the existing technology, such as the 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 with reference to the existing technology.
  • the ultrasonic probe can be a ring array probe, a linear array probe, an area array probe, and the like.
  • the circular array probe Take the circular array probe as an example, use the mode of circular emission to collect transmission data.
  • the transmitted wave is mainly received by the element opposite to the position of the transmitting element, so the present invention uses the data received by the element opposite to the transmitting element to reconstruct the sound velocity. If other linear array probes and area array probes are used, similar treatments can be applied.
  • the first is to use the AIC method to extract the transit time of pure water data, which is recorded as tof water .
  • the algorithm flow of the AIC method is as follows: First, select an appropriate window near the estimated transit time on the data, and the window length is N. Taking the current retrieval point as the demarcation point, divide the window into two sections, calculate the AIC value of the current retrieval point, search the points in the window in turn, record the point with the smallest AIC value as the transit time point.
  • the initial estimated transit time point can refer to conventional operations in the prior art.
  • the approximate transit time point can be calculated according to the theoretical sound velocity of water, or the transit time point can be observed on the waveform with the naked eye.
  • the approximate location The selection of the window length N is also based on the actual shape of the waveform and empirically selected. If it is not suitable, it can be appropriately increased or decreased.
  • the calculation method of the AIC value can directly refer to the related prior art.
  • the matching window is determined, taking tof water as the starting point of the window, taking the time t water_max at the maximum amplitude of the pure water data as the end of the window, the window length is denoted as w, and the matching window data is denoted as W water .
  • the imaging area is then divided.
  • the size of the grid needs to be determined according to the number of effective transit time differences, that is, the number of effective equations in the path-slowness-time equations. Invalid transit time difference, for example, some channel signal is missing, it is invalid.
  • the number of unknowns and the number of equations need to be consistent. After extracting the transit time difference, the wrong channel should be eliminated. Assuming that the effective transit time difference after elimination is n t , the number of grids is If If it is a decimal, then round up (round up, round down, or rounding up to the nearest integer can be used).
  • the present invention is based on the premise that the propagation medium is relatively uniform, so the propagation path of the sound wave can be replaced by a straight line approximation.
  • the problem is transformed into the problem of obtaining the intersection point of the line of the two array elements passing through the grid. After the intersection is obtained, the distance between the two intersections is sequentially calculated, and the distance is the length of the path in a single grid (if there is only one intersection in a grid, the path length in the grid is 0).
  • Figure 3 it is a schematic diagram of a straight path, S is the transmitting array element, and R is the receiving array element.
  • Matrix vectorization that is, the The two-dimensional matrix becomes a one-dimensional column vector.
  • the blank grid has a zero value, which means that the grid is not crossed by a path.
  • the gray grid is a non-zero value, and the non-zero value is the length of the path in the network. It means that the grid is crossed by the path.
  • the path n t paths matrix L size is n t ⁇ n t.
  • the processing methods in the prior art can be used to correspond each array element to a coordinate system, and the corresponding array element coordinates can be obtained.
  • ⁇ S is the amount of slowness change
  • ⁇ T is the difference in transit time
  • L is the path of n t ⁇ n t Matrix
  • the dimension of the equation system (3) is the number of effective transit time differences.
  • the quasi-Newton method uses the quasi-Newton method to solve the equations, output ⁇ S (the obtained ⁇ S is a vector, and the value can be positive or negative), the slowness of water plus this slowness change, and the reciprocal is the speed reconstruction value of the phantom .
  • the obtained velocity reconstruction value of the phantom is also in vector form.
  • Imaging steps Reconstruct the velocity value of the phantom in vector form into two dimensions (that is, the above After the inverse process of matrix vectorization), a two-dimensional pixel map is obtained, and each pixel in the two-dimensional pixel map is the sound velocity value.
  • the attenuation coefficient reconstruction method based on ray theory in the present invention includes:
  • the attenuation coefficient reconstruction method based on the ray theory in the present invention is to appropriately process the transmission data collected by the ultrasound CT system, and finally reconstruct the ultrasound attenuation image.
  • the steps include the extraction of signal energy parameters, the calculation of the ray path, and the inverse problem. Solve, display the attenuation coefficient image.
  • the energy parameter can be the amplitude, intensity, and energy of the signal.
  • the path-attenuation-energy parameter equation set is constructed, as shown in the following formula (4), where ⁇ A is the variation of the attenuation coefficient, and ⁇ P is the ratio of the energy parameter, that is, the energy parameter of the phantom data and the water data The ratio of the energy parameters.
  • logarithmic compression and gray scale mapping may be performed in sequence according to the method in the prior art to obtain an ultrasound image with a gray value within a specified range.
  • the invention is suitable for various commercial ultrasonic CT system probes, such as ring probes, linear array probes, area array probes, concave arrays, etc.

Landscapes

  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Molecular Biology (AREA)
  • Veterinary Medicine (AREA)
  • Public Health (AREA)
  • General Health & Medical Sciences (AREA)
  • Animal Behavior & Ethology (AREA)
  • Surgery (AREA)
  • Medical Informatics (AREA)
  • Biomedical Technology (AREA)
  • Radiology & Medical Imaging (AREA)
  • Pathology (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Biophysics (AREA)
  • General Physics & Mathematics (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Acoustics & Sound (AREA)
  • Mathematical Optimization (AREA)
  • Algebra (AREA)
  • Mathematical Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Ultra Sonic Daignosis Equipment (AREA)

Abstract

一种基于射线理论的超声CT图像重建方法及系统，方法包括基于射线理论的超声CT声速重建方法及超声CT衰减系数重建方法；基于射线理论的超声CT声速重建方法包括：(1)渡越时间之差的提取；(2)计算声波从发射阵元到接收阵元经过的射线路径；(3)反问题的求解：通过采用拟牛顿法解路径-慢度-时间方程组，即可得到待测对象的速度重建值向量。通过对方法的整体流程进行改进，尤其是通过对射线理论的优化，利用特定的射线计算处理方式，可实现快速、稳定的超声CT声速重建及超声CT衰减系数重建，进而实现了基于射线理论的超声CT图像重建。

Description

一种基于射线理论的超声CT图像重建方法及系统 【技术领域】
本发明涉及超声断层成像中透射式超声成像模式，属于功能成像技术领域，更具体地，涉及一种基于射线理论的超声CT图像重建方法及系统。
【背景技术】
超声CT是指通过超声探头对物体发射超声波并接收反射数据或者透射数据，利用这些数据重建出超声断层图像，以便观测物体内部的三维信息。超声检测具有价格低廉、对人体无害等优点，随着探头加工工艺和计算机高性能运算的快速发展，超声断层成像技术近些年来又再次成为了研究热点。
超声CT成像有两种成像模式，反射式成像和透射式成像。由于采集多方位的反射信息，超声CT的反射图像具有更高的图像分辨率，以便辅助医生看到更微小的病变组织。而通过透射数据可以重建出声速、衰减系数等功能参数，属于功能像的领域。有研究指出，在病变早期，病变组织先有功能参数的变化，后有结构变化。因此，透射式成像对病变的早期成像和诊断具有重要意义，能更早地诊断出病变。
超声CT声速重建和衰减系数重建的重建方法大致相同，以超声CT声速重建方法为例，超声CT声速重建方法包括基于射线理论的重建方法和基于波动理论的重建方法。基于波动理论的重建方法成像分辨率更高，但是受微小误差扰动影响较大，因而不够稳定，并且运算量极大，目前还不适用于实际应用。基于射线理论的重建方法模型更简单，鲁棒性更高，并且运算量更小，目前来看是一类高效的、稳定的、实用的声速重建方法。
目前国内外对超声CT声速重建方法的研究甚少。主要原因是该技术存在几个难点：超声CT系统的搭建难度较大，数据获取困难；超声CT声速 重建方法涉及大尺度的矩阵运算，需要大量的计算开销；超声CT声速重建过程需要解决一个反问题，因而反问题的求解是一个难点。
基于射线理论的超声CT声速重建方法存在几个难点：渡越时间的准确提取难度较大，受噪声影响和系统误差影响较大；射线路径的计算方法众多，针对不同应用场景适用性不同；同样需要解决一个反问题，反问题的求解是一个难点。
【发明内容】
针对现有技术的以上缺陷或改进需求，本发明的目的在于提供一种基于射线理论的超声CT图像重建方法及系统，其中通过对方法的整体流程进行改进，尤其是通过对射线理论的优化，利用特定的射线计算处理方式，可实现快速、稳定的超声CT声速重建及超声CT衰减系数重建，进而实现了基于射线理论的超声CT图像重建。
为实现上述目的，按照本发明的一个方面，提供了一种基于射线理论的超声CT声速重建方法，其特征在于，包括以下步骤：
(1)渡越时间之差的提取：
首先基于同一发射阵元和接收阵元，发射超声波后分别采集来自纯水和待测对象的超声透射波的数据，分别得到按通道一一对应的纯水数据、以及待测对象数据；
然后采用AIC法提取纯水数据的渡越时间，记为tof water；接着，确定匹配窗，以tof water为窗起点，以纯水数据的最大幅值处时间t water_max为窗终点，则，窗长度记为w，匹配窗内的窗数据记为W water；
随后在对应通道的待测对象数据上寻找窗长度保持为w的滑动窗，该滑动窗内的窗数据记为W object；接着，将W object和W water互相关，计算得到互相关系数；调整所述滑动窗的窗起点，从而滑动所述滑动窗得到一系列互相关系数，选取这些互相关系数中数值最大的互相关系数对应的滑动窗作 为滑动窗寻找结果，记该滑动窗寻找结果的窗起点为tof object，则渡越时间之差Δtof＝tof object-tof water；
接着，调整通道，重复提取处理，最终得到一系列通道对应的一系列渡越时间之差Δtof；
(2)计算声波从发射阵元到接收阵元经过的射线路径：
记所述步骤(1)得到的一系列有效渡越时间之差Δtof的总数量为n t，对成像区域进行剖分，使成像区域的网格数满足Σ×Σ；其中Σ为正整数，当
Figure PCTCN2019084712-appb-000001
为整数时，则Σ等于
Figure PCTCN2019084712-appb-000002
当
Figure PCTCN2019084712-appb-000003
为非整数时，Σ为
Figure PCTCN2019084712-appb-000004
按向上取整、向下取整或四舍五入取整得到的整数；
所述成像区域是建立在发射阵元和接收阵元对应的二维平面直角阵元坐标系中，然后对于每一组发射阵元-接收阵元组，从发射阵元在阵元坐标系中的坐标位置向接收阵元在阵元坐标系中的坐标位置作直线连线，得到该连线在该成像区域中每一个网格中的路径长度，进而得到关于每个网格中路径长度的、Σ×Σ的二维矩阵，接着将该矩阵向量化得到关于每个网格中路径长度的向量；最后，将每一组发射阵元-接收阵元组得到的关于每个网格中路径长度的向量排列形成关于全部发射阵元-接收阵元组的Σ 2×Σ 2的二维矩阵路径矩阵L；
(3)反问题的求解：
将所述步骤(1)得到的有效渡越时间之差Δtof进行向量化，得到ΔT；接着构建如式(3)所示的路径-慢度-时间方程组：
LΔS＝ΔT   (3)
其中，ΔS为待求解的慢度变化量；然后，采用拟牛顿法解方程组，求解得到一维含有Σ 2个元素的向量ΔS；接着，用水中超声波的慢度加上慢度变化量ΔS，取倒数，即可得到待测对象的速度重建值向量。
按照本发明的另一方面，本发明提供了一种基于射线理论的超声CT衰 减系数重建方法，其特征在于，包括以下步骤：
(1)首先基于同一发射阵元和接收阵元，发射超声波后分别采集来自纯水和待测对象的超声透射波的能量参数，分别得到按通道一一对应的纯水能量参数、以及待测对象能量参数，然后计算得到待测对象能量参数与纯水能量参数之间的比值；
接着，调整通道，重复提取与计算处理，最终得到一系列通道对应的能量参数比值；
(2)计算声波从发射阵元到接收阵元经过的射线路径：
记所述步骤(1)得到的一系列能量参数比值的总数量为n t，对成像区域进行剖分，使成像区域的网格数满足Σ×Σ；其中Σ为正整数，当
Figure PCTCN2019084712-appb-000005
为整数时，则Σ等于
Figure PCTCN2019084712-appb-000006
当
Figure PCTCN2019084712-appb-000007
为非整数时，Σ为
Figure PCTCN2019084712-appb-000008
按向上取整、向下取整或四舍五入取整得到的整数；
所述成像区域是建立在发射阵元和接收阵元对应的二维平面直角阵元坐标系中，然后对于每一组发射阵元-接收阵元组，从发射阵元在阵元坐标系中的坐标位置向接收阵元在阵元坐标系中的坐标位置作直线连线，得到该连线在该成像区域中每一个网格中的路径长度，进而得到关于每个网格中路径长度的、Σ×Σ的二维矩阵，接着将该矩阵向量化得到关于每个网格中路径长度的向量；最后，将每一组发射阵元-接收阵元组得到的关于每个网格中路径长度的向量排列形成关于全部发射阵元-接收阵元组的Σ 2×Σ 2的二维矩阵路径矩阵L；
(3)反问题的求解：
将所述步骤(1)得到的一系列能量参数比值进行向量化，得到ΔP；接着构建如式(4)所示的路径-衰减-能量参数方程组：
LΔA＝ΔP   (4)
其中，ΔA为待求解的衰减系数变化量；然后，采用拟牛顿法解方程组，求 解得到一维含有Σ 2个元素的向量ΔA；接着，用水中超声波的衰减系数加上所述衰减系数变化量ΔA，即可得到待测对象的衰减系数重建值向量。
按照本发明的又一方面，本发明提供了一种利用上述基于射线理论的超声CT声速重建方法的超声CT图像重建方法，其特征在于，该方法利用如权利要求1所述基于射线理论的超声CT声速重建方法，还包括以下步骤：
(4)成像：将得到的待测对象的速度重建值向量二维化，形成Σ×Σ矩阵；然后基于该Σ×Σ矩阵得到二维像素图，该二维像素图中的每个像素与声速值相对应。
作为本发明的进一步优选，所述二维像素图是对声速值进行对数压缩、灰度映射、并最终显示得到。
按照本发明的再一方面，本发明提供了一种利用如权利要求2所述基于射线理论的超声CT衰减系数重建方法的超声CT图像重建方法，其特征在于，该方法利用如权利要求2所述基于射线理论的超声CT衰减系数重建方法，还包括以下步骤：
(4)成像：将得到的待测对象的衰减系数重建值向量二维化，形成Σ×Σ矩阵；然后基于该Σ×Σ矩阵得到二维像素图，该二维像素图中的每个像素与衰减系数值相对应。
作为本发明的进一步优选，所述步骤(4)中，所述二维像素图是对衰减系数值进行对数压缩、灰度映射、并最终显示得到。
按照本发明的再一方面，本发明提供了一种基于射线理论的超声CT声速重建系统，其特征在于，该系统包括：
渡越时间之差提取模块，用于：基于同一发射阵元和接收阵元，发射超声波后分别采集来自纯水和待测对象的超声透射波的数据，分别得到按通道一一对应的纯水数据、以及待测对象数据；采用AIC法提取纯水数据的渡越时间，记为tof water；确定匹配窗，以tof water为窗起点，以纯水数据的最大幅值处时间t water_max为窗终点，则，窗长度记为w，匹配窗内的窗 数据记为W water；在对应通道的待测对象数据上寻找窗长度保持为w的滑动窗，该滑动窗内的窗数据记为W object；将W object和W water互相关，计算得到互相关系数；调整所述滑动窗的窗起点，从而滑动所述滑动窗得到一系列互相关系数，选取这些互相关系数中数值最大的互相关系数对应的滑动窗作为滑动窗寻找结果，记该滑动窗寻找结果的窗起点为tof object，则渡越时间之差Δtof＝tof object-tof water；调整通道，重复提取处理，最终得到一系列通道对应的一系列渡越时间之差Δtof；
声波从发射阵元到接收阵元经过的射线路径计算模块，用于：记得到的一系列有效渡越时间之差Δtof的总数量为n t，对成像区域进行剖分，使成像区域的网格数满足Σ×Σ；其中Σ为正整数，当
Figure PCTCN2019084712-appb-000009
为整数时，则Σ等于
Figure PCTCN2019084712-appb-000010
当
Figure PCTCN2019084712-appb-000011
为非整数时，Σ为
Figure PCTCN2019084712-appb-000012
按向上取整、向下取整或四舍五入取整得到的整数；所述成像区域是建立在发射阵元和接收阵元对应的二维平面直角阵元坐标系中，对于每一组发射阵元-接收阵元组，从发射阵元在阵元坐标系中的坐标位置向接收阵元在阵元坐标系中的坐标位置作直线连线，得到该连线在该成像区域中每一个网格中的路径长度，得到关于每个网格中路径长度的、Σ×Σ的二维矩阵，将该矩阵向量化得到关于每个网格中路径长度的向量；将每一组发射阵元-接收阵元组得到的关于每个网格中路径长度的向量排列形成关于全部发射阵元-接收阵元组的Σ 2×Σ 2的二维矩阵路径矩阵L；
反问题求解模块，用于：将得到的有效渡越时间之差Δtof进行向量化，得到ΔT；构建如式(3)所示的路径-慢度-时间方程组：
LΔS＝ΔT   (3)
其中，ΔS为待求解的慢度变化量；
采用拟牛顿法解方程组，求解得到一维含有Σ 2个元素的向量ΔS；用水中超声波的慢度加上慢度变化量ΔS，取倒数，得到待测对象的速度重建值 向量。
按照本发明的最后一方面，本发明提供了一种基于射线理论的超声CT衰减系数重建系统，其特征在于，该系统包括：
能量参数提取模块，用于：基于同一发射阵元和接收阵元，发射超声波后分别采集来自纯水和待测对象的超声透射波的能量参数，分别得到按通道一一对应的纯水能量参数、以及待测对象能量参数，计算得到待测对象能量参数与纯水能量参数之间的比值；调整通道，重复提取与计算处理，最终得到一系列通道对应的能量参数比值；
声波从发射阵元到接收阵元经过的射线路径计算模块，用于：记得到的一系列能量参数比值的总数量为n t，对成像区域进行剖分，使成像区域的网格数满足Σ×Σ；其中Σ为正整数，当
Figure PCTCN2019084712-appb-000013
为整数时，则Σ等于
Figure PCTCN2019084712-appb-000014
当
Figure PCTCN2019084712-appb-000015
为非整数时，Σ为
Figure PCTCN2019084712-appb-000016
按向上取整、向下取整或四舍五入取整得到的整数；所述成像区域是建立在发射阵元和接收阵元对应的二维平面直角阵元坐标系中，对于每一组发射阵元-接收阵元组，从发射阵元在阵元坐标系中的坐标位置向接收阵元在阵元坐标系中的坐标位置作直线连线，得到该连线在该成像区域中每一个网格中的路径长度，得到关于每个网格中路径长度的、Σ×Σ的二维矩阵，将该矩阵向量化得到关于每个网格中路径长度的向量；将每一组发射阵元-接收阵元组得到的关于每个网格中路径长度的向量排列形成关于全部发射阵元-接收阵元组的Σ 2×Σ 2的二维矩阵路径矩阵L；
反问题求解模块，用于：将得到的一系列能量参数比值进行向量化，得到ΔP；构建如式(4)所示的路径-衰减-能量参数方程组：
LΔA＝ΔP   (4)
其中，ΔA为待求解的衰减系数变化量；
采用拟牛顿法解方程组，求解得到一维含有Σ 2个元素的向量ΔA；用水 中超声波的衰减系数加上所述衰减系数变化量ΔA，得到待测对象的衰减系数重建值向量。
通过本发明所构思的以上技术方案，与现有技术相比，由于基于射线理论，并且在重建过程中采用渡越时间之差(Δtof)或能量参数的比值等参量，可实现快速、稳定的超声CT声速重建及超声CT衰减系数重建。所谓基于射线理论，即假设传播介质较为均匀，超声在均匀介质中的传播路径可近似为射线。
以基于射线理论的声速重建方法为例，其主要步骤包括渡越时间的提取，射线路径的计算，及反问题的求解：
首先对数据进行渡越时间之差(Δtof)的提取。渡越时间即波形到达接收位置的时间，该时间是接收波形中首达波形的起跳点时间。渡越时间之差是指有体模的数据和纯水数据的渡越时间的差值。本发明采用互相关法(CC)与互信息法(AIC)相结合的方法，并引入最大幅度位置信息(MAX)，提取体模数据和纯水数据的渡越时间之差Δtof，简称AIC-MAX-CC法。
接着计算声波从发射阵元到接收阵元经过的射线路径。本发明基于传播介质较为均匀的前提，因此路径用直线路径近似。问题转换为计算两点连线经过网格的交点获取问题。获得交点之后，依次计算出两两交点之间的距离，即为单个网格中路径的长度。
最后是建立并求解路径-慢度-时间方程组，即反问题的建立与求解。慢度是声速的倒数，本发明中的慢度指体模相对于纯水的慢度增量，路径是传播路径在网格中跨越的长度，时间是体模和纯水数据的渡越时间之差。本发明采用拟牛顿法来求解该方程组，拟牛顿法作为一种迭代求解方程组的方法，介于梯度下降法和牛顿法之间，是一种良好的迭代优化方法。
而以基于射线理论的衰减系数重建方法，除了采用能量参数的比值替代渡越时间之差(Δtof)，并用路径-衰减-能量参数方程组替代路径-慢度-时间方程组外，与基于射线理论的声速重建方法大致相似。
总体而言，本发明提出的基于射线理论的声速重建方法及系统的优点有以下几点：1、方程组的构建采用渡越时间之差，而不是渡越时间本身，减小了系统误差的影响；2、渡越时间之差的提取结合了互相关法和互信息法，并引入最大幅度位置信息，简称AIC-MAX-CC法，具有较强的抗噪能力；3、路径计算简单，基于传播介质声速均匀的假设，路径计算转化为计算两点连线经过网格交点的问题，简单高效；4、采用拟牛顿法求解反问题，避免了海森矩阵的计算，运算量小，求解精度高。而基于射线理论的衰减系数重建方法及系统的优点则有以下几点：1、方程组的构建采用能量参数之比，而不是能量本身，减小了系统误差的影响；2、路径计算简单，基于传播介质声速均匀的假设，路径计算转化为计算两点连线经过网格交点的问题，简单高效；3、采用拟牛顿法求解反问题，避免了海森矩阵的计算，运算量小，求解精度高。
【附图说明】
图1是AIC法示意图。
图2是AIC-MAX-CC法示意图，其中，上图对应纯水数据，下图对应体模数据。
图3是直线路径示意图。
图4是超声CT环阵(即UCT环阵)和乳腺定制体模1768-00(CIRSINC,USA)。
图5是对该体模利用超声CT扫描得到的某一层数据重建的反射图像。
图6是本发明提出的声速重建算法对该层数据重建的结果(即，本发明算法重建出的声速图像)。
图7是本发明实施例2中声速重建算法重建得到的衰减系数图像。
【具体实施方式】
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体 实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。
实施例1
本发明中基于射线理论的声速重建方法是将超声CT系统采集到的透射数据进行适当的处理，最后重建得到超声声速图像，其步骤包括渡越时间的提取，射线路径的计算，反问题的求解，声速图像的显示。
采集数据来自于超声CT硬件系统，硬件系统主要包括超声探头、发射电路模块、接收电路模块、数据采集模块，以及计算机系统。它们之间的连接关系、配合工作方式等可直接参考现有技术，如文献(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)。发射电路模块、接收电路模块、数据采集模块等硬件模块，可参考现有技术进行构建。超声探头可以是环形阵列探头、线阵探头、面阵探头等。以环形阵列探头为例，利用环形发射的模式，采集透射数据。透射波主要由发射阵元位置对面的阵元接收，因此本发明采用发射阵元对面的阵元接收到的数据重建声速。若采用其他线阵探头、面阵探头，也可相似处理。
首先是采用AIC法提取纯水数据的渡越时间，记为tof water。
AIC法的算法流程如下：首先在数据上预估的渡越时间点附近选取合适的窗，窗长度为N。以当前检索点为分界点，将窗分为两段，计算当前检索点的AIC值，对窗内的点依次检索，记AIC值最小的点为渡越时间点。
其中，最初预估的渡越时间点，可参考现有技术中的常规操作，例如，可以根据水的理论声速来计算大概的渡越时间点，或者利用肉眼在波形上观察渡越时间点所在的大概位置。窗长度N的选择也是根据波形的实际形状，经验选取，如果不合适可以适当增大或减小。AIC值的计算方法可直接参考相关现有技术。
接着确定匹配窗，以tof water为窗起点，以纯水数据的最大幅值处时间t water_max为窗终点，窗长度记为w，匹配窗窗数据记为W water。
随后在对应通道的体模数据上寻找窗长度为w的滑动窗，滑动窗窗数据记为W object。将W object和W water互相关，计算得到互相关系数，选取互相关系数最大的滑动点数记为p(滑动窗起点在tof water之前则p为负，滑动窗起点在tof water之后则p为正)，则Δtof＝p。如图2所示。也就是说，将滑动窗的窗数据长度保持为w，然后进行滑动，每滑动一个点，计算一个互相关系数，从中选取互相关系数最大时的对应滑动点数记为p。
由于有多个通道，每个通道的渡越时间之差均需要分别求解。
随后要对成像区域进行剖分。网格大小的确定需要根据有效的渡越时间差个数，也即路径-慢度-时间方程组的有效方程数目来确定。无效的渡越时间差，例如有些通道信号缺失，则属于无效。
为了保持方程的正定性，未知数的个数和方程的个数需要保持一致。在提取渡越时间差之后，要剔除掉错误的通道，假设剔除之后有效的渡越时间差个数是n t，则网格数为
Figure PCTCN2019084712-appb-000017
若
Figure PCTCN2019084712-appb-000018
为小数，则取整即可(采用向上取整、向下取整或四舍五入取整中的任意一种取整规则均可)。
接着计算声波从发射阵元到接收阵元经过的射线路径。本发明基于传播介质较为均匀的前提，因此声波的传播路径可以用直线近似代替。问题转换为获取两阵元位置连线经过网格的交点问题。获得交点之后，依次计算出两两交点之间的距离，该距离即为单个网格中路径的长度(若某个网格中只有一个交点，则在该网格中的路径长度为0)。如附图3所示，是直线路径的示意图，S为发射阵元，R为接收阵元。将
Figure PCTCN2019084712-appb-000019
矩阵向量化(即，将
Figure PCTCN2019084712-appb-000020
二维矩阵变成一维列向量)，在向量化之后，空白格是零值，代表该格子不被路径跨过，灰色格子是非零值，该非零值即为该网络中的路径长度，代表该格子被路径跨过。则n t条路径的路径矩阵L大小是 n t×n t。
对于上述网格，可以根据探头出厂给定的尺寸参数，利用现有技术中的处理方法，将各个阵元对应一个坐标系当中，并得到相应的阵元坐标。
在路径计算完成之后，构建路径-慢度-时间方程组，如下式(3)所示，其中，ΔS是慢度变化量，ΔT是渡越时间之差，L是n t×n t的路径矩阵；方程组(3)的维度是有效的渡越时间之差的数目。
LΔS＝ΔT   (3)
采用拟牛顿法解方程组，输出ΔS(得到的ΔS是个向量，取值中可以有正有负)，用水的慢度加上此慢度变化量，取倒数，则为体模的速度重建值。得到的该体模的速度重建值，同样是向量形式。
成像步骤：将向量形式的体模的速度重建值，二维化(即上述
Figure PCTCN2019084712-appb-000021
矩阵向量化的逆过程)后得到二维像素图，该二维像素图中的每个像素即为声速值。
实施例2
与实施例1相似，本发明中基于射线理论的衰减系数重建方法，包括：
本发明中基于射线理论的衰减系数重建方法是将超声CT系统采集到的透射数据进行适当的处理，最后重建得到超声衰减图像，其步骤包括信号能量参数的提取，射线路径的计算，反问题的求解，衰减系数图像的显示。
能量参数可以是信号的幅值、强度、能量。
射线路径的计算同实施例1。
反问题的求解同实施例1。
在路径计算完成之后，构建路径-衰减-能量参数方程组，如下式(4)所示，其中，ΔA是衰减系数变化量，ΔP是能量参数的比值，即仿体数据的能量参数与水数据的能量参数的比值。
LΔA＝ΔP  (4)
采用拟牛顿法求解(4)，得到ΔA，用水的衰减系数加上此衰减系数变化量， 则为体模的衰减系数重建值。
最后是衰减系数的图像显示，如图7所示。
上述实施例中，对于得到的超声图像，还可以按照现有技术中的方法依次进行对数压缩和灰度映射，得到灰度值在规定范围内的超声图像。
本发明适用于各种商业的超声CT系统探头，如环形探头、线阵探头、面阵探头、凹阵等。
本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (8)

  1. 一种基于射线理论的超声CT声速重建方法，其特征在于，包括以下步骤：
    (1)渡越时间之差的提取：
    首先基于同一发射阵元和接收阵元，发射超声波后分别采集来自纯水和待测对象的超声透射波的数据，分别得到按通道一一对应的纯水数据、以及待测对象数据；
    然后采用AIC法提取纯水数据的渡越时间，记为tof water；接着，确定匹配窗，以tof water为窗起点，以纯水数据的最大幅值处时间t water_max为窗终点，则，窗长度记为w，匹配窗内的窗数据记为W water；
    随后在对应通道的待测对象数据上寻找窗长度保持为w的滑动窗，该滑动窗内的窗数据记为W object；接着，将W object和W water互相关，计算得到互相关系数；调整所述滑动窗的窗起点，从而滑动所述滑动窗得到一系列互相关系数，选取这些互相关系数中数值最大的互相关系数对应的滑动窗作为滑动窗寻找结果，记该滑动窗寻找结果的窗起点为tof object，则渡越时间之差Δtof＝tof object-tof water；
    接着，调整通道，重复提取处理，最终得到一系列通道对应的一系列渡越时间之差Δtof；
    (2)计算声波从发射阵元到接收阵元经过的射线路径：
    记所述步骤(1)得到的一系列有效渡越时间之差Δtof的总数量为n t，对成像区域进行剖分，使成像区域的网格数满足Σ×Σ；其中Σ为正整数，当
    Figure PCTCN2019084712-appb-100001
    为整数时，则Σ等于
    Figure PCTCN2019084712-appb-100002
    当
    Figure PCTCN2019084712-appb-100003
    为非整数时，Σ为
    Figure PCTCN2019084712-appb-100004
    按向上取整、向下取整或四舍五入取整得到的整数；
    所述成像区域是建立在发射阵元和接收阵元对应的二维平面直角阵元坐标系中，然后对于每一组发射阵元-接收阵元组，从发射阵元在阵元坐标 系中的坐标位置向接收阵元在阵元坐标系中的坐标位置作直线连线，得到该连线在该成像区域中每一个网格中的路径长度，进而得到关于每个网格中路径长度的、Σ×Σ的二维矩阵，接着将该矩阵向量化得到关于每个网格中路径长度的向量；最后，将每一组发射阵元-接收阵元组得到的关于每个网格中路径长度的向量排列形成关于全部发射阵元-接收阵元组的Σ 2×Σ 2的二维矩阵路径矩阵L；
    (3)反问题的求解：
    将所述步骤(1)得到的有效渡越时间之差Δtof进行向量化，得到ΔT；接着构建如式(3)所示的路径-慢度-时间方程组：
    LΔS＝ΔT  (3)
    其中，ΔS为待求解的慢度变化量；然后，采用拟牛顿法解方程组，求解得到一维含有Σ 2个元素的向量ΔS；接着，用水中超声波的慢度加上慢度变化量ΔS，取倒数，即可得到待测对象的速度重建值向量。
  2. 一种基于射线理论的超声CT衰减系数重建方法，其特征在于，包括以下步骤：
    (1)首先基于同一发射阵元和接收阵元，发射超声波后分别采集来自纯水和待测对象的超声透射波的能量参数，分别得到按通道一一对应的纯水能量参数、以及待测对象能量参数，然后计算得到待测对象能量参数与纯水能量参数之间的比值；
    接着，调整通道，重复提取与计算处理，最终得到一系列通道对应的能量参数比值；
    (2)计算声波从发射阵元到接收阵元经过的射线路径：
    记所述步骤(1)得到的一系列能量参数比值的总数量为n t，对成像区域进行剖分，使成像区域的网格数满足Σ×Σ；其中Σ为正整数，当
    Figure PCTCN2019084712-appb-100005
    为整数时，则Σ等于
    Figure PCTCN2019084712-appb-100006
    当
    Figure PCTCN2019084712-appb-100007
    为非整数时，Σ为
    Figure PCTCN2019084712-appb-100008
    按向上取整、向下取 整或四舍五入取整得到的整数；
    所述成像区域是建立在发射阵元和接收阵元对应的二维平面直角阵元坐标系中，然后对于每一组发射阵元-接收阵元组，从发射阵元在阵元坐标系中的坐标位置向接收阵元在阵元坐标系中的坐标位置作直线连线，得到该连线在该成像区域中每一个网格中的路径长度，进而得到关于每个网格中路径长度的、Σ×Σ的二维矩阵，接着将该矩阵向量化得到关于每个网格中路径长度的向量；最后，将每一组发射阵元-接收阵元组得到的关于每个网格中路径长度的向量排列形成关于全部发射阵元-接收阵元组的Σ 2×Σ 2的二维矩阵路径矩阵L；
    (3)反问题的求解：
    将所述步骤(1)得到的一系列能量参数比值进行向量化，得到ΔP；接着构建如式(4)所示的路径-衰减-能量参数方程组：
    LΔA＝ΔP  (4)
    其中，ΔA为待求解的衰减系数变化量；然后，采用拟牛顿法解方程组，求解得到一维含有Σ 2个元素的向量ΔA；接着，用水中超声波的衰减系数加上所述衰减系数变化量ΔA，即可得到待测对象的衰减系数重建值向量。
  3. 利用如权利要求1所述基于射线理论的超声CT声速重建方法的超声CT图像重建方法，其特征在于，该方法利用如权利要求1所述基于射线理论的超声CT声速重建方法，还包括以下步骤：
    (4)成像：将得到的待测对象的速度重建值向量二维化，形成Σ×Σ矩阵；然后基于该Σ×Σ矩阵得到二维像素图，该二维像素图中的每个像素与声速值相对应。
  4. 如权利要求3所述超声CT图像重建方法，其特征在于，所述步骤(4)中，所述二维像素图是对声速值进行对数压缩、灰度映射、并最终显示得到。
  5. 利用如权利要求2所述基于射线理论的超声CT衰减系数重建方法 的超声CT图像重建方法，其特征在于，该方法利用如权利要求2所述基于射线理论的超声CT衰减系数重建方法，还包括以下步骤：
    (4)成像：将得到的待测对象的衰减系数重建值向量二维化，形成Σ×Σ矩阵；然后基于该Σ×Σ矩阵得到二维像素图，该二维像素图中的每个像素与衰减系数值相对应。
  6. 如权利要求5所述超声CT图像重建方法，其特征在于，所述步骤(4)中，所述二维像素图是对衰减系数值进行对数压缩、灰度映射、并最终显示得到。
  7. 一种基于射线理论的超声CT声速重建系统，其特征在于，该系统优选包括：
    渡越时间之差提取模块，用于：基于同一发射阵元和接收阵元，发射超声波后分别采集来自纯水和待测对象的超声透射波的数据，分别得到按通道一一对应的纯水数据、以及待测对象数据；采用AIC法提取纯水数据的渡越时间，记为tof water；确定匹配窗，以tof water为窗起点，以纯水数据的最大幅值处时间t water_max为窗终点，则，窗长度记为w，匹配窗内的窗数据记为W water；在对应通道的待测对象数据上寻找窗长度保持为w的滑动窗，该滑动窗内的窗数据记为W object；将W object和W water互相关，计算得到互相关系数；调整所述滑动窗的窗起点，从而滑动所述滑动窗得到一系列互相关系数，选取这些互相关系数中数值最大的互相关系数对应的滑动窗作为滑动窗寻找结果，记该滑动窗寻找结果的窗起点为tof object，则渡越时间之差Δtof＝tof object-tof water；调整通道，重复提取处理，最终得到一系列通道对应的一系列渡越时间之差Δtof；
    声波从发射阵元到接收阵元经过的射线路径计算模块，用于：记得到的一系列有效渡越时间之差Δtof的总数量为n t，对成像区域进行剖分，使成像区域的网格数满足Σ×Σ；其中Σ为正整数，当
    Figure PCTCN2019084712-appb-100009
    为整数时，则Σ等 于
    Figure PCTCN2019084712-appb-100010
    当
    Figure PCTCN2019084712-appb-100011
    为非整数时，Σ为
    Figure PCTCN2019084712-appb-100012
    按向上取整、向下取整或四舍五入取整得到的整数；所述成像区域是建立在发射阵元和接收阵元对应的二维平面直角阵元坐标系中，对于每一组发射阵元-接收阵元组，从发射阵元在阵元坐标系中的坐标位置向接收阵元在阵元坐标系中的坐标位置作直线连线，得到该连线在该成像区域中每一个网格中的路径长度，得到关于每个网格中路径长度的、Σ×Σ的二维矩阵，将该矩阵向量化得到关于每个网格中路径长度的向量；将每一组发射阵元-接收阵元组得到的关于每个网格中路径长度的向量排列形成关于全部发射阵元-接收阵元组的Σ 2×Σ 2的二维矩阵路径矩阵L；
    反问题求解模块，用于：将得到的有效渡越时间之差Δtof进行向量化，得到ΔT；构建如式(3)所示的路径-慢度-时间方程组：
    LΔS＝ΔT  (3)
    其中，ΔS为待求解的慢度变化量；
    采用拟牛顿法解方程组，求解得到一维含有Σ 2个元素的向量ΔS；用水中超声波的慢度加上慢度变化量ΔS，取倒数，得到待测对象的速度重建值向量。
  8. 一种基于射线理论的超声CT衰减系数重建系统，其特征在于，该系统优选包括：
    能量参数提取模块，用于：基于同一发射阵元和接收阵元，发射超声波后分别采集来自纯水和待测对象的超声透射波的能量参数，分别得到按通道一一对应的纯水能量参数、以及待测对象能量参数，计算得到待测对象能量参数与纯水能量参数之间的比值；调整通道，重复提取与计算处理，最终得到一系列通道对应的能量参数比值；
    声波从发射阵元到接收阵元经过的射线路径计算模块，用于：记得到的一系列能量参数比值的总数量为n t，对成像区域进行剖分，使成像区域 的网格数满足Σ×Σ；其中Σ为正整数，当
    Figure PCTCN2019084712-appb-100013
    为整数时，则Σ等于
    Figure PCTCN2019084712-appb-100014
    当
    Figure PCTCN2019084712-appb-100015
    为非整数时，Σ为
    Figure PCTCN2019084712-appb-100016
    按向上取整、向下取整或四舍五入取整得到的整数；所述成像区域是建立在发射阵元和接收阵元对应的二维平面直角阵元坐标系中，对于每一组发射阵元-接收阵元组，从发射阵元在阵元坐标系中的坐标位置向接收阵元在阵元坐标系中的坐标位置作直线连线，得到该连线在该成像区域中每一个网格中的路径长度，得到关于每个网格中路径长度的、Σ×Σ的二维矩阵，将该矩阵向量化得到关于每个网格中路径长度的向量；将每一组发射阵元-接收阵元组得到的关于每个网格中路径长度的向量排列形成关于全部发射阵元-接收阵元组的Σ 2×Σ 2的二维矩阵路径矩阵L；
    反问题求解模块，用于：将得到的一系列能量参数比值进行向量化，得到ΔP；构建如式(4)所示的路径-衰减-能量参数方程组：
    LΔA＝ΔP  (4)
    其中，ΔA为待求解的衰减系数变化量；
    采用拟牛顿法解方程组，求解得到一维含有Σ 2个元素的向量ΔA；用水中超声波的衰减系数加上所述衰减系数变化量ΔA，得到待测对象的衰减系数重建值向量。
PCT/CN2019/084712 2019-04-11 2019-04-28 一种基于射线理论的超声ct图像重建方法及系统 WO2020206755A1 (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
US16/968,151 US20210236095A1 (en) 2019-04-11 2019-04-28 Ultrasound ct image reconstruction method and system based on ray theory

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN201910286906.7 2019-04-11
CN201910286906.7A CN110051387B (zh) 2019-04-11 2019-04-11 一种基于射线理论的超声ct图像重建方法及系统

Publications (1)

Publication Number Publication Date
WO2020206755A1 true WO2020206755A1 (zh) 2020-10-15

Family

ID=67318659

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2019/084712 WO2020206755A1 (zh) 2019-04-11 2019-04-28 一种基于射线理论的超声ct图像重建方法及系统

Country Status (3)

Country Link
US (1) US20210236095A1 (zh)
CN (1) CN110051387B (zh)
WO (1) WO2020206755A1 (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114577322A (zh) * 2022-02-15 2022-06-03 南京大学 一种基于六边形网格路径计算的声速成像方法

Families Citing this family (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110575202B (zh) * 2019-08-27 2020-06-02 华中科技大学 一种基于费马原理的超声ct图像重建方法及系统
CN110772281B (zh) * 2019-10-23 2022-03-22 哈尔滨工业大学（深圳） 基于改进射线追踪法的超声ct成像系统
CN112764040B (zh) * 2019-11-01 2022-06-14 复旦大学 一种基于射线理论相位修正的合成孔径波束形成方法
CN112363038B (zh) * 2020-09-30 2021-04-27 温州大学 一种宽禁带半导体异质结渡越时间二极管噪声检测方法及系统
CN112674794B (zh) * 2020-12-21 2023-02-10 苏州二向箔科技有限公司 一种结合深度学习与吉洪诺夫正则化反演的超声ct声速重建方法
CN115267673B (zh) * 2022-03-07 2024-05-07 清华大学 考虑重建网格偏移的稀疏声源成像方法、系统
CN114886469B (zh) * 2022-05-11 2024-06-14 中国科学院声学研究所 一种超声ct阵列探头的阵元定位方法和装置
CN117949542B (zh) * 2024-03-26 2024-06-04 中建安装集团有限公司 一种工业管廊抱柱脚手架牢固度检测方法、介质及系统

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5603326A (en) * 1993-03-22 1997-02-18 Siemens Aktiengesellschaft Method and apparatus for displaying an image obtained by echo signals
US20050080334A1 (en) * 2003-10-08 2005-04-14 Scimed Life Systems, Inc. Method and system for determining the location of a medical probe using a reference transducer array
US20160109589A1 (en) * 2014-10-20 2016-04-21 Jonathan Liu Velocity tomography using property scans
CN105997153A (zh) * 2016-07-15 2016-10-12 华中科技大学 一种超声ct的成像方法
CN108348221A (zh) * 2015-09-01 2018-07-31 戴尔菲纳斯医疗科技公司 使用超声波波形断层成像的组织成像和分析

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106153647B (zh) * 2015-04-08 2021-04-13 清华大学 能谱ct成像系统及数据采集和重建能谱ct图像的方法
WO2018022565A1 (en) * 2016-07-25 2018-02-01 The General Hospital Corporation System and method for tomographic image reconstruction

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5603326A (en) * 1993-03-22 1997-02-18 Siemens Aktiengesellschaft Method and apparatus for displaying an image obtained by echo signals
US20050080334A1 (en) * 2003-10-08 2005-04-14 Scimed Life Systems, Inc. Method and system for determining the location of a medical probe using a reference transducer array
US20160109589A1 (en) * 2014-10-20 2016-04-21 Jonathan Liu Velocity tomography using property scans
CN108348221A (zh) * 2015-09-01 2018-07-31 戴尔菲纳斯医疗科技公司 使用超声波波形断层成像的组织成像和分析
CN105997153A (zh) * 2016-07-15 2016-10-12 华中科技大学 一种超声ct的成像方法

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114577322A (zh) * 2022-02-15 2022-06-03 南京大学 一种基于六边形网格路径计算的声速成像方法
CN114577322B (zh) * 2022-02-15 2022-11-22 南京大学 一种基于六边形网格路径计算的声速成像方法

Also Published As

Publication number Publication date
CN110051387B (zh) 2020-05-19
US20210236095A1 (en) 2021-08-05
CN110051387A (zh) 2019-07-26

Similar Documents

Publication Publication Date Title
WO2020206755A1 (zh) 一种基于射线理论的超声ct图像重建方法及系统
WO2021004076A1 (zh) 基于人工智能芯片的适形穿戴式生物信息监测设备及系统
CN102920478B (zh) 一种合成聚焦的便携式b型超声成像方法
CN110164550B (zh) 一种基于多视角协同关系的先天性心脏病辅助诊断方法
US11759177B2 (en) Three-dimensional ultrasound tomography method and system based on spiral scanning
CN109875606B (zh) 基于先验反射成像的超声ct声速成像的方法
CN104688271A (zh) 合成聚焦超声成像方法和装置
CN113850883A (zh) 基于注意力机制的磁粒子成像重建方法
Awasthi et al. Sinogram super-resolution and denoising convolutional neural network (SRCN) for limited data photoacoustic tomography
CN105997153B (zh) 一种超声ct的成像方法
CN110070612A (zh) 一种基于生成对抗网络的ct图像层间插值方法
CN114119362A (zh) 利用神经网络提高超声图像分辨率的系统和方法
CN110974293B (zh) 一种基于c型探头的合成孔径成像方法
CN111248858A (zh) 一种基于频域波数域的光声断层成像重建方法
CN112465924B (zh) 一种基于多特征融合的快速医学图像重构方法
Noda et al. Ultrasound imaging with a flexible probe based on element array geometry estimation using deep neural network
US11191519B2 (en) Device, system, and method for hemispheric breast imaging
CN102688071A (zh) 超声浅表组织与器官容积扫描断层成像方法
CN109766646B (zh) 一种基于稀疏通道回波数据重建的超声成像方法及装置
CN111223157A (zh) 一种基于深度残差网络的超声ct声速成像方法
CN110575202B (zh) 一种基于费马原理的超声ct图像重建方法及系统
Yuan et al. A modified ray tracing method for ultrasound computed tomography in breast imaging
CN114359327A (zh) 一种基于ai的超声心动锚定点动态测量模型及应用
CN115500817A (zh) 一种融合毫米波雷达和深度学习模型的无接触式心震信号检测方法及其装置
JP5959880B2 (ja) 超音波診断装置

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 19924116

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 19924116

Country of ref document: EP

Kind code of ref document: A1