WO2016106959A1 - 多参数模型指导的迭代提取血管造影图像运动参数的方法 - Google Patents
多参数模型指导的迭代提取血管造影图像运动参数的方法 Download PDFInfo
- Publication number
- WO2016106959A1 WO2016106959A1 PCT/CN2015/072681 CN2015072681W WO2016106959A1 WO 2016106959 A1 WO2016106959 A1 WO 2016106959A1 CN 2015072681 W CN2015072681 W CN 2015072681W WO 2016106959 A1 WO2016106959 A1 WO 2016106959A1
- Authority
- WO
- WIPO (PCT)
- Prior art keywords
- fourier transform
- frequency
- iteration
- module
- point
- 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
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING OR CALCULATING; COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/20—Analysis of motion
-
- 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
- G06—COMPUTING OR CALCULATING; COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/20—Analysis of motion
- G06T7/262—Analysis of motion using transform domain methods, e.g. Fourier domain methods
-
- G—PHYSICS
- G06—COMPUTING OR CALCULATING; COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10116—X-ray image
-
- G—PHYSICS
- G06—COMPUTING OR CALCULATING; COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30004—Biomedical image processing
- G06T2207/30101—Blood vessel; Artery; Vein; Vascular
Definitions
- the invention belongs to the field of intersection of digital signal processing and medical imaging, and more particularly relates to a method for iteratively extracting motion parameters of an angiographic image guided by a multi-parameter model.
- the image of the coronary blood vessel on the contrast surface undergoes a two-dimensional motion due to the influence of respiratory motion.
- the physician may move the catheter bed in order to enable the contrast sequence to cover the entire coronary system, resulting in a two-dimensional translation of the entire frame between different frames of an imaging sequence. Therefore, the coronary angiography image records on the one hand the projection of the heart's own motion on a two-dimensional plane, and also superimposes the two-dimensional motion of the coronary artery on the contrast surface and the two-dimensional translational motion of the patient caused by the respiratory motion of the human body. Then, to reconstruct a three-dimensional cardiovascular tree from a dynamic two-dimensional angiogram, it is necessary to separate the heart, respiratory and translational movements of the human body.
- the external registration method is to set some markers before the angiography, and track them during the angiography, but this marking method is usually invasive; the internal registration method is divided into marking method and segmentation method, etc.
- the marking method is in the figure. Select some anatomical points for registration. However, not every figure has such an anatomical point, and it is also necessary to be familiar with human anatomy.
- the segmentation method is registered by the segmented anatomical line, which is applicable to both the rigid body model and the deformation model, but it can only extract the motion of the anatomical line, and cannot separate the whole rigid translational motion separately.
- the apparent vascular motion in the angiogram contains not only the patient's translation, but also the physiological and respiratory movements of the blood vessel itself.
- Another method is to obtain respiratory and cardiac motion parameters under dual-arm x-ray contrast.
- the method must first reconstruct the three-dimensional motion sequence of the coronary artery and then decompose to obtain respiratory and cardiac motion.
- the processing is very complicated.
- the two different angles of contrast images involved in reconstruction have different respiratory motion phases, so the literature method is not suitable for single-arm imaging systems.
- the contrast time is usually very short, and the length of the contrast image obtained by us is limited, so that the motion of the structural feature points in the contrast image is finitely long.
- This finite-length motion contains periodic components, non-periodic components, and incomplete periodic components, so conventional discrete Fourier transforms cannot be separated.
- the invention provides a method for automatically extracting multiple motion parameters such as translational motion, respiratory motion and cardiac motion from an X-ray contrast sequence image by automatically and continuously optimizing the iteration, by selecting some feature points on the coronary vessels and tracking in the sequence.
- the motion curves of them are obtained with time, and then under the constraint of multi-motion parameter model, the discrete Fourier forward-reverse transform and the estimated reconstructed signal and the original signal have the minimum local mean square error of each frequency point in the frequency domain.
- the motion components are optimized by loop iteration to minimize the mean squared difference between the original signal and the reconstructed signal, and the optimal motion curves such as two-dimensional heart, tremor, respiration and translation are obtained. Has a good Pro Bed suitability.
- the present invention provides a method for iteratively extracting motion parameters of an angiographic image guided by a multi-parameter model, the purpose of which is to solve the problem that the existing method cannot automatically extract the translational motion, and Technical problems such as breathing and heart movement cannot be accurately separated.
- a multi-parameter model-guided iterative method for extracting angiographic image motion parameters including the following steps:
- step (7) determining whether the estimated minimum mean square error is greater than a preset threshold, if it is greater, then proceeding to step (7), otherwise the process ends;
- step (1) adopts the following formula:
- L(n) represents the translational motion
- r i (n) represents the respiratory motion of the i-th feature point
- c i (n) represents the cardiac motion of the i-th feature point
- h i (n) represents the i-th feature
- t i (n) represents the other motion of the ith feature point.
- step (2) adopts the following formula:
- L(k), C(k), R(k), and H(k) represent the respective harmonics of L(n), c(n), r(n), and h(n), respectively.
- the estimated minimum mean square error of the frequency point in step (5) It is calculated using the following formula:
- step (7) comprises the following sub-steps:
- the heart signal spectrum of the j+1th iteration is calculated, and after the discrete Fourier transform is performed in the frequency range of the frequency point, the j+1th optimal cardiac time domain signal is obtained by the following formula.
- ⁇ is the frequency obtained by the discrete Fourier transform of the signal
- M represents the size of the window, generally set to 3;
- the minimum mean square error in step (9) It is calculated using the following formula:
- a multi-parameter model-directed iterative method for extracting angiographic image motion parameters comprising:
- the first module is configured to automatically select one vascular structural feature point from the medical image of the contrast sequence, and automatically track the feature points in the entire angiographic sequence to obtain a tracking sequence of each feature point.
- i (n), i 1, ..., I ⁇ , where n represents the frame number of the medical image in the sequence of contrast pictures:
- a fourth module configured to perform a Fourier transform on the tracking sequence of the frequency point in the amplitude and frequency ranges of each frequency point to obtain a Fourier transform result
- a fifth module configured to perform an inverse Fourier transform on the Fourier transform result obtained by the fourth module, and obtain an estimated minimum mean square error of the frequency point;
- the sixth module is configured to determine whether the estimated minimum mean square error is greater than a preset threshold, and if it is greater, transfer to the seventh module, otherwise the process ends;
- the seventh module is configured to process the frequency spectrum of each frequency point by using a multi-parameter iterative optimization algorithm to obtain a time domain signal after the j+1th iteration;
- the eighth module is configured to process the residual signal by using a translation model to obtain a translation signal after the j+1th iteration;
- the ninth module is configured to superimpose the time domain signal after the j+1th iteration and the translation signal after the j+1th iteration to obtain the estimated mixed signal after the j+1th iteration, and calculate the first formula by using the following formula Minimum mean square error after j+1 iterations:
- the tenth module is configured to determine whether the minimum mean square error after the j+1th iteration is greater than the sixth module.
- the threshold in the middle if it is greater, returns to the seventh module, otherwise the process ends.
- step (2) the method for automatically acquiring structural feature points to obtain a motion curve overcomes the problem of setting markers on the surface of the heart for tracking and acquiring feature points.
- step (7) can solve the problem that the translation movement cannot be automatically extracted in the literature; and solve the problem that other periodic movements (such as respiratory movement, cardiac motion, etc.) cannot be separated;
- step (7) Due to the multiple loop iterative optimization method of step (7), the problem that the motion signals extracted in the patent are inaccurate and the robustness is not solved is solved.
- FIG. 1 is a flow chart of a method for iteratively extracting motion parameters of an angiographic image guided by a multi-parameter model of the present invention
- Figure 2 is an original contrast image obtained at an angle
- FIG. 3 is a feature point selected from the blood vessel skeleton image of FIG. 2;
- 5(a) and (b) are the X-axis component spectrum of the fourth feature point of FIG. 2 and the Y-axis component spectrum of the fourth feature point;
- 6(a) and (b) are the heart motion curves of the X-axis and the Y-axis separated by the respective feature points of FIG. 2;
- Figures 9(a) and (b) show the X-axis and Y-axis translation curves separated by the feature points of Figure 2 and A comparison of the translation curves obtained by manually tracking the bone points.
- Figure 11 is an original contrast image of another angle obtained
- Figure 12 is a feature point of the blood vessel skeleton image selected in Figure 11;
- 16(a) and (b) are respiratory motion curves of the X-axis and the Y-axis separated by the respective feature points of FIG. 11;
- 17(a) and (b) are residual curves of the X-axis and Y-axis separation separated by the respective feature points of Fig. 11.
- Minimum mean square error that is, the mean squared summation of the estimated signal and the original signal difference is the smallest.
- the multi-iterative parametric model takes advantage of the properties of discrete Fourier transform, and the method of minimum mean square error is used to continuously optimize the discrete signals in the time domain and frequency domain to minimize the residual.
- the present invention proposes multiple motion parameters of structural feature points.
- Number model The vascular structural feature points are automatically selected on the coronary angiogram, and they are automatically tracked in the contrast sequence to obtain their displacement curves with time.
- the reconstructed signal and the original signal are estimated based on the minimum initial value of the local mean square error of each frequency point in the frequency domain.
- the loop iteration is used to make the difference signal between the original signal and the reconstructed signal.
- the smallest global optimization results in a final estimated translational motion curve, respiratory motion profile, and cardiac motion parameters.
- the method has wide applicability and flexibility, and can be applied to almost all contrast sequence images (of course, it needs to have a clear blood vessel distribution and contains two or more cardiac cycles, which is not difficult to do), including X-ray contrast images of the arms.
- the method for iteratively extracting angiographic image motion parameters guided by the multi-parameter model of the present invention includes the following steps:
- L(n) represents the translational motion
- r i (n) represents the respiratory motion of the i-th feature point
- c i (n) represents the cardiac motion of the i-th feature point
- h i (n) represents the i-th feature
- t i (n) represents the other motion of the i-th feature point
- L(k), C(k), R(k), and H(k) represent the respective harmonics of L(n), c(n), r(n), and h(n), respectively.
- This step includes the following substeps:
- the heart signal spectrum of the j+1th iteration is calculated, and after the discrete Fourier transform is obtained in the frequency range of the frequency point, the j+1th optimal cardiac time domain signal is obtained by the following formula.
- ⁇ is the frequency obtained by the discrete Fourier transform of the signal
- M represents the size of the window, generally set to 3;
- Figure 2 is a contrast image of a patient 1 at an angle (50.8 °, 30.2 °) and automatically extracts the image
- the motion component experiment of the bifurcation feature points (Fig. 3) of each structure of the sequence is 10 frames, and the contrast ratio of the contrast sequence is 12.5 frames/s. Due to the short contrast time, structural feature points A1 to A4 with longer residence time in the contrast sequence are selected for the experiment.
- Figure 11 is a contrast image of a patient 4 at an angle (30.8°, 15.3°) and an exercise component for automatically extracting the bifurcation feature points of each structure of the image sequence.
- the algorithm can extract the two cardiac motion cycles of the patient in the cardiac cycle range of 9 frames and 12 frames respectively, and the sampling rate of the contrast sequence is 12.5 frames/s.
- the size of the peak reflects the distribution of each feature point on the surface of the blood vessel. If the feature point is closer to the lung, the peak value of the feature point is larger.
- Figures (a) and (b) of Figure 8 show the tremor signal after separation, indicating that the patient has significant jitter during the treatment, and the amplitude of the jitter varies with the distribution of the feature points.
- Figure (a) and Figure (b) of Figure 9 show the comparison between the translation curve of each feature point separation and the manual tracking translation curve. It can be seen that the two are basically identical except for the error caused by manual tracking.
- Figures (a) and (b) of Figure 10 and Figure (a) and Figure of Figure 18 (b) Separate the residual curves of each motion for each feature point, the amplitude of which is less than 2 pixels, which can be regarded as the result of algorithmic errors such as segmentation and skeleton extraction, and the motion components are completely separated.
- a patient 4 in Fig. 11 does not have a high frequency component of tremor, but two different frequency components in the frequency range of the heart signal appear, and the appearance of such a signal is The doctor provides a great reference value for diagnosing the patient's condition.
- the method of the present invention considers that under the condition of actual one-arm X-ray angiography, suitable feature points (such as the intersection of ribs, etc.) can not be found in all angiograms, so the selection of vascular structural feature points and Fourier frequency are adopted.
- the domain filtering phase and the multivariate optimization combination approach to extract multiple motion components have wider applicability and flexibility than simple manual tracking to obtain respiratory or translational motion, and are applicable to almost all contrast sequence images.
- this method has higher safety and operability than the method of directly setting the marker point near the heart and then tracking by the relevant imaging means, because the marker added in the tissue in the tissue is generally invasive. It will cause more or less damage to the human body itself, and the process of adding, imaging, eliminating, and extracting the respiratory motion of the marker is complicated, which brings inevitable troubles and errors to the actual operation.
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Multimedia (AREA)
- Apparatus For Radiation Diagnosis (AREA)
- Health & Medical Sciences (AREA)
- General Health & Medical Sciences (AREA)
- Medical Informatics (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Radiology & Medical Imaging (AREA)
- Quality & Reliability (AREA)
- Mathematical Physics (AREA)
- Image Analysis (AREA)
Abstract
公开了一种多参数模型指导的迭代提取血管造影图像运动参数的方法,包括:从造影图序列的医学图像中自动选取I个血管结构特征点,并分别对这些特征点在整个造影图序列中进行自动跟踪,以得到每个特征点的跟踪序列,利用离散傅里叶变换对每个特征点的跟踪序列进行处理,以生成离散傅里叶变换结果S i(k),初始化迭代参数j=0,并获取离散傅里叶变换结果S i(k)中各频点的幅度和频率范围,在各频点的幅度和频率范围中对该频点的跟踪序列进行傅里叶变换,以得到傅里叶变换结果,将获得的傅里叶变换结果进行逆傅里叶变换,并获取该频点的估计最小均方误差。能够解决现有方法中存在的不能自动提取平移运动、以及不能准确分离呼吸和心脏等运动的技术问题。
Description
本发明属于数字信号处理与医学成像的交叉领域,更具体地,涉及一种多参数模型指导的迭代提取血管造影图像运动参数的方法。
在X射线造影系统中,由于呼吸运动的影响,冠状动脉血管在造影面上的图像会发生二维运动。在造影过程中,医生为了能够使造影序列覆盖整个冠脉系统,可能会移动导管床,会导致一个成像序列的不同帧之间整体存在二维平移。因此,冠脉造影图像一方面记录有心脏自身运动在二维平面上的投影,同时也叠加有人体的呼吸运动引起的冠状动脉在造影面上的二维运动和病人的二维平移运动。那么,要从动态二维血管造影图重建三维心血管树,则需将人体的心脏运动,呼吸运动和平移运动等分离开来。
在造影过程中病人平移将直接影响运动分离的结果,而将平移运动分离实际上是进行帧间图像配准,现有的医学图像配准分外部和内部。外部配准方法是在造影前设置一些标记,在造影过程中跟踪它们,但是这种标记法通常是具有侵入性的;内部配准方法分为标记法和分割法等,标记法是在图中选取一些解剖结构点来进行配准。然而并不是每幅图都具有这样的解剖结构点,而且也需要熟悉人体解剖结构。分割法通过分割出来的解剖结构线来配准,既适用于刚体模型也适用于形变模型,但是它只能提取出解剖结构线的运动,不能单独分离出整体的刚性平移运动。然而造影图中表观的血管运动不仅包含了病人平移,还有血管自身的生理运动和呼吸运动。
也有文献中报导了将心脏运动和呼吸运动分离的方法,其中一种方法是预先在体内的一些器官如心脏、膈肌等设置标记点,然后对它们进行跟
踪,把跟踪这些结构特征点所得到的运动曲线近似为呼吸运动。无论是直接手动跟踪造影图中的非心脏结构特征点,还是在造影同时就记录这些点的运动,都是有缺陷的。前者主要在于它的适用性很差;后者的实现需要大量实验控制,对一般的临床应用不合适。
另外一种方法是在双臂x射线造影条件下获得呼吸和心脏运动参数,该方法必须先重构出冠脉的三维运动序列而后分解得到呼吸和心脏运动,处理过程十分复杂。在单臂三维重建中,参与重建的两幅不同角度的造影图像呼吸运动时相不同,所以文献的方法不适用于单臂造影系统。特别地,实际临床应用中,由于造影剂对人体的毒害性,造影时间通常是很短的,我们得到的造影图像的长度有限,从而导致造影图像中的结构特征点的运动是有限长的,这种有限长的运动中包含有周期分量、非周期分量和周期不完整的分量,因此,常规的离散傅里叶变换不能分离。
在申请号为201310750294.5的发明申请“单臂X射线血管造影图像多运动参数分解估计方法”中,虽然也是使用的傅里叶频域分离,但是它不能自动的将各运动参数通过循环迭代的方法自动提取,这样得到的各运动参数不够准确,此方法的鲁棒性不高。申请号为201310752751.4的发明申请“一种X射线血管造影图像中多运动参数的分解估计方法”中采用的是经验模态分解(Empirical Mode Decomposed,简称EMD)的方法估计各参数,同样得到的各运动参数不够准确,鲁棒性不高。
本发明提出了一种自动不断循环优化迭代从X射线造影序列图像中提取平移运动,呼吸运动和心脏运动等多运动参数的方法,通过在冠脉血管上选取一些特征点并在序列中进行跟踪得到它们随时间的运动曲线,然后在多运动参数模型约束下,以离散傅里叶正-反变换和估计的重建信号与原始信号在频域中各频点的局部均方误差最小为准则,通过循环迭代使原始信号与重建信号的差值信号均方最小的优化方法对这些运动分量进行优化处理,得到二维心脏、震颤、呼吸和平移等最优运动曲线。具有很好的临
床适用性。
【发明内容】
针对现有技术的以上缺陷或改进需求,本发明提供了一种多参数模型指导的迭代提取血管造影图像运动参数的方法,其目的在于,解决现有方法中存在的不能自动提取平移运动、以及不能准确分离呼吸和心脏等运动的技术问题。
为实现上述目的,按照本发明的一个方面,提供了一种多参数模型指导的迭代提取血管造影图像运动参数的方法,包括以下步骤:
(1)从造影图序列的医学图像中自动选取I个血管结构特征点,并分别对这些特征点在整个造影图序列中进行自动跟踪,以得到每个特征点的跟踪序列{si(n),i=1,…,I},其中n表示医学图像在造影图序列中的帧号:
(2)利用离散傅里叶变换对步骤(1)获得的每个特征点的跟踪序列{si(n),i=1,…,I}进行处理,以生成离散傅里叶变换结果Si(k):
(3)初始化迭代参数j=0,并获取步骤(2)中获得的离散傅里叶变换结果Si(k)中各频点的幅度和频率范围;
(4)在各频点的幅度和频率范围中对该频点的跟踪序列进行傅里叶变换,以得到傅里叶变换结果;
(5)将步骤(4)获得的傅里叶变换结果进行逆傅里叶变换,并获取该频点的估计最小均方误差;
(6)判断估计最小均方误差是否大于预设的阈值,如果大于,则转入步骤(7),否则过程结束;
(7)利用多参数迭代优化算法对各频点的频谱进行处理,以得到第j+1次迭代后的时域信号;
(8)通过平移模型对剩余信号进行处理,以得到第j+1次迭代后的平移信号;
(9)将第j+1次迭代后的时域信号和第j+1次迭代后的平移信号叠加后得到第j+1次迭代后的估计混合信号,并利用以下公式计算第j+1次迭代后的最小均方误差:
(10)判断第j+1次迭代后的最小均方误差是否大于步骤(6)中的阈值,如果大于则返回步骤(7),否则过程结束。
优选地,步骤(1)是采用以下公式:
si(n)=L(n)+ri(n)+ci(n)+hi(n)+ti(n),i∈[1,I]
其中,L(n)表示平移运动;ri(n)表示第i个特征点的呼吸运动;ci(n)表示第i个特征点的心脏运动;hi(n)表示第i个特征点的震颤运动;ti(n)表示第i个特征点的其他运动。
优选地,步骤(2)是采用以下公式:
Si(k)=L(k)+Ri(k)+Ci(k)+Hi(k)
其中k表示频点,L(k)、C(k)、R(k)和H(k)分别表示L(n)、c(n)、r(n)和h(n)的各次谐波系数。
优选地,步骤(7)包括以下子步骤:
按照本发明的另一方面,提供了一种多参数模型指导的迭代提取血管造影图像运动参数的系统,包括:
第一模块,用于从造影图序列的医学图像中自动选取I个血管结构特征点,并分别对这些特征点在整个造影图序列中进行自动跟踪,以得到每个特征点的跟踪序列{si(n),i=1,…,I},其中n表示医学图像在造影图序列中的帧号:
第二模块,用于利用离散傅里叶变换对第一模块获得的每个特征点的跟踪序列{si(n),i=1,…,I}进行处理,以生成离散傅里叶变换结果Si(k):
第三模块,用于初始化迭代参数j=0,并获取第二模块中获得的离散傅里叶变换结果Si(k)中各频点的幅度和频率范围;
第四模块,用于在各频点的幅度和频率范围中对该频点的跟踪序列进行傅里叶变换,以得到傅里叶变换结果;
第五模块,用于将第四模块获得的傅里叶变换结果进行逆傅里叶变换,并获取该频点的估计最小均方误差;
第六模块,用于判断估计最小均方误差是否大于预设的阈值,如果大于,则转入第七模块,否则过程结束;
第七模块,用于利用多参数迭代优化算法对各频点的频谱进行处理,以得到第j+1次迭代后的时域信号;
第八模块,用于通过平移模型对剩余信号进行处理,以得到第j+1次迭代后的平移信号;
第九模块,用于将第j+1次迭代后的时域信号和第j+1次迭代后的平移信号叠加后得到第j+1次迭代后的估计混合信号,并利用以下公式计算第j+1次迭代后的最小均方误差:
第十模块,用于判断第j+1次迭代后的最小均方误差是否大于第六模块
中的阈值,如果大于则返回第七模块,否则过程结束。
总体而言,通过本发明所构思的以上技术方案与现有技术相比,能够取得下列有益效果:
1、由于采用了步骤(2),本发明自动获取结构特征点得到运动曲线的方法克服了在心脏的表面设置标记进行跟踪获取特征点的运动。
2、采用步骤(7)的多次循环迭代寻优的方法既可以解决文献中不能自动提取平移运动的缺陷;又解决了其他周期运动(如呼吸运动、心脏运动等)不能分离的问题;
3、由于采用了步骤(7)的多次循环迭代寻优的方法解决了专利中所提取的各运动信号不准确且鲁棒性不好的问题。
图1是本发明多参数模型指导的迭代提取血管造影图像运动参数的方法的流程图;
图2为获得的一个角度的原始造影图像;
图3为图2中血管骨架图像选取的特征点;
图4(a)和(b)为图2的特征点A1-A4的X轴和Y轴的原始运动曲线;
图5(a)和(b)为图2的第四个特征点的X轴分量频谱和第四个特征点的Y轴分量频谱;
图6(a)和(b)为图2的各特征点分离出的X轴和Y轴的心脏运动曲线;
图7(a)和(b)为图2的各特征点分离出的X轴和Y轴的呼吸运动曲线;
图8(a)和(b)为图2的各特征点分离出的X轴和Y轴的高频运动曲线;
图9(a)和(b)为图2的各特征点分离出的X轴和Y轴平移曲线与
手动跟踪骨骼点得到的平移曲线的对比。
图10(a)和(b)为图2的各特征点分离出的X轴和Y轴分离的残余曲线。
图11为获得的另一个角度的原始造影图像;
图12为图11中血管骨架图像选取的特征点;
图13(a)和(b)为图11的特征点A1-A5的X轴和Y轴的原始运动曲线;
图14(a)和(b)为图11的第四个特征点的X轴分量频谱和第四个特征点的Y轴分量频谱;
图15(a)和(b)为图11的各特征点分离出的X轴和Y轴的心脏运动曲线;
图16(a)和(b)为图11的各特征点分离出的X轴和Y轴的呼吸运动曲线;
图17(a)和(b)为图11的各特征点分离出的X轴和Y轴分离的残余曲线。
[根据细则26改正13.03.2015]
图18(a)和(b)为图11的各特征点分离各个运动的剩余曲线。
图18(a)和(b)为图11的各特征点分离各个运动的剩余曲线。
为了使本发明的目的、技术方案及优点更加清楚明白,以下结合附图及实施例,对本发明进行进一步详细说明。应当理解,此处所描述的具体实施例仅仅用以解释本发明,并不用于限定本发明。此外,下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。
以下首先就本发明的技术术语进行解释和说明:
最小均方误差:即估计信号与原始信号差值均方累加求和最小。本文多迭代参数模型利用离散傅里叶变换的性质,采取最小均方误差的方法不断在时域和频域迭代使得残差最小的原则求得各分离信号的最优。
为解决现有技术中存在的问题,本发明提出了结构特征点的多运动参
数模型。在冠脉血管造影图上自动选取血管结构特征点,在造影图序列中对其进行自动跟踪,得到它们随时间变化的位移曲线。通过正-反傅里叶变换,重建信号与原始信号在频域中各频点的局部均方误差最小初值估计为基础,结合平移模型,采用循环迭代使原始信号与重建信号的差值信号最小的全局优化得到最终估计的平移运动曲线、呼吸运动曲线以及心脏运动参数。该方法具有广泛的适用性和灵活性,几乎可适用于所有造影序列图(当然需满足有较清晰的血管分布和包含两个或两个以上心动周期,这点不难做到),也包括双臂X射线造影图像。
如图1所示,本发明多参数模型指导的迭代提取血管造影图像运动参数的方法包括以下步骤:
(1)从造影图序列的医学图像中自动选取I个血管结构特征点(其中I为自然数),并分别对这些特征点在整个造影图序列中进行自动跟踪,以得到每个特征点的跟踪序列{si(n),i=1,…,I},其中n表示医学图像在造影图序列中的帧号:
si(n)=L(n)+ri(n)+ci(n)+hi(n)+ti(n),i∈[1,I]
其中,L(n)表示平移运动;ri(n)表示第i个特征点的呼吸运动;ci(n)表示第i个特征点的心脏运动;hi(n)表示第i个特征点的震颤运动;ti(n)表示第i个特征点的其他运动;
(2)利用离散傅里叶变换对步骤(1)获得的每个特征点的跟踪序列{si(n),i=1,…,I}进行处理,以生成离散傅里叶变换结果Si(k):
Si(k)=L(k)+Ri(k)+Ci(k)+Hi(k)
其中k表示频点,L(k)、C(k)、R(k)和H(k)分别表示L(n)、c(n)、r(n)和h(n)的各次谐波系数;
(3)初始化迭代参数j=0,并获取步骤(2)中获得的离散傅里叶变换结果Si(k)中各频点的幅度和频率范围;
本步骤包括以下子步骤:
图2是某病人1在角度(50.8°,30.2°)下的造影图像以及自动提取该图像
序列各结构分叉特征点(图3)的运动分量实验。其中,该病人的心脏运动周期为10帧,造影图序列的采样率为12.5帧/s。由于造影时间较短,本实验中选取在造影序列图中停留时间较长的结构特征点A1~A4用于实验。
图11是某病人4在角度(30.8°,15.3°)下的造影图像以及自动提取该图像序列各结构分叉特征点的运动分量实验。其中,该病人患有心率不齐疾病,利用本文算法可以提取出该病人在心脏周期范围内的两种心脏运动周期分别为9帧和12帧,造影图序列的采样率为12.5帧/s。
从图4中的图(a)和图(b)以及图13中的图(a)和图(b)中可以看到,特征点的原始运动曲线受呼吸运动影响较大,比较凌乱,而且不同特征点的运动幅度不同,这是因为不同特征点位于心脏的不同部位,而心脏不同部位的运动情况是不一样的。图5的图(a)和图(b)以及图14的图(a)和图(b)表示混合信号经过离散傅里叶后的频谱图,从频谱图中可以看出,各频率点呈现的峰值很明显,峰值点的多少说明了该混合信号包含的各频率信号的种类,特别地,如果在0频率点处呈现的峰值很大,则该混合信号一定存在平移信号。图6的图(a)和图(b)、图15的图(a)图(b)以及图16的图(a)和图(b)为分离后的心脏运动曲线,从图中可以看出各特征点运动曲线显示了良好的周期性。图7的图(a)和图(b)以及图17的图(a)和图(b)表示分离后的呼吸运动曲线,从图中可以看出各特征点的呼吸运动曲线不仅具有良好的周期性,而且波峰所在位置基本一致,这是因为造影时间短,各特征点短时间的相位变化非常小。波峰的大小反映了各特征点在血管表面的分布,若特征点离肺部越近,则该特征点的波峰值越大。图8的图(a)和图(b)为分离后的震颤信号,说明了该病人在治疗的过程中出现了比较明显的抖动,抖动的幅度随着各特征点的分布而不一样。图9的图(a)和图(b)分别为各特征点分离的平移曲线与手动跟踪的平移曲线的对比,可以看出两者除了因为手动跟踪引起的误差外,基本一致。图10的图(a)和图(b)以及图18的图(a)和图
(b)分别为各特征点分离各个运动的剩余曲线,其幅度都少于2个像素,可以看作是因为分割,骨架提取等算法误差引起的,运动分量基本完全分离。
不同于图2中的某病人1的是,图11中的某病人4没有出现震颤的高频成分,却出现了在心脏信号频率范围内的两种不同频率成分,这类信号的出现,为医生诊断病人病情提供了极大的参考价值。
总而言之,本发明方法考虑到在实际单臂X射线造影的条件下,并非在所有造影图中都能找到合适的特征点(像肋骨的交点等),所以采用血管结构特征点的选取和傅立叶频域滤波相以及多变量寻优结合的途径来提取多个运动分量,比单纯的用手动跟踪来得到呼吸运动或者平移运动具有更广泛的适用性和灵活性,几乎可适用于所有造影序列图,同时此方法相较于直接在心脏附近组织设置标识点再通过相关成像手段进行跟踪的方法拥有更高安全性和可操作性,这是因为在体内组织添加的标记物一般是可侵入性的,会对人体自身产生或多或少的损害,并且其标记物的添加、成像、排除、呼吸运动提取整个过程都是繁杂的,为实际操作中带来不可避免的麻烦与误差。
本领域的技术人员容易理解,以上所述仅为本发明的较佳实施例而已,并不用以限制本发明,凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等,均应包含在本发明的保护范围之内。
Claims (7)
- 一种多参数模型指导的迭代提取血管造影图像运动参数的方法,其特征在于,包括以下步骤:(1)从造影图序列的医学图像中自动选取I个血管结构特征点,并分别对这些特征点在整个造影图序列中进行自动跟踪,以得到每个特征点的跟踪序列{si(n),i=1,…,I},其中n表示医学图像在造影图序列中的帧号:(2)利用离散傅里叶变换对步骤(1)获得的每个特征点的跟踪序列{si(n),i=1,…,I}进行处理,以生成离散傅里叶变换结果Si(k):(3)初始化迭代参数j=0,并获取步骤(2)中获得的离散傅里叶变换结果Si(k)中各频点的幅度和频率范围;(4)在各频点的幅度和频率范围中对该频点的跟踪序列进行傅里叶变换,以得到傅里叶变换结果;(5)将步骤(4)获得的傅里叶变换结果进行逆傅里叶变换,并获取该频点的估计最小均方误差;(6)判断估计最小均方误差是否大于预设的阈值,如果大于,则转入步骤(7),否则过程结束;(7)利用多参数迭代优化算法对各频点的频谱进行处理,以得到第j+1次迭代后的时域信号;(8)通过平移模型对剩余信号进行处理,以得到第j+1次迭代后的平移信号;(9)将第j+1次迭代后的时域信号和第j+1次迭代后的平移信号叠加后得到第j+1次迭代后的估计混合信号,并利用以下公式计算第j+1次迭代后的最小均方误差:(10)判断第j+1次迭代后的最小均方误差是否大于步骤(6)中的阈值,如果大于则返回步骤(7),否则过程结束。
- 根据权利要求1所述的方法,其特征在于,步骤(1)是采用以下公式:si(n)=L(n)+ri(n)+ci(n)+hi(n)+ti(n),i∈[1,I]其中,L(n)表示平移运动;ri(n)表示第i个特征点的呼吸运动;ci(n)表示第i个特征点的心脏运动;hi(n)表示第i个特征点的震颤运动;ti(n)表示第i个特征点的其他运动。
- 根据权利要求2所述的方法,其特征在于,步骤(2)是采用以下公式:Si(k)=L(k)+Ri(k)+Ci(k)+Hi(k)其中k表示频点,L(k)、C(k)、R(k)和H(k)分别表示L(n)、c(n)、r(n)和h(n)的各次谐波系数。
- 根据权利要求4所述的方法,其特征在于,步骤(7)包括以下子步骤:
- 一种多参数模型指导的迭代提取血管造影图像运动参数的系统,其特征在于,包括:第一模块,用于从造影图序列的医学图像中自动选取I个血管结构特征点,并分别对这些特征点在整个造影图序列中进行自动跟踪,以得到每个特征点的跟踪序列{si(n),i=1,…,I},其中n表示医学图像在造影图序列中的帧号:第二模块,用于利用离散傅里叶变换对第一模块获得的每个特征点的跟踪序列{si(n),i=1,…,I}进行处理,以生成离散傅里叶变换结果Si(k):第三模块,用于初始化迭代参数j=0,并获取第二模块中获得的离散傅里叶变换结果Si(k)中各频点的幅度和频率范围;第四模块,用于在各频点的幅度和频率范围中对该频点的跟踪序列进行傅里叶变换,以得到傅里叶变换结果;第五模块,用于将第四模块获得的傅里叶变换结果进行逆傅里叶变换,并获取该频点的估计最小均方误差;第六模块,用于判断估计最小均方误差是否大于预设的阈值,如果大于,则转入第七模块,否则过程结束;第七模块,用于利用多参数迭代优化算法对各频点的频谱进行处理,以得到第j+1次迭代后的时域信号;第八模块,用于通过平移模型对剩余信号进行处理,以得到第j+1次迭代后的平移信号;第九模块,用于将第j+1次迭代后的时域信号和第j+1次迭代后的平移信号叠加后得到第j+1次迭代后的估计混合信号,并利用以下公式计算第j+1次迭代后的最小均方误差:第十模块,用于判断第j+1次迭代后的最小均方误差是否大于第六模块 中的阈值,如果大于则返回第七模块,否则过程结束。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| US14/960,461 US20160189394A1 (en) | 2014-12-30 | 2015-12-07 | Method for iteratively extracting motion parameters from angiography images |
Applications Claiming Priority (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN2014108441394 | 2014-12-30 | ||
| CN201410844139.4A CN104517301B (zh) | 2014-12-30 | 2014-12-30 | 多参数模型指导的迭代提取血管造影图像运动参数的方法 |
Related Child Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| US14/960,461 Continuation-In-Part US20160189394A1 (en) | 2014-12-30 | 2015-12-07 | Method for iteratively extracting motion parameters from angiography images |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| WO2016106959A1 true WO2016106959A1 (zh) | 2016-07-07 |
Family
ID=52792546
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| PCT/CN2015/072681 Ceased WO2016106959A1 (zh) | 2014-12-30 | 2015-02-10 | 多参数模型指导的迭代提取血管造影图像运动参数的方法 |
Country Status (2)
| Country | Link |
|---|---|
| CN (1) | CN104517301B (zh) |
| WO (1) | WO2016106959A1 (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114170192A (zh) * | 2021-12-10 | 2022-03-11 | 武汉工程大学 | 星载光学相机焦面震颤精细检测方法 |
Families Citing this family (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107569251B (zh) * | 2017-08-29 | 2020-12-22 | 上海联影医疗科技股份有限公司 | 医学成像方法和系统及非暂态计算机可读存储介质 |
| CN108852385B (zh) * | 2018-03-13 | 2022-03-04 | 中国科学院上海应用物理研究所 | 一种x射线造影方法及基于x射线造影的动态阅片方法 |
| US12020136B2 (en) * | 2018-09-28 | 2024-06-25 | Samsung Electronics Co., Ltd. | Operating method and training method of neural network and neural network thereof |
| CN112716509B (zh) * | 2020-12-24 | 2023-05-02 | 上海联影医疗科技股份有限公司 | 一种医学设备的运动控制方法及系统 |
| CN112998853B (zh) * | 2021-02-25 | 2023-05-23 | 四川大学华西医院 | 腹部血管动态造影2d建模方法、3d建模方法及检测系统 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101773395A (zh) * | 2009-12-31 | 2010-07-14 | 华中科技大学 | 一种从单臂x射线造影图中提取呼吸运动参数的方法 |
| US20110208678A1 (en) * | 2010-02-19 | 2011-08-25 | Oracle International Corporation | Mechanical shock feature extraction for overstress event registration |
| CN103747742A (zh) * | 2011-04-14 | 2014-04-23 | 明尼苏达大学评议会 | 使用超声成像的脉管表征 |
| CN103810721A (zh) * | 2013-12-30 | 2014-05-21 | 华中科技大学 | 单臂x射线血管造影图像多运动参数分解估计方法 |
| CN103886615A (zh) * | 2013-12-31 | 2014-06-25 | 华中科技大学 | 一种x射线血管造影图像中多运动参数的分离估计方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7209777B2 (en) * | 2000-11-30 | 2007-04-24 | General Electric Company | Method and apparatus for automated tracking of non-linear vessel movement using MR imaging |
-
2014
- 2014-12-30 CN CN201410844139.4A patent/CN104517301B/zh active Active
-
2015
- 2015-02-10 WO PCT/CN2015/072681 patent/WO2016106959A1/zh not_active Ceased
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101773395A (zh) * | 2009-12-31 | 2010-07-14 | 华中科技大学 | 一种从单臂x射线造影图中提取呼吸运动参数的方法 |
| US20110208678A1 (en) * | 2010-02-19 | 2011-08-25 | Oracle International Corporation | Mechanical shock feature extraction for overstress event registration |
| CN103747742A (zh) * | 2011-04-14 | 2014-04-23 | 明尼苏达大学评议会 | 使用超声成像的脉管表征 |
| CN103810721A (zh) * | 2013-12-30 | 2014-05-21 | 华中科技大学 | 单臂x射线血管造影图像多运动参数分解估计方法 |
| CN103886615A (zh) * | 2013-12-31 | 2014-06-25 | 华中科技大学 | 一种x射线血管造影图像中多运动参数的分离估计方法 |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114170192A (zh) * | 2021-12-10 | 2022-03-11 | 武汉工程大学 | 星载光学相机焦面震颤精细检测方法 |
| CN114170192B (zh) * | 2021-12-10 | 2024-11-15 | 武汉工程大学 | 星载光学相机焦面震颤精细检测方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN104517301A (zh) | 2015-04-15 |
| CN104517301B (zh) | 2017-07-07 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| KR102114415B1 (ko) | 의료 영상 정합 방법 및 장치 | |
| KR101428005B1 (ko) | 소수의 저선량 ct 영상을 이용하여 pet 영상을 움직임 보상 및 감쇠 보정하는 방법 | |
| US8515146B2 (en) | Deformable motion correction for stent visibility enhancement | |
| CN103886615B (zh) | 一种x射线血管造影图像中多运动参数的分离估计方法 | |
| Groher et al. | Deformable 2D-3D registration of vascular structures in a one view scenario | |
| KR102233427B1 (ko) | 의료 영상 정합 방법 및 장치 | |
| JP2019500146A (ja) | 体部の3次元モデル | |
| CN101773395B (zh) | 一种从单臂x射线造影图中提取呼吸运动参数的方法 | |
| WO2016106959A1 (zh) | 多参数模型指导的迭代提取血管造影图像运动参数的方法 | |
| US20100189337A1 (en) | Method for acquiring 3-dimensional images of coronary vessels, particularly of coronary veins | |
| US20130101187A1 (en) | Coronary artery motion modeling | |
| US20040087850A1 (en) | Method and apparatus for medical intervention procedure planning | |
| JP6271097B2 (ja) | デジタルサブトラクション血管造影 | |
| CN104519793A (zh) | 用于内腔的线性映射的系统 | |
| JP2014140742A (ja) | 動いている器官の目標エリア内のオブジェクトを追跡するための方法および装置 | |
| KR101900679B1 (ko) | 혈관 특징 정보를 기반으로 하는 삼차원 심혈관 정합 방법, 이를 수행하기 위한 기록 매체 및 장치 | |
| CN108294768B (zh) | 序列图像多参数配准的x射线心血管造影减影方法及系统 | |
| Banerjee et al. | Point-cloud method for automated 3D coronary tree reconstruction from multiple non-simultaneous angiographic projections | |
| CN120599292B (zh) | 一种基于目标检测的血管通路病变识别方法及系统 | |
| Zhang et al. | 3D ultrasound centerline tracking of abdominal vessels for endovascular navigation | |
| Klugmann et al. | Deformable respiratory motion correction for hepatic rotational angiography | |
| KR101785293B1 (ko) | 3d cta 영상 정보를 활용한 2d xa 영상에서의 자동 혈관 구조 해석 방법, 이를 수행하기 위한 기록 매체 및 장치 | |
| CN113538419A (zh) | 一种图像处理方法和系统 | |
| Zhang et al. | A novel structural features-based approach to automatically extract multiple motion parameters from single-arm X-ray angiography | |
| CN114040707B (zh) | 诊断支持程序 |
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: 15874652 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: 15874652 Country of ref document: EP Kind code of ref document: A1 |
|
| 122 | Ep: pct application non-entry in european phase |
Ref document number: 15874652 Country of ref document: EP Kind code of ref document: A1 |