WO2012171388A1 - 一种弹性成像中的位移检测方法及装置 - Google Patents

一种弹性成像中的位移检测方法及装置 Download PDF

Info

Publication number
WO2012171388A1
WO2012171388A1 PCT/CN2012/073135 CN2012073135W WO2012171388A1 WO 2012171388 A1 WO2012171388 A1 WO 2012171388A1 CN 2012073135 W CN2012073135 W CN 2012073135W WO 2012171388 A1 WO2012171388 A1 WO 2012171388A1
Authority
WO
WIPO (PCT)
Prior art keywords
cross
point
data
correlation
correlation phase
Prior art date
Legal status (The legal status 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 status listed.)
Ceased
Application number
PCT/CN2012/073135
Other languages
English (en)
French (fr)
Inventor
李双双
樊睿
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Shenzhen Mindray Bio Medical Electronics Co Ltd
Original Assignee
Shenzhen Mindray Bio Medical Electronics Co Ltd
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 Shenzhen Mindray Bio Medical Electronics Co Ltd filed Critical Shenzhen Mindray Bio Medical Electronics Co Ltd
Priority to US14/126,222 priority Critical patent/US9607405B2/en
Publication of WO2012171388A1 publication Critical patent/WO2012171388A1/zh
Anticipated expiration legal-status Critical
Ceased legal-status Critical Current

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/60—Analysis of geometric attributes
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N29/00—Investigating or analysing materials by the use of ultrasonic, sonic or infrasonic waves; Visualisation of the interior of objects by transmitting ultrasonic or sonic waves through the object
    • G01N29/04—Analysing solids
    • G01N29/06—Visualisation of the interior, e.g. acoustic microscopy
    • G01N29/0654—Imaging
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/48—Diagnostic techniques
    • A61B8/485—Diagnostic techniques involving measuring strain or elastic properties
    • 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
    • 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
    • G01S7/52042—Details of receivers using analysis of echo signal for target characterisation determining elastic properties of the propagation medium or of the reflective target
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/0002—Inspection of images, e.g. flaw detection
    • G06T7/0012—Biomedical image inspection
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2291/00—Indexing codes associated with group G01N29/00
    • G01N2291/02—Indexing codes associated with the analysed material
    • G01N2291/024—Mixtures
    • G01N2291/02475—Tissue characterisation
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2291/00—Indexing codes associated with group G01N29/00
    • G01N2291/02—Indexing codes associated with the analysed material
    • G01N2291/028—Material parameters
    • G01N2291/02827—Elastic parameters, strength or force

Definitions

  • the present invention relates to an elastic imaging technique, and more particularly to a displacement detecting method and apparatus in elastic imaging. Background technique
  • Ultrasound elastography is an important aid to the detection of B-mode sonograms in cancer detection, especially benign malignant differentiation of breast cancer, and is rapidly applied to clinical practice.
  • Ultrasound elastography mainly obtains ultrasonic echo information of the target tissue by means of ultrasonic imaging, and then detects the tissue elasticity information through a specific algorithm, and visually displays it in the form of an image to assist the doctor in diagnosis or treatment.
  • the traditional ultrasound elastography method requires the probe to slightly compress the tissue or use the body's own breathing, vascular pulse and other processes to obtain two consecutive frames of ultrasonic echo signals, and then obtain a displacement between the two frames by a specific displacement detection method.
  • the spatial position change information of the target tissue at two different times and the axial strain of the tissue can be obtained by obtaining an axial gradient of the displacement.
  • This strain information can reflect the elasticity of the tissue. Under the same external force compression, the greater the strain, the softer the tissue, and the smaller the strain, the harder the tissue.
  • the strain information of the target tissue area is expressed in the form of an image, which can directly reflect the soft and hard difference or the elastic difference between different tissues, that is, a strain image. This method is also called strain imaging.
  • the tissue Doppler-based method achieves the most simple, but depends on the Doppler signal-to-noise ratio and processing method, the same There is an aliasing problem, and its lateral displacement is ignored, and the image quality is relatively poor.
  • the main technical problem to be solved by the present invention is to provide a displacement detecting method and apparatus in elastic imaging, which can reduce the calculation amount of the imaging process.
  • the present invention provides a displacement detecting method in elastic imaging, which includes:
  • a gradient is obtained from the displacement result to obtain a strain result.
  • the invention also provides a displacement detecting device in elastic imaging, comprising:
  • a target point obtaining device configured to acquire a target point
  • a search device configured to acquire a cross-correlation phase calculation position of the target point in the second frame image
  • a cross-correlation phase calculation device configured to calculate a cross-correlation phase according to the cross-correlation phase calculation position
  • the longitudinal displacement result calculating means is configured to calculate a longitudinal displacement result according to the cross-correlation phase; and the strain result calculating means is configured to obtain a strain result by obtaining a gradient of the displacement result.
  • the beneficial effects of the present invention are as follows:
  • the elastic imaging method and apparatus according to the embodiments of the present invention acquires the I/Q two-channel echo baseband signals after two frames of down-conversion before and after compression, and uses guided phase estimation to quickly detect the displacement between two frames.
  • Information, and then the axial gradient calculation to obtain the strain information not only can obtain higher quality strain images, but also greatly reduce the amount of calculation, to meet the clinical real-time requirements.
  • FIG. 1 is a perspective view of an embodiment of an elastic imaging system of the present invention
  • FIG. 2 is a flowchart of an embodiment of a method for detecting displacement in elastic imaging according to the present invention
  • FIG. 3 is a schematic diagram of frame data gridization in an embodiment of the present invention
  • FIG. 4 is a schematic diagram of a displacement search strategy in an embodiment of the present invention.
  • Fig. 5 is a block diagram showing an embodiment of a displacement detecting device in elastic imaging according to the present invention. Detailed ways
  • the ultrasonic probe transmits ultrasonic waves and receives echo information according to the preset scanning rules of the system.
  • the radio frequency (RF) signal is output, and then the quadrature demodulation phase is generated to generate I.
  • the downsampling rate is preset by the system, and then the displacement phase detection is performed by the guided phase zero estimation (GPZE)
  • GPZE guided phase zero estimation
  • the GPZE algorithm mainly detects the phase of the cross-correlation function between the two frames before and after, and derives the correspondence between the phase and the longitudinal displacement by guiding the search thinking, and calculates the longitudinal displacement between the two frames.
  • the quantity while greatly reducing the amount of calculation, ensures the quality of the displacement estimation.
  • the range of displacement detection has been extended.
  • the GPZE algorithm is not only suitable for small displacement but also for large displacement.
  • the cross-correlation function (., x .) can be calculated by using the baseband data of the corresponding offset positions of the two frames before and after, due to the cross-correlation function and R 'The function only differs in phase, so that the above cross-correlation phase can be further calculated.
  • the time shift or the longitudinal shift between the two frames is finally expressed as:
  • a displacement detecting method in the elastic imaging provided by the present invention includes:
  • the longitudinal displacement result is calculated from the cross-correlation phase.
  • the baseband signal data before acquiring the target point, the baseband signal data may be acquired, where the baseband signal data is meshed, the nodes of the mesh are displacement detection estimation points, and the target points are from the displacement detection estimation points.
  • the cross-correlation phase calculation position of the acquisition target point in the second frame picture may include the data longitudinal position of the cross-correlation phase calculation.
  • the step of acquiring the cross-correlation phase calculation position of the target point in the second frame image comprises: obtaining a longitudinal displacement result of the previous calculation point of the target point; in the second frame image, calculating the position of the target point and the previous one
  • the sum of the longitudinal displacement results of the points is the point at which the central region searches for the greatest correlation with the nuclear data, and the longitudinal position of the point is the longitudinal position of the cross-correlation phase calculation of the target point in the second frame image.
  • the cross-correlation phase calculation position of the acquisition target point in the second frame picture may also include the data lateral position of the cross-correlation phase calculation.
  • the step of acquiring the cross-correlation phase calculation position of the target point in the second frame image comprises: laterally searching for the point with the highest data correlation in the vicinity of the longitudinal position of the cross-correlation phase calculation, and the lateral position of the point is the data of the cross-correlation phase calculation. Lateral position.
  • the step of laterally searching for the point at which the data correlation is most near the longitudinal position of the cross-correlation phase calculation is calculated once every interval of a certain number of points.
  • the position obtained in the step of acquiring the cross-correlation phase calculation position of the target point in the second frame image includes the data longitudinal position calculated by the cross-correlation phase and the data lateral position calculated by the cross-correlation phase;
  • the steps of calculating the position of the cross-correlation phase in the second frame image include:
  • the lateral position of the point is the lateral position of the cross-correlation calculation;
  • the region centered on the sum of the position of the target point and the longitudinal displacement result of the previous calculation point of the target point searches for the point with the greatest correlation with the kernel data, and the longitudinal position of the point is the target point in the second frame.
  • the longitudinal position of the cross-correlation phase calculation in the image is the longitudinal position of the cross-correlation phase calculation in the image.
  • FIG. 2 is a flow chart of an embodiment of a displacement detecting method in elastography according to the present invention, including:
  • the calculation of the displacement data for each frame requires the use of two frames of I/Q baseband signal data, and the resulting displacement data refers to the spatial relative displacement between the two frames of signals.
  • the two frames of baseband signal data may be two consecutive frames of data, and may be two frames of data with a certain frame interval, and the number of frame intervals is preset by the system. Using two frames of data at a certain interval, the amount of displacement between two frames of data used for calculation can be effectively adjusted. Small, making the resulting strain image better.
  • the baseband signal data of each frame is divided into two data of I and Q.
  • the sampling rate of the baseband signal may be different from the sampling rate of the RF data, and the downsampling rate is preset by the system. The larger the downsampling rate, the smaller the amount of computation will be, but it will affect the quality of the displacement estimate and the spatial resolution of the final image.
  • the envelope data needs to be further calculated for the search calculation of the subsequent envelope offset.
  • the calculation of the envelope is performed for each sample point of the I/Q data, and the envelope data at the sample point ' is calculated as:
  • the displacement estimation does not necessarily calculate the position of each sampling point of the baseband signal data, and the position of the displacement detection estimation point is pre-divided, which can effectively avoid or reduce the amount of redundant calculation.
  • FIG. 3 is a schematic diagram of frame data meshing in an embodiment of the present invention.
  • a black dot ie, a mesh node
  • a black line in the figure indicates baseband signal data or envelope data.
  • the meshing is based on the position of the data sample point of the first frame of the two frame signals.
  • the processes of offset search, phase shift detection, displacement calculation, strain estimation, and the like are preferably performed with each node in the above-described mesh (i.e., displacement detection estimation point) as a target.
  • the cross-correlation phase calculation position it is necessary to find the longitudinal position of the data of the cross-correlation phase calculation of the target point in the second frame image, that is, to determine the data longitudinal direction of the cross-correlation phase calculation.
  • the offset that is, the vertical offset when the envelope cross-correlation function between the two frames of data is maximized.
  • the search for the longitudinal offset is based on the idea of block-matching and searches for the envelope data.
  • the core data is selected as the center, and the position with the greatest correlation with the core data is searched longitudinally in the second frame data.
  • the offset between the position and the target point is the desired longitudinal offset.
  • the most relevant discrimination in the search can be based on the SAD method, the NCC method, and the like. The smallest position of the SAD or the largest position of the NCC is the position with the most correlation. Other similar criteria can also be used.
  • the invention also incorporates a guiding thinking on the basis of block-matching to improve the calculation speed.
  • ROI region of interest
  • the guided search can have two different implementation methods, the difference is that the depth of the data calculation starts differently, and the corresponding calculation amount is also different.
  • Method 1 Regardless of the depth of the ROI (region of interest) selected by the user, the calculation always starts from the data of the probe surface (ie, the node with a depth of 0 or closest to 0), and the longitudinal displacement of the node at each depth above
  • a small vertical search range is set to search based on this initial value, thereby obtaining the maximum longitudinal offset of the envelope cross-correlation function.
  • the system sets its corresponding "previous depth node” longitudinal displacement result to be fixed at zero.
  • the initial position of the envelope longitudinal offset corresponding to the current node is ", that is, the envelope of the second frame is required.
  • the data is set with a small vertical search range centered on the depth + , and the position where the correlation between the two is the largest is found, which is the longitudinal offset of the envelope corresponding to the current node.
  • the vertical search range is preset by the system. The smaller the range, the smaller the amount of calculation.
  • Method 2 The calculation always starts from the shallowest depth node in the user-selected R0I, and the displacement result of each time the depth of the node is the initial value of the vertical data position offset of the node of the current depth, in this Based on the initial value, a small vertical search range is set for searching, thereby obtaining the maximum longitudinal offset of the envelope cross-correlation function. For the shallowest depth node in the ROI, special processing is required because the final longitudinal displacement result of the "node of the previous depth" is unknown. It is assumed that the depth position is a large vertical search range in the second frame data in the search, and the position with the highest correlation is found as the envelope longitudinal offset corresponding to the node.
  • the vertical search range is preset by the system and generally needs to be set larger than when there is a guided search (that is, when there is a search initial value). If the setting is too small, the search result may be inaccurate.
  • the displacement result of the node of the previous depth is u y
  • the depth position of the current node is the initial offset of the envelope corresponding to the current node.
  • the initial position of the search is + that is required in the second frame envelope.
  • the data is set with a small vertical search range centered on the depth + , and the position where the correlation between the two is the largest is found, which is the longitudinal offset of the envelope corresponding to the current node.
  • the vertical search range is preset by the system, and the smaller the range, the smaller the calculation amount. Assuming that the result of the longitudinal offset of the envelope obtained by the current node is "., the longitudinal position of the data in the first frame and the second frame corresponding to the current node in the cross-correlation phase calculation is respectively ⁇ and ". .
  • the longitudinal position After determining the longitudinal position, preferably, it is also possible to continue to find the lateral position, i.e., to find the lateral offset when the two-frame data envelope cross-correlation function is maximized.
  • the search of the lateral position is guided by the result of the above longitudinal position, that is, after the longitudinal position is obtained, the longitudinal position is found in the second frame, and a certain horizontal search range (or search area) is set centering on the longitudinal position.
  • the lateral offset of the position where the cross-correlation of the two-frame envelope data is the largest relative to the original position within the range is obtained.
  • the lateral position of the second frame data at the final cross-correlation phase calculation is selected based on the lateral offset.
  • This horizontal search range is preset by the system. The smaller the search range, the smaller the amount of calculation.
  • the correlation discrimination in the search is based on a search similar to the vertical position, and SAD, NCC or other discriminating methods can be used.
  • the two-dimensional position of the current target node in the first frame data is , where is the amount of time shift, and ⁇ is the position of the scan line.
  • Set the longitudinal offset of the envelope obtained in the previous step to "., when searching, set a small horizontal search area centered on the position, and find the position where the two-frame envelope data has the largest correlation in the region.
  • the offset from the original position is the obtained lateral offset.
  • the data position of the first frame and the second frame corresponding to the current node in the cross-correlation phase calculation is ⁇ and ⁇ + ". , ).
  • the search operation for the above lateral position does not need to be continuously performed along the nodes of each depth, and the horizontal position can be updated by searching for the nodes at intervals of several depths. This can further reduce the amount of calculation while ensuring that the image quality does not change significantly.
  • the size of the depth interval of the horizontal position update is preset by the system. If the interval is too small, the update is frequent, which may increase the amount of redundancy calculation. If the interval is too large, the update is less, which may affect the cross-correlation phase estimation quality.
  • cross-correlation phase estimation can be performed.
  • the phase correlation calculation is performed using two frames of I/Q data before and after.
  • the size of the block data is preset by the system, and the size of the block has an influence on the cross-correlation phase estimation result.
  • the data blocks involved in the calculation are called kernels and can be one-dimensional or two-dimensional.
  • the position coordinates of the current target point (ie, the current node) in the first frame data are (n), and the longitudinal and lateral offsets of the envelopes searched in the first two steps are respectively ". and .
  • the position coordinate corresponding to the second frame data is + ". , ⁇ + ), in Fig. 4, with ( ⁇ ) and + “. ⁇ ) as the center or the reference point extending a certain distance in the longitudinal and lateral directions, the area in the figure is obtained, and the data in this area is the nuclear data.
  • the phase calculation uses I and Q data.
  • the I and Q core data of the corresponding positions in the first frame and the second frame are taken out, and the phase is obtained as follows:
  • the relative positions indicated in the data are the same, for example, assuming (0, 0) indicates the point in the lower left corner of the core data in the core data of the first frame, Then (0, 0) also indicates the point in the lower left corner of the core data in the core data of the second frame.
  • the calculation method is related to the final displacement estimation result of the previous depth node, assuming that the final longitudinal displacement result of the previous depth node is , B' J :
  • the longitudinal displacement result between the final two frames is: Wherein, it is the signal period, which corresponds to the angular frequency of the signal center.
  • the introduction of " ⁇ / 2 term compensates for the aliasing of the phase calculation, so that the algorithm is suitable for small displacements and for large displacements.
  • the above longitudinal displacements are expressed in units of sampling time. , can also be converted into physical length units to represent, the two are one-to-one correspondence.
  • the displacement data is graded along the longitudinal direction to obtain the strain result, that is, the strain value.
  • the resulting strain results can be corrected for errors, such as detecting abnormal jump points or obvious error points for correction; or spatially smoothing them to improve image display; or Map with different grayscale or color maps, To enhance image contrast. Other operations that increase image quality can also be performed.
  • strain results of the desired area are output and displayed as strain images, which can reflect the difference in tissue elasticity of the region.
  • the GPZE displacement detection algorithm uses the displacement result of the previous depth to guide the displacement calculation of the next depth, which reduces the amount of search calculation; on the other hand, the phase estimation method is used to calculate the displacement, and the sampling rate of the original data is required. Not high, which greatly reduces the amount of calculation.
  • lateral search can also be introduced to improve the estimation quality of the longitudinal displacement, so that the final image signal-to-noise ratio is higher.
  • phase aliasing is compensated, which extends the applicable range of displacement estimation.
  • the longitudinal position of the data of the cross-correlation phase calculation is searched for, and then the lateral position is searched.
  • the lateral position may be searched first, and then the longitudinal position may be searched; or, the above-mentioned search for the lateral position and the longitudinal position is performed.
  • the search method of the lateral position may be centered on the target point, and set a certain horizontal search range, in which the lateral offset of the position where the cross-correlation of the two-frame envelope data is the largest relative to the original position is begging.
  • the method for finding the longitudinal position may be, for example, searching for a point having the greatest correlation with the nuclear data, centering on the sum of the position of the target point and the longitudinal displacement result of the earlier calculated point, and the longitudinal position of the point is the target The longitudinal position of the cross-correlation phase calculation in the second frame image.
  • an embodiment of the present invention further provides a displacement detecting apparatus for elastic imaging, which includes:
  • a target point obtaining device 501 configured to acquire a target point
  • a searching device 503 configured to acquire a cross-correlation phase calculation position of the target point in the second frame image
  • a cross-correlation phase calculation device 505 configured to calculate a cross-correlation phase according to the cross-correlation phase calculation position
  • a longitudinal displacement result calculating means 507 configured to calculate a longitudinal displacement result according to the cross-correlation phase
  • the strain result calculating means 509 is configured to obtain a strain result by obtaining a gradient of the displacement result.
  • the method further comprises:
  • a baseband signal acquiring device configured to acquire baseband signal data
  • a meshing device configured to divide the mesh in the baseband signal data, wherein the node of the mesh is a displacement detection estimation point, and the target point is from the displacement detection estimation point.
  • the cross-correlation phase calculation position of the acquisition target point in the second frame image includes a data longitudinal position calculated by the cross-correlation phase
  • the searching device is specifically configured to:
  • the longitudinal position of the point is the longitudinal position of the cross-correlation phase calculation of the target point in the second frame image.
  • the cross-correlation phase calculation position of the acquisition target point in the second frame image further includes a data lateral position calculated by the cross-correlation phase
  • the search device is also used to:
  • a point at which the data correlation is most correlated is searched in the vicinity of the longitudinal position of the cross-correlation phase calculation, and the lateral position of the point is the data lateral position calculated by the cross-correlation phase.
  • the search means calculates a point at a certain interval every time when the point at which the data correlation is greatest is searched laterally near the longitudinal position of the cross-correlation phase calculation.
  • the cross-correlation phase calculation device is specifically configured to:
  • the first core and the second core are obtained by respectively calculating the position of the target point in the position of the first frame and the cross-correlation phase in the second frame as the center position.
  • the cross correlation phase is:
  • T 2 is the signal period, which corresponds to the angular frequency of the signal center; it is the result of the longitudinal displacement of the calculated point.
  • the longitudinal displacement result calculation device is specifically configured to:
  • the method further comprises:
  • a calibration device configured to perform error correction on the strain result
  • a smoothing processing device configured to perform spatial smoothing on the strain result
  • Mapping means for mapping the strain results using different grayscale or color maps The elastic imaging method and apparatus according to the embodiment of the present invention acquires the I/Q two-channel echo baseband signal after two frames of compression before and after compression, and adopts guided phase estimation (GPZE, guided phase zero estimation) to quickly detect two.
  • GPZE guided phase estimation
  • the displacement information between the frames, and the axial gradient calculation to obtain the strain information not only can obtain the higher quality strain image, but also greatly reduce the calculation amount and meet the clinical real-time requirements.

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • General Health & Medical Sciences (AREA)
  • Pathology (AREA)
  • Medical Informatics (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Radiology & Medical Imaging (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Theoretical Computer Science (AREA)
  • Surgery (AREA)
  • Immunology (AREA)
  • Biomedical Technology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Molecular Biology (AREA)
  • Biophysics (AREA)
  • Animal Behavior & Ethology (AREA)
  • Public Health (AREA)
  • Veterinary Medicine (AREA)
  • Acoustics & Sound (AREA)
  • Chemical & Material Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Biochemistry (AREA)
  • Quality & Reliability (AREA)
  • Geometry (AREA)
  • Ultra Sonic Daignosis Equipment (AREA)
  • Investigating Or Analyzing Materials By The Use Of Ultrasonic Waves (AREA)

Abstract

一种弹性成像中的位移检测方法和装置，该方法包括：获取目标点；获取目标点在第二帧图像中的互相关相位计算位置；根据所述互相关相位计算位置计算互相关相位；根据所述互相关相位计算纵向位移结果；对所述位移结果求梯度得到应变结果。通过该弹性成像方法和装置，获取压缩前后两帧降采样后的I/Q两路回波基带信号采用引导式相位估计快速检测两帧之间位移信息，再进行轴向梯度计算得到应变信息，不仅可获得较高质量的应变图像，而且大大降低了计算量，满足临床实时要求。

Description

一种弹性成像中的位移检测方法 ^置
技术领域
本发明涉及弹性成像技术， 尤其涉及一种弹性成像中的位移检测方法 及装置。 背景技术
超声弹性成像作为癌症检测， 尤其是乳腺癌良性恶性判别中， 对 B模 式声像图检测的重要辅助手段， 快速应用于临床。 超声弹性成像主要是通 过超声成像手段， 获取目标组织的超声回波信息， 再通过特定的算法检测 出组织弹性信息， 并以图像形式直观显示出来， 以辅助医生诊断或治疗。
传统的超声弹性成像方法需要探头轻微压缩组织或者借助人体自身的 呼吸、 血管搏动等过程， 获取先后两帧超声回波信号， 然后通过特定的位 移检测方法， 获得两帧信号之间的位移 （displacement ), 即为目标组织在 两个不同时刻空间位置变化信息， 通过对位移求轴向 （axial )梯度， 即可 得到组织的轴向应变（strain )信息。 该应变信息即可反映组织的弹性。 在 相同外力压缩下， 应变越大， 表示组织越软， 应变越小， 则表示组织越硬。 将目标组织区域的应变信息以图像形式表现出来， 可直观反映不同组织间 的软硬差别或弹性差别， 即为应变图像（strain image )。 这种方式又叫做 应变成像 ( strain imaging )。
在上述处理过程中， 位移检测是否准确、 计算是否快速， 共同影响着 最终应变图像的对比度噪声比 (contrast-noise-ratio, CNR) , 是否能实时实 现、 帧率是否满足临床需求等。 传统的位移检测方法一般分为三类： 基于 互相关 ( cross-correlation ) 的时域检测方法、 基于相移 ( phase-shift ) 的 方法、 以及基于组织多普勒（TDI )的方法。 基于互相关的方法所得图像质 量好， 但计算量较大， 且对回波数据采样率有一定要求； 基于相移的方法 计算快速， 但存在相位混叠问题， 主要适用于位移较小的情况， 且通常忽 略其侧向（lateral )位移量， 图像质量相对较差； 基于组织多普勒的方法实 现最筒单， 但依赖于多普勒（Doppler )信号的信噪比和处理方法， 同样存 在混叠问题， 且忽略其侧向 （lateral )位移量， 图像质量相对较差。 发明内容
本发明要解决的主要技术问题是， 提供一种弹性成像中的位移检测方 法及装置， 可以降低成像过程的计算量。
为解决上述技术问题， 本发明提供了一种弹性成像中的位移检测方法， 包括：
获取目标点；
获取目标点在第二帧图像中的互相关相位计算位置；
根据所述互相关相位计算位置计算互相关相位；
根据所述互相关相位计算纵向位移结果；
对所述位移结果求梯度得到应变结果。
本发明还提供了一种弹性成像中的位移检测装置， 包括：
目标点获取装置， 用于获取目标点；
搜索装置， 用于获取目标点在第二帧图像中的互相关相位计算位置； 互相关相位计算装置， 用于根据所述互相关相位计算位置计算互相关 相位；
纵向位移结果计算装置， 用于根据所述互相关相位计算纵向位移结果； 应变结果计算装置， 用于对所述位移结果求梯度得到应变结果。
本发明的有益效果是： 通过本发明实施例提出的弹性成像方法和装置， 获取压缩前后两帧降采样后的 I/Q 两路回波基带信号采用引导式相位估计 快速检测两帧之间位移信息， 再进行轴向梯度计算得到应变信息， 不仅可 获得较高质量的应变图像， 而且大大降低了计算量， 满足临床实时要求。 附图说明
图 1为本发明一种弹性成像系统的一实施例的筒图；
图 2为本发明一种弹性成像中的位移检测方法的一实施例的流程图； 图 3为本发明实施例中的帧数据网格化示意图；
图 4为本发明实施例中的位移搜索策略示意图；
图 5为本发明一种弹性成像中的位移检测装置的一实施例的模块图。 具体实施方式
下面通过具体实施方式结合附图对本发明作进一步详细说明。
如图 1 为本发明一种弹性成像系统的一实施例的筒图。 弹性成像模式 下， 超声探头以系统预先设定好的扫描规则进行超声发射并接收回波信息， 经过波束合成环节后输出射频（ radio frequency, RF)信号， 再经过正交解 调阶段， 生成 I/Q两路基带信号， 同时进行降采样过程， 使得 I/Q信号的采 样率下降， 降采样率由系统预先设定， 接着经过引导零相位估计（GPZE, guided phase zero estimation )算法位移检测环节， 每次利用一对 I/Q信 号帧计算出一帧位移结果， 然后经过应变计算环节， 基于位移数据计算得 到目标组织的应变信号。 此后， 再经过后处理环节对其进行错误校正、 平 滑等操作， 最后显示输出一幅应变图像。
GPZE 算法主要通过引导搜索的思维来更准确快速的检测出前后两帧 信号间的互相关函数的相位， 并推导该相位与纵向位移量之间的对应关系， 计算得到两帧信号间的纵向位移量， 在大大降低了计算量的同时保证了位 移估计的质量。 此外， 位移检测的范围也得以拓展， GPZE 算法不仅适用 于小位移情况， 也适用于大位移的情况。
假设前后两帧信号分别表示为：
fu(t,x) = Au(t,x)ei(ev+9>
fc (t, x) = A (t-uy,x- χ)βηωΛ'-^)+θ] =A (t-uy,x- Ux)eJi c 其中， 是信号中心频率， 为信号初始相位， 为信号周期， 与中心 频率 相对应。 为压缩导致的空间横向偏移量， 为压缩所导致的纵向 偏移量。 总是可表示为 =" /2的形式， 其中 n为整数， r总是分布于
_7：/2~7 /2范围内。 经过正交解调以后成为基带信号， 基带信号的复数形式可表示为：
Figure imgf000005_0001
fcb (t,x) = Ac(t-uy,x- Ux、e]'- -W、 于是， 基带信号之间的互相关函数表示为： ("。， )= j fub (t,x)- fcb*(t + u0,x + x0)dt
t0-Ai -Ac(t-uy+u0,x-ux+x0)e-jla,A-nT Z-I>+e]dt χ) -Ac(t -uy + uQ,x-ux + xQ)dt
Figure imgf000006_0001
其中， "。和 代表计算互相关函数时的两帧目标数据之间的纵向与横向 偏移量， R 是两帧数据的包络信号之间的互相关函数。
将上述等式两边同时乘以 则有：
R'(u0,x0) = Rb(u0,x0)e—Jm
= (- i) (w 0)
Figure imgf000006_0002
于是， 只要找到令包络互相关函数 K"。，X 最大时的"。和 。， 即找到两 帧包络数据之间的纵向、 横向偏移量， 即可利用前后两帧相应偏移位置的 基带数据计算出互相关函数 ("。， x。）， 由于互相关函数与 R'函数之间只相差 关系，从而可以进一步计算得到上述互相关相位 = 最终得到两帧 信号间的时移量或纵向位移量表达为：
uy = nTc 12 + τ = nT c 12 + φ I c 注意： 由于时移量（时间单位） 与纵向位移量（长度单位）是一一对 应的， 所以本发明实施例的描述中， "纵向位移量" 也即表达是 "时移量" 的意思， 不刻意加以区别。
本发明所提供的一种弹性成像中的位移检测方法， 包括：
获取目标点；
获取目标点在第二帧图像中的互相关相位计算位置；
根据互相关相位计算位置计算互相关相位；
根据互相关相位计算纵向位移结果。
具体的， 在一个实施例中， 获取目标点之前， 可获取基带信号数据， 在该基带信号数据划分网格， 网格的节点为位移检测估计点， 而目标点则 来自位移检测估计点。 在一个实施例中， 获取目标点在第二帧图片中的互相关相位计算位置 可包括互相关相位计算的数据纵向位置。 则获取目标点在第二帧图像中的 互相关相位计算位置的步骤包括： 获取目标点的较前计算点的纵向位移结 果； 在第二帧图像中， 在以目标点的位置与较前计算点的纵向位移结果的 和为中心的区域搜索与该核数据相关性最大的点， 该点的纵向位置即为目 标点在第二帧图像中的互相关相位计算的纵向位置。
在一个实施例中， 获取目标点在第二帧图片中的互相关相位计算位置 还可包括互相关相位计算的数据横向位置。 则获取目标点在第二帧图像中 的互相关相位计算位置的步骤包括： 在互相关相位计算的纵向位置附近横 向搜索数据相关性最大的点， 该点的横向位置为互相关相位计算的数据横 向位置。
在一个实施例中， 在互相关相位计算的纵向位置附近横向搜索数据相 关性最大的点的步骤是每间隔一定数量的点计算一次。
在另一个实施例中， 获取目标点在第二帧图像中的互相关相位计算位 置的步骤中获取的位置包括互相关相位计算的数据纵向位置与互相关相位 计算的数据横向位置； 则获取目标点在第二帧图像中的互相关相位计算位 置的步骤包括：
在第二帧图像中， 以目标点为中心， 在预置的横向搜索范围内搜索两 帧包络数据互相关最大的点， 该点的横向位置为互相关计算的横向位置； 在第二帧图像中， 以目标点的位置与目标点的较前计算点的纵向位移 结果的和为中心的区域搜索与该核数据相关性最大的点， 该点的纵向位置 即为目标点在第二帧图像中的互相关相位计算的纵向位置。
如图 2所示为本发明一种弹性成像中的位移检测方法的一实施例的流 程图， 包括：
201、 获取数据；
具体地， 获取一对 I/Q基带信号帧数据。
每一帧位移数据计算需要利用两帧 I/Q基带信号数据，所得位移数据指 两帧信号之间的空间相对位移。 上述两帧基带信号数据， 可以是连续两帧 数据， 可以是有一定帧间隔的两帧数据， 帧间隔数量由系统预先设定。 使 用一定间隔的两帧数据， 可以有效调整用于计算的两帧数据间的位移量大 小， 使得最终获得的应变图像质量更好。
每一帧基带信号数据， 都分为 I、 Q两路数据。 基带信号的采样率可能 与 RF数据的采样率不同， 降采样率由系统预先设定。 降采样率越大， 后面 的计算量会越小， 但会影响位移估计的质量以及最终图像的空间分辨率。
得到 I、 Q两路数据后， 还需要进一步计算出其包络数据， 以用于后面 包络偏移量的搜索计算。 包络的计算对 I/Q数据的每个采样点都要进行，采 样点' '处的包络数据计算为：
envelop(i) =
Figure imgf000008_0001
203、 划分位移估计点网格；
由于基带信号的采样率仍然较高， 而位移一般较小， 相邻采样点之间 的位移差别非常微小。 因此， 位移估计不一定要对基带信号数据的每个采 样点位置进行计算， 预先划分好位移检测估计点的位置， 可以有效避免或 减少冗余计算量。
如图 3所示为本发明实施例中的帧数据网格化示意图， 黑点 （即网格 节点处） 即为可能进行位移估计的位置， 图中黑线表示基带信号数据或包 络数据， 网格划分以两帧信号中第一帧的数据采样点位置为基准。
纵向， 从最浅深度的数据开始， 每隔一定数量的数据采样点 （或者每 隔一定深度）后取点， 该点所在位置需要进行位移估计， 该纵向间隔数量 由系统预先设定。
横向， 从探头中心扫描数据线开始， 每隔一定数量的采样线 （或者每 隔一定宽度）后取点， 该点所在位置需要进行位移估计， 该横向间隔数量 由系统预先设定。
本发明中偏移量搜索、 相移检测、 位移计算、 应变估计等过程， 优选 地， 都是以上述网格中的各节点 （即位移检测估计点） 为目标来进行。
划分网格之后， 获取应变图像所需的计算量将会大大减少， 特别是当 横向和纵向间隔很大的时候。 但是间隔太大对成像质量有一定的影响， 会 影响位移估计结果以及最终图像的空间分辨力。
205、 寻找互相关相位计算的数据纵向位置；
为了获得互相关相位计算位置， 所以要寻找目标点在第二帧图像中的 互相关相位计算的数据纵向位置， 即就是确定互相关相位计算的数据纵向 偏移量， 也就是找到令两帧数据之间的包络互相关函数最大时的纵向偏移 量。
该纵向偏移量的搜索以块匹配（block-matching )的思路为基础， 针对 包络数据进行搜索。 每次搜索以包络数据网格划分后的第一帧中某个节点 位置为目标点， 以此为中心选择核数据， 在第二帧数据中沿纵向搜索与该 核数据相关性最大的位置， 该位置与目标点之间的偏移量就是所求的纵向 偏移量。搜索中相关性最大的判别依据可以采用 SAD法、 NCC法等。 SAD 最小的位置或者 NCC最大的位置即为相关性最大的位置。也可以采用其他 类似的判别依据。
搜索时， 本发明在 block-matching的基础上还加入了引导的思维， 以 提高计算速度。 在用户选定了感兴趣的区域（ ROI , region of interest )之 后， 引导搜索又可以有两种不同的实现方法， 区别在于数据计算开始的深 度不同， 相应的计算量也有所不同
方法一： 无论用户选定的 ROI ( region of interest )深度是多少， 计算 总是从探头表面的数据 (即深度为 0或者最接近 0的节点）开始， 每次以 上一个深度的节点的纵向位移结果作为当前节点的纵向数据位置偏移量的 搜索初始值， 在这个初始值的基础上设定一个小的纵向搜索范围进行搜索， 从而得到令上述包络互相关函数最大的纵向偏移量。
对于最浅深度的节点， 由于实际上该节点没有 "上一个深度的节点" 了， 系统设置其对应的 "上一深度节点" 的纵向位移结果固定为 0。
对于其他深度的节点， 假设当前节点深度位置为 而上一个深度的节 点所得最终纵向位移结果为 , 则当前节点对应的包络纵向偏移量搜索初 始位置是" ， 即需要在第二帧包络数据以深度 + 为中心设置一个小的 纵向搜索范围， 找到该范围中令两者相关性最大的位置， 即为当前节点对 应的包络纵向偏移量。 该纵向搜索范围由系统预先设定， 范围越小， 计算 量越小。
方法二： 计算总是从用户选定的 R0I中的最浅深度的节点开始， 每次 以上一个深度的节点的位移结果作为当前深度的节点的纵向数据位置偏移 量的搜索初始值， 在这个初始值的基础上设定一个小的纵向搜索范围进行 搜索， 从而得到令上述包络互相关函数最大的纵向偏移量。 对于 ROI中最浅深度的节点， 由于其 "上一深度的节点" 的最终纵向 位移结果不可知， 因此需要特殊处理。 假设其深度位置为 搜索时需要在 第二帧数据中纵向位置 上下设置一个较大的纵向搜索范围，找到相关性最 大的位置作为该节点对应的包络纵向偏移量。 该纵向搜索范围由系统预先 设定， 一般需要比有引导搜索时 （即有搜索初始值时） 的搜索范围设置得 更大， 如果设置太小则可能导致搜索结果不准确。
对于 ROI中其他深度的节点， 假设上一个深度的节点所得位移结果为 uy , 当前节点深度位置为 则当前节点对应的包络纵向偏移量搜索初始位 置是 + 即需要在第二帧包络数据以深度 + 为中心设置一个小的纵向 搜索范围， 找到该范围中令两者相关性最大的位置， 即为当前节点对应的 包络纵向偏移量。 该纵向搜索范围由系统预先设定， 范围越小， 计算量越 小。 假设当前节点所得包络纵向偏移量结果为"。， 则当前节点在互相关相 位计算时对应的第一帧和第二帧中的数据纵向位置分别是和 ^ + "。。
207、 寻找互相关相位计算的数据横向位置；
确定好纵向位置后， 优选地， 还可以继续寻找横向位置， 即寻找令两 帧数据包络互相关函数最大时的横向偏移量。
横向位置的搜索在上述纵向位置的结果上引导进行， 即得到上述纵向 位置后， 在第二帧中找到该纵向位置， 以该纵向位置为中心， 设置一定的 横向搜索范围 （或叫搜索区域）， 在该范围内两帧包络数据互相关最大的位 置相对于原位置的横向偏移量即为所求。 最终互相关相位计算时第二帧数 据的横向位置即根据该横向偏移量来选取。 该横向搜索范围由系统预先设 定。 搜索范围越小， 计算量越小。 搜索中的相关性判别依据类似纵向位置 的搜索， 可以采用 SAD、 NCC或其他判别方法。
假设第一帧数据中当前目标节点的二维位置为 , 其中 是时移量， ^是扫描线位置。 设上一步中所得包络纵向偏移量为"。， 则搜索时， 以""。 位置为中心设置一个小的横向搜索区域， 寻找该区域中两帧包络数据互相 关最大的位置， 其相对于原位置的偏移量即为所求横向偏移量。 假设该横 向偏移量结果为 。， 则当前节点在互相关相位计算时对应的第一帧与第二 帧的数据位置分别是 ^ 和 ^ + "。， )。
由于组织的横向位移相对较小， 随着深度的变化其改变也较小， 因此， 对上述横向位置的搜索操作不需要沿着每个深度的节点连续进行， 可以每 间隔几个深度的节点后再横向搜索以更新横向位置。 这样可以在保证图像 质量不发生明显变化的情况下进一步减少计算量。 该横向位置更新的深度 间隔的大小由系统预先设定， 间隔太小则更新较频繁， 可能会增加冗余计 算量， 间隔太大则更新较少， 可能会影响互相关相位估计质量。
还可以根本不进行横向位置的搜索， 只进行纵向位置的搜索， 横向位 置默认不变， 则计算会进一步筒化。 当组织实际横向位移不大时， 对互相 关相位估计质量影响不大； 不过当组织实际横向位移较大时， 估计质量会 有所降低。
209、 互相关相位计算；
得到互相关相位计算的数据位置后， 可以进行互相关相位估计。 互相 关相位计算利用前后两帧 I/Q数据进行。 计算时， 以该位置为中心（或者为 基准）， 取出附近或周围一块数据进行计算， 如图 4所示。 该块数据的大小 由系统预先设定， 数据块的大小对互相关相位估计结果有影响。
参与计算的数据块称作核（kernel ), 可以是一维的， 也可以是二维的。 在图 4 中， 假设当前目标点 （即当前节点）在第一帧数据里的位置坐标是 (n) , 而前两步搜索出来的包络纵向和横向偏移量分别为"。和 。， 则其对应 第二帧数据里的位置坐标就是 + "。，^+ ) ,图 4中分别以 (Π)与 + "。Ή ) 为中心或基准点向纵向和横向附近延伸一定距离后得到图中的区域， 这个 区域内的数据就是核数据。 目标点在两帧图像中产生核的方法有多种， 可 以以该点作为基准点向各方向延伸一定的点数得到， 也可以用其他方式得 到， 但是在两帧图像中产生核的方法是一样的， 所以两帧图像中的核数据 大小相同， 用于计算的每个点都在另一帧图像中的核区域中有对应的点。
相位计算使用 I、 Q数据。 将第一帧和第二帧中相应位置的 I、 Q核数 据取出， 则求得相位为：
n y^ 2、
φ - arctan[(-l) 其中, l与 βΐ是第一帧的数据， 与 β2是第二帧的数据， 表示核中各 个数据采样点在核中的相对坐标或位置，同一 在两帧数据中表示的相对 位置相同， 例如假设（0, 0 )在第一帧的核数据中表示核数据左下角的点， 那么 （0, 0 )在第二帧的核数据中也表示核数据左下角的点。
上式中， "的计算方法与上一深度节点的最终位移估计结果有关，假设 上一深度节点的最终纵向位移结果为 ， 贝' J :
u
n = roundf ~―)
T /2 其中， round表示四舍五入取整（也可以采用向上取整或向下取整的方 式）。
π π
此外， 由于 arctan函数计算的结果范围在— ϊ ~ ϊ之间， 还需要根据其上 式
Figure imgf000012_0001
中分子分母的正负情况将结果校正到
- Γ〜 Γ范围中去。按三角函数的相位分布规律，上述分子的符号与 si 对应， 分母的符号与 ec)s^对应。
211、 计算纵向位移结果； 经过上述计算后， 对于当前节点 （或位移估计点）， 最终两帧之间的纵 向位移结果为：
Figure imgf000012_0002
其中， 为信号周期， 与信号中心角频率 相对应。 "^ / 2项的引入， 补偿了相位计算的混叠性， 从而使得本算法即适用于小位移的情况， 也适 用于大位移的情况。 上述纵向位移均以采样时间为单位的方式来表示， 也可以换算成物理 长度单位来表示， 二者是一一对应的。
得到网格中所有节点 （或位移估计点）处的纵向位移结果后， 对位移 数据沿着纵向求梯度， 即可得到应变结果， 即应变值。
在应变后处理中， 可以对所得的应变结果进行一定的错误校正， 比如 检测其中的异常跳变点或者明显错误点进行校正； 还可以对其进行空间平 滑， 以改善图像显示效果； 或者对其使用不同的灰阶或彩色图谱进行映射， 以增强图像对比度。 还可以进行其他增加图像质量的操作。
最终， 将所需区域的应变结果进行输出显示成为应变图像， 即可反映 该区域组织弹性差异。
GPZE位移检测算法， 一方面由于使用了上一深度的位移结果来引导 下一深度的位移计算， 减少了搜索计算量； 另一方面使用相位估计的方法 来计算位移， 对原始数据采样率的要求不高， 从而大大降低了计算量。
同时， 还可以引入横向搜索， 提高了纵向位移量的估计质量， 使得最 终的图像信噪比更高。
此外， 还补偿了相位混叠带来的干扰， 从而扩展了位移估计的适用范 围。
在上述实施例中， 首先寻找互相关相位计算的数据纵向位置， 然后再 寻找横向位置， 实际上， 也可以先寻找横向位置， 然后再寻找纵向位置； 或者， 上述对横向位置与纵向位置的搜索可以无确定的时间先后关系， 可 以分别独立进行； 寻找的方法同上述实施例的相应的方法。 具体地， 例如 横向位置的搜索方法可以是以目标点为中心， 设置一定的横向搜索范围， 在该范围内两帧包络数据互相关最大的位置相对于原位置的横向偏移量即 为所求。 纵向位置的寻找方法例如可以是： 以目标点的位置与所述较前计 算点的纵向位移结果的和为中心的区域搜索与该核数据相关性最大的点， 该点的纵向位置即为目标点在第二帧图像中的互相关相位计算的纵向位 置。
如图 5所示， 本发明实施例还提出了一种弹性成像中的位移检测装置， 包括：
目标点获取装置 501 , 用于获取目标点；
搜索装置 503, 用于获取目标点在第二帧图像中的互相关相位计算位 置；
互相关相位计算装置 505, 用于根据所述互相关相位计算位置计算互 相关相位；
纵向位移结果计算装置 507, 用于根据所述互相关相位计算纵向位移 结果；
应变结果计算装置 509, 用于对所述位移结果求梯度得到应变结果。 优选地， 还包括：
基带信号获取装置， 用于获取基带信号数据；
网格划分装置， 用于在所述基带信号数据划分网格， 网格的节点为位 移检测估计点， 所述目标点来自位移检测估计点。
优选地， 所述获取目标点在第二帧图像中的互相关相位计算位置包括 互相关相位计算的数据纵向位置；
所述搜索装置具体用于：
获取所述目标点的较前计算点的纵向位移结果； 在第二帧图像中， 在 以所述目标点的位置与所述较前计算点的纵向位移结果的和为中心的区域 搜索与该核数据相关性最大的点， 该点的纵向位置即为目标点在第二帧图 像中的互相关相位计算的纵向位置。
优选地， 所述获取目标点在第二帧图像中的互相关相位计算位置还包 括互相关相位计算的数据横向位置；
所述搜索装置还用于：
在所述互相关相位计算的纵向位置附近横向搜索数据相关性最大的 点， 所述点的横向位置为互相关相位计算的数据横向位置。
优选地， 所述搜索装置在互相关相位计算的纵向位置附近横向搜索数 据相关性最大的点时， 每间隔一定数量的点计算一次。
优选地， 所述互相关相位计算装置具体用于：
以目标点在第一帧的位置以及第二帧中的互相关相位计算位置分别为 中心位置得到第一核与第二核，
所述互相关相位 为：
φ = arctan[( -1)" .∑W，潭
rtil (i ) (. i) + Ql (i )Q1(.i ) 其中, l与 βΐ是第一帧的数据， 与 β2是第二帧的数据， 表示核中各 个数据采样点在核中的相对位置；
其中， n = round( y
T 2 其中， 为信号周期， 与信号中心角频率 相对应； 为较前计算点 的纵向位移结果。
优选地， 所述纵向位移结果计算装置具体用于：
、上
计弄： u y' = ηΤ c 12 + φ r Ι ω c 其中， "'为所述纵向位移结果， ^为所述互相关相位；
优选地， 还包括：
校正装置， 用于对所述应变结果进行错误校正；
或，
平滑处理装置， 用于对所述应变结果进行空间平滑处理；
或，
映射装置， 用于对所述应变结果使用不同的灰阶或彩色图谱进行映射。 本发明实施例提出的弹性成像方法和装置， 获取压缩前后两帧降采样 后的 I/Q 两路回波基带信号 （baseband signal ) 采用引导式相位估计 ( GPZE, guided phase zero estimation )快速检测两帧之间位移信息， 再 进行轴向梯度计算得到应变信息， 不仅可获得较高质量的应变图像， 而且 大大降低了计算量， 满足临床实时要求。
以上内容是结合具体的实施方式对本发明所作的进一步详细说明， 不能认 定本发明的具体实施只局限于这些说明。 对于本发明所属技术领域的普通 技术人员来说， 在不脱离本发明构思的前提下， 还可以做出若干筒单推演 或替换， 都应当视为属于本发明的保护范围。

Claims

权利要求书
1、 一种弹性成像中的位移检测方法， 其特征在于， 包括：
获取目标点；
获取目标点在第二帧图像中的互相关相位计算位置；
根据所述互相关相位计算位置计算互相关相位；
根据所述互相关相位计算纵向位移结果。
2、 如权利要求 1所述的弹性成像中的位移检测方法， 其特征在于： 获取目标点之前还包括：
获取基带信号数据；
在所述基带信号数据划分网格， 网格的节点为位移检测估计点；
所述目标点来自位移检测估计点。
3、 如权利要求 1所述的弹性成像中的位移检测方法， 其特征在于： 所述获取目标点在第二帧图像中的互相关相位计算位置的步骤中获取的位 置包括互相关相位计算的数据纵向位置；
所述获取目标点在第二帧图像中的互相关相位计算位置的步骤包括： 获取所述目标点的较前计算点的纵向位移结果；
在第二帧图像中， 在以所述目标点的位置与所述较前计算点的纵向位移结 果的和为中心的区域搜索与该核数据相关性最大的点， 该点的纵向位置即 为目标点在第二帧图像中的互相关相位计算的纵向位置。
4、 如权利要求 3所述的弹性成像中的位移检测方法， 其特征在于： 所述获取目标点在第二帧图像中的互相关相位计算位置的步骤中获取的位 置还包括互相关相位计算的数据横向位置；
所述获取目标点在第二帧图像中的互相关相位计算位置的步骤还包括： 在所述互相关相位计算的纵向位置附近横向搜索数据相关性最大的点， 所 述点的横向位置为互相关相位计算的数据横向位置。
5、 如权利要求 4所述的弹性成像中的位移检测方法， 其特征在于： 所述在互相关相位计算的纵向位置附近横向搜索数据相关性最大的点的步 骤是每间隔一定数量的点计算一次。
6、 如权利要求 1所述的弹性成像中的位移检测方法， 其特征在于： 所述获取目标点在第二帧图像中的互相关相位计算位置的步骤中获取的位 置包括互相关相位计算的数据纵向位置与互相关相位计算的数据横向位 置；
所述获取目标点在第二帧图像中的互相关相位计算位置的步骤包括： 在第二帧图像中， 以目标点为中心， 在预置的横向搜索范围内搜索两帧包 络数据互相关最大的点， 该点的横向位置为互相关计算的横向位置； 在第二帧图像中， 以目标点的位置与目标点的较前计算点的纵向位移结果 的和为中心的区域搜索与该核数据相关性最大的点， 该点的纵向位置即为 目标点在第二帧图像中的互相关相位计算的纵向位置。
7、 如权利要求 1所述的弹性成像中的位移检测方法， 其特征在于： 根据所述互相关相位计算位置计算互相关相位的步骤包括：
以目标点在第一帧的位置以及第二帧中的互相关相位计算位置分别为中心 位置得到第一核与第二核，
所述互相关相位 为： φ = arctan[( -1)" .∑W，潭 其中， ι与 βι是第一帧的数据， 与 是第二帧的数据， （^')表示核中各个数 据采样点在核中的相对位置；
其中， n
Figure imgf000017_0001
其中， 为信号周期， 与信号中心角频率 相对应； ^为较前计算点的纵向 位移结果。
8、 如权利要求 7所述的弹性成像中的位移检测方法， 其特征在于： 根据所述互相关相位计算纵向位移结果的步骤包括： 计算： uy， = nTc l 1 + φ I coc
其中， " 'y为所述纵向位移结果， ^为所述互相关相位；
9、如权利要求 1所述的弹性成像中的位移检测方法，其特征在于，还包括： 对所述应变结果进行错误校正；
或，
对所述应变结果进行空间平滑处理；
或，
对所述应变结果使用不同的灰阶或彩色图谱进行映射。
10、 一种弹性成像中的位移检测装置， 其特征在于， 包括：
目标点获取装置， 用于获取目标点；
搜索装置， 用于获取目标点在第二帧图像中的互相关相位计算位置； 互相关相位计算装置， 用于根据所述互相关相位计算位置计算互相关相位； 纵向位移结果计算装置， 用于根据所述互相关相位计算纵向位移结果； 应变结果计算装置， 用于对所述位移结果求梯度得到应变结果。
11、 如权利要求 10所述的弹性成像中的位移检测装置， 其特征在于， 还包 括：
基带信号获取装置， 用于获取基带信号数据；
网格划分装置， 用于在所述基带信号数据划分网格， 网格的节点为位移检 测估计点， 所述目标点来自位移检测估计点。
12、 如权利要求 1 0所述的弹性成像中的位移检测装置， 其特征在于： 所述获取目标点在第二帧图像中的互相关相位计算位置包括互相关相位计 算的数据纵向位置；
所述搜索装置用于：
获取所述目标点的较前计算点的纵向位移结果； 在第二帧图像中， 在以所 述目标点的位置与所述较前计算点的纵向位移结果的和为中心的区域搜索 与该核数据相关性最大的点， 该点的纵向位置即为目标点在第二帧图像中 的互相关相位计算的纵向位置。
13、 如权利要求 12所述的弹性成像中的位移检测装置， 其特征在于： 所述获取目标点在第二帧图像中的互相关相位计算位置还包括互相关相位 计算的数据横向位置；
所述搜索装置还用于：
在所述互相关相位计算的纵向位置附近横向搜索数据相关性最大的点， 所 述点的横向位置为互相关相位计算的数据横向位置。
14、 如权利要求 13所述的弹性成像中的位移检测装置， 其特征在于： 所述搜索装置在互相关相位计算的纵向位置附近横向搜索数据相关性最大 的点时， 每间隔一定数量的点计算一次。
15、 如权利要求 10所述的弹性成像中的位移检测装置， 其特征在于： 所述获取目标点在第二帧图像中的互相关相位计算位置的步骤中获取的位 置包括互相关相位计算的数据纵向位置与互相关相位计算的数据横向位 置；
所述搜索装置用于：
在第二帧图像中， 以目标点为中心， 在预置的横向搜索范围内搜索两帧包 络数据互相关最大的点， 该点的横向位置为互相关计算的横向位置； 在第二帧图像中， 以目标点的位置与目标点的较前计算点的纵向位移结果 的和为中心的区域搜索与该核数据相关性最大的点， 该点的纵向位置即为 目标点在第二帧图像中的互相关相位计算的纵向位置。
16、 如权利要求 10所述的弹性成像中的位移检测装置， 其特征在于： 所述互相关相位计算装置用于：
以目标点在第一帧的位置以及第二帧中的互相关相位计算位置分别为中心 位置得到第一核与第二核，
所述互相关相位 为： φ = arctan[(-l)" .∑ '·，巢 -A ( ·)¾ ( ·)] 其中， l与 2l是第一帧的数据， ^与 是第二帧的数据， 表示核中各个数 据采样点在核中的相对位置； 其中'
n = round(
Tc / 2
其中， 为信号周期， 与信号中心角频率 相对应； 为较前计算点的纵 向位移结果。
17、 如权利要求 16所述的弹性成像中的位移检测装置， 其特征在于： 所述纵向位移结果计算装置用于： 计算： uy， = nTc l 1 + φ I coc
其中， 为所述纵向位移结果， ^为所述互相关相位；
18、 如权利要求 10所述的弹性成像中的位移检测装置， 其特征在于， 还包 括：
校正装置， 用于对所述应变结果进行错误校正；
或，
平滑处理装置， 用于对所述应变结果进行空间平滑处理；
或，
映射装置， 用于对所述应变结果使用不同的灰阶或彩色图谱进行映射。
19、 一种超声成像系统， 其特征在于： 包括如权利要求 1 0-18任一项所述 的弹性成像中的位移检测装置。
PCT/CN2012/073135 2011-06-14 2012-03-27 一种弹性成像中的位移检测方法及装置 Ceased WO2012171388A1 (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
US14/126,222 US9607405B2 (en) 2011-06-14 2012-03-27 Method and device for detecting displacement in elastography

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN201110159194.6A CN102824194B (zh) 2011-06-14 2011-06-14 一种弹性成像中的位移检测方法及装置
CN201110159194.6 2011-06-14

Publications (1)

Publication Number Publication Date
WO2012171388A1 true WO2012171388A1 (zh) 2012-12-20

Family

ID=47327648

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2012/073135 Ceased WO2012171388A1 (zh) 2011-06-14 2012-03-27 一种弹性成像中的位移检测方法及装置

Country Status (3)

Country Link
US (1) US9607405B2 (zh)
CN (1) CN102824194B (zh)
WO (1) WO2012171388A1 (zh)

Families Citing this family (15)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103211625B (zh) * 2013-01-11 2015-08-19 深圳市恩普电子技术有限公司 基于弹性成像的生物位移计算方法
CN103735287B (zh) * 2013-12-05 2015-11-18 中国科学院苏州生物医学工程技术研究所 一种血管内超声弹性成像二维多级混合位移估计方法
CN104739442B (zh) * 2013-12-25 2017-06-16 深圳迈瑞生物医疗电子股份有限公司 压力弹性成像位移检测方法、装置和超声成像设备
US9782152B2 (en) * 2014-08-18 2017-10-10 Vanderbilt University Method and system for real-time compression correction for tracked ultrasound and applications of same
CN110432926B (zh) * 2014-09-03 2022-06-07 深圳迈瑞生物医疗电子股份有限公司 弹性测量检测方法及系统
CN105092595B (zh) * 2015-08-31 2018-03-02 哈尔滨工业大学（威海） 应用于钢轨探伤的光声弹性成像方法及装置
CN105232087B (zh) * 2015-11-05 2018-01-09 无锡祥生医疗科技股份有限公司 超声弹性成像实时处理系统
JP6728767B2 (ja) * 2016-02-29 2020-07-22 コニカミノルタ株式会社 超音波診断装置及び超音波情報処理方法
CN107970043B (zh) * 2017-12-28 2021-01-19 深圳开立生物医疗科技股份有限公司 一种剪切波的检测方法及装置
CN109745073B (zh) * 2019-01-10 2021-08-06 武汉中旗生物医疗电子有限公司 弹性成像位移的二维匹配方法及设备
CN110477948B (zh) * 2019-08-21 2022-05-10 东软医疗系统股份有限公司 弹性成像方法及装置、成像设备、存储介质
CN113397588B (zh) * 2020-03-16 2024-08-30 深圳市理邦精密仪器股份有限公司 弹性成像方法、装置及医疗设备
CN111833383B (zh) * 2020-07-27 2021-02-09 西南石油大学 二维联合局部位移拟合的区域增长贝叶斯运动追踪方法
CN114680936A (zh) * 2020-12-25 2022-07-01 深圳迈瑞生物医疗电子股份有限公司 血管超声数据处理方法、设备及存储介质
CN112674799B (zh) * 2021-01-05 2022-11-25 青岛海信医疗设备股份有限公司 超声弹性成像方法、电子设备及存储介质

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1586408A (zh) * 2004-08-20 2005-03-02 清华大学 一种多尺度的生物组织位移估计方法
CN1964670A (zh) * 2004-10-12 2007-05-16 株式会社日立医药 超声波探头以及超声波成像装置
CN101530333A (zh) * 2002-07-31 2009-09-16 株式会社日立医药 超声诊断系统和应变分布及弹性模量分布显示方法
US20090247871A1 (en) * 2008-03-25 2009-10-01 Tomy Varghese Rapid two/three-dimensional sector strain imaging
US20100106018A1 (en) * 2008-10-27 2010-04-29 Jingfeng Jiang Ultrasonic strain imaging device with selectable cost-function

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP3182479B2 (ja) * 1993-08-12 2001-07-03 淑 中山 弾性計測装置
JP4644819B2 (ja) * 2004-04-22 2011-03-09 国立大学法人電気通信大学 微小変位計測法及び装置
US8100831B2 (en) * 2006-11-22 2012-01-24 General Electric Company Direct strain estimator for measuring elastic properties of tissue

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101530333A (zh) * 2002-07-31 2009-09-16 株式会社日立医药 超声诊断系统和应变分布及弹性模量分布显示方法
CN1586408A (zh) * 2004-08-20 2005-03-02 清华大学 一种多尺度的生物组织位移估计方法
CN1964670A (zh) * 2004-10-12 2007-05-16 株式会社日立医药 超声波探头以及超声波成像装置
US20090247871A1 (en) * 2008-03-25 2009-10-01 Tomy Varghese Rapid two/three-dimensional sector strain imaging
US20100106018A1 (en) * 2008-10-27 2010-04-29 Jingfeng Jiang Ultrasonic strain imaging device with selectable cost-function

Also Published As

Publication number Publication date
CN102824194A (zh) 2012-12-19
US20140254869A1 (en) 2014-09-11
CN102824194B (zh) 2016-08-03
US9607405B2 (en) 2017-03-28

Similar Documents

Publication Publication Date Title
WO2012171388A1 (zh) 一种弹性成像中的位移检测方法及装置
JP4932984B2 (ja) 超音波撮像において組織変形の実時間計算および表示を実現する方法
US8265358B2 (en) Ultrasonic image processing apparatus and method for processing ultrasonic image
US9612142B2 (en) Method and system for measuring flow through a heart valve
US7837625B2 (en) Ultrasonic image processor and ultrasonic diagnostic instrument
JP2010525850A (ja) ひずみ画像表示システム
CN104739442B (zh) 压力弹性成像位移检测方法、装置和超声成像设备
JP5865463B1 (ja) 超音波画像処理装置
CN110215233A (zh) 一种基于超声平面波扫描的分段式脉搏波成像方法
JP6515095B2 (ja) 解剖学的にインテリジェントな心エコー検査法における肋骨妨害物描出
WO2010083469A1 (en) Dynamic ultrasound processing using object motion calculation
JP6381972B2 (ja) 医用画像処理装置および医用画像診断装置
CN102824193B (zh) 一种弹性成像中的位移检测方法、装置及系统
EP3508132A1 (en) Ultrasound system and method for correcting motion-induced misalignment in image fusion
EP2924656A1 (en) Diagnostic image generation apparatus and diagnostic image generation method
CN106580372A (zh) 一种超声彩色血流成像的脉冲重复频率调整方法及装置
JP2021525619A (ja) 胎児体重推定を実施するための方法およびシステム
JPH07178086A (ja) 超音波診断装置及び超音波診断方法
US20160063742A1 (en) Method and system for enhanced frame rate upconversion in ultrasound imaging
US20210033440A1 (en) Ultrasonic system for detecting fluid flow in an environment
CN101926657A (zh) 一种超声图像特征追踪方法及其系统
KR101656127B1 (ko) 계측 장치 및 그 제어 프로그램
CN106604683B (zh) 超声波诊断装置
JPH08131436A (ja) 超音波診断装置
US20100081932A1 (en) Ultrasound Volume Data Processing

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: 12800218

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

WWE Wipo information: entry into national phase

Ref document number: 14126222

Country of ref document: US

32PN Ep: public notification in the ep bulletin as address of the adressee cannot be established

Free format text: NOTING OF LOSS OF RIGHTS PURSUANT TO RULE 112(1) EPC (EPO FORM 1205 DATED 07.05.2014)

122 Ep: pct application non-entry in european phase

Ref document number: 12800218

Country of ref document: EP

Kind code of ref document: A1