WO2015101059A1 - 一种x射线血管造影图像中多运动参数的分离估计方法 - Google Patents

一种x射线血管造影图像中多运动参数的分离估计方法 Download PDF

Info

Publication number
WO2015101059A1
WO2015101059A1 PCT/CN2014/085724 CN2014085724W WO2015101059A1 WO 2015101059 A1 WO2015101059 A1 WO 2015101059A1 CN 2014085724 W CN2014085724 W CN 2014085724W WO 2015101059 A1 WO2015101059 A1 WO 2015101059A1
Authority
WO
WIPO (PCT)
Prior art keywords
motion
curve
respiratory
sequence
cardiac
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/CN2014/085724
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.)
Huazhong University of Science and Technology
Original Assignee
Huazhong University of Science and Technology
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 Huazhong University of Science and Technology filed Critical Huazhong University of Science and Technology
Priority to US14/698,894 priority Critical patent/US9286676B2/en
Publication of WO2015101059A1 publication Critical patent/WO2015101059A1/zh
Anticipated expiration legal-status Critical
Ceased legal-status Critical Current

Links

Images

Classifications

    • 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/246—Analysis of motion using feature-based methods, e.g. the tracking of corners or segments
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/50—Image enhancement or restoration using two or more images, e.g. averaging or subtraction
    • 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/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/20—Special algorithmic details
    • G06T2207/20212—Image combination
    • G06T2207/20224—Image subtraction
    • 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/30048—Heart; Cardiac
    • 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
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2211/00—Image generation
    • G06T2211/40—Computed tomography
    • G06T2211/404—Angiography

Definitions

  • the invention belongs to the technical field of contrast image processing, and more particularly to a method for separating and estimating multi-motion parameters in an X-ray angiography image.
  • the thoracic rhythm is enlarged and reduced to complete inhalation and exhalation. This is the breathing movement.
  • Respiratory motion causes an overall translation of the heart in three dimensions.
  • two-dimensional motion of the coronary vessels occurs on the contrast surface due to 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.
  • 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 translational 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, in order to get a two-dimensional angiogram closer to the real situation, it is necessary to separate the heart, respiratory and translational movements of the human body.
  • the general method of separating translational motion is to perform inter-frame image registration.
  • Many medically documented medical image registrations are usually external and internal.
  • the external registration method is to track the contrast during the angiography by setting some obvious marks before the angiography, but this method is usually invasive; the internal registration method is divided into the marking method and the segmentation method after the angiography is completed. Wait.
  • the notation is to select some anatomical points in the figure for registration. However, not every figure has such anatomical points, and it is also necessary to understand the human anatomy.
  • the segmentation method can be registered by the segmented anatomical line. It can be applied to the rigid body model or the deformation model, but it can only extract the motion of the anatomical line, but can not separate the whole rigid translational motion.
  • the vascular motion extracted from the angiogram includes not only the patient's translational motion, but also the physiological and respiratory movements of the blood vessel itself.
  • some angiograms contain bones such as the spine that do not have physiological motion, they can be registered as internal markers by internal markers, but not all angiograms are stored. Registration can be performed on such features.
  • the translational motion is difficult to extract.
  • the extraction of such translational motion is rarely mentioned in the field of medical imaging.
  • some organs in the body such as the heart and diaphragm, are synchronized with the movement of the lungs to perform three-dimensional motion. Therefore, tracking the structural feature points that do not move with the heart can obtain the motion curve, which can be approximated as respiratory motion.
  • more used in X-ray angiography is to directly track non-cardiac structural feature points in the angiogram or record the motion of these points while angiography. Obviously, both options are flawed.
  • the former is mainly because its applicability is very poor, because we can't guarantee that there are markers in this frame that meet this condition, and it is necessary to understand the human anatomy. Therefore, when the above feature points are not present in the angiogram, the respiratory motion is difficult to extract.
  • the latter implementation requires a lot of experimental control, which is not suitable for general clinical applications.
  • the present invention provides an X-ray vessel
  • the method for separating and estimating multi-motion parameters in angiographic images is based on frequency domain filtering and multivariate optimization, which can accurately separate each motion parameter, which is simple, fast and harmless.
  • the invention provides a method for separating and estimating multi-motion parameters in an X-ray angiography image, comprising the following steps:
  • the above technical solution conceived by the present invention is compared with the prior art, since it is considered that the X-ray angiography image reflects various motion signals such as a translational motion signal, a cardiac motion signal, a respiratory motion signal, and a high frequency.
  • the combination of motion signals and other factors proposes the separation and estimation of multiple motion parameters, which can make the separated movements more accurate and provide accurate and effective help for doctors' diagnosis.
  • FIG. 1 is a block diagram showing an implementation flow of a method for separating and estimating multi-motion parameters in an X-ray angiography image according to an embodiment of the present invention
  • Fig. 2(a) is a mark of a contrast image and a bone point in the Fourier transform extraction, and the corresponding projection angle is (-26.5°, -20.9°);
  • Fig. 2(b) is a marker of the angiogram and the skeletal points in the Fourier transform extraction, and the corresponding projection angle is (42.3°, 26.8°);
  • Fig. 3(a) is a skeleton diagram and a feature point mark of Fig. 2(a).
  • Fig. 3(b) is a skeleton diagram and feature point marks of Fig. 2(b).
  • Fig. 4(a) is the original curve obtained by tracking all the marked points in the contrast sequence in Fig. 3(a), wherein the solid line is the change of the lateral coordinate (X-axis coordinate) of each feature point, and the broken line is the characteristic. The change in the longitudinal coordinate of the point (Y-axis coordinate).
  • Fig. 4(b) is the original curve obtained by tracking all the marker points in the contrast sequence in Fig. 3(b), wherein the solid line is the change of the lateral coordinate (X-axis coordinate) of each feature point, and the broken line is The change in the longitudinal coordinate (Y-axis coordinate) of each feature point.
  • Fig. 4(e) is a cardiac motion curve obtained by decomposing the original curve in Fig. 3(a), wherein the solid line is the change of the lateral coordinate (X-axis coordinate) of each feature point, and the broken line is the longitudinal coordinate of each feature point ( Change in the Y-axis coordinate).
  • Fig. 4(f) is a curve of the heart motion obtained by decomposing the original curve in Fig. 3(b), wherein the solid line is the change of the lateral coordinate (X-axis coordinate) of each feature point, and the broken line is the longitudinal coordinate of each feature point ( Change in the Y-axis coordinate).
  • Fig. 4(g) is a respiratory motion signal curve obtained by decomposing the original curve in Fig. 3(a), wherein the solid line is the change of the lateral coordinate (X-axis coordinate) of each feature point, and the broken line is the longitudinal coordinate of each feature point. (Y-axis coordinate) change.
  • Fig. 4(h) is a respiratory motion signal curve obtained by decomposing the original curve in Fig. 3(b), wherein the solid line is the change of the lateral coordinate (X-axis coordinate) of each feature point, and the broken line is the longitudinal coordinate of each feature point. (Y-axis coordinate) change.
  • Fig. 4(i) is a high-frequency component motion signal curve obtained by decomposing the original curve in Fig. 3(a), wherein the solid line is the change of the lateral coordinate (X-axis coordinate) of each feature point, and the broken line is the characteristic point of each feature point. The change in the longitudinal coordinate (Y-axis coordinate).
  • Fig. 4(j) is a high-frequency component motion signal curve obtained by decomposing the original curve in Fig. 3(b), wherein the solid line is the change of the lateral coordinate (X-axis coordinate) of each feature point, and the broken line is the characteristic point of each feature point. The change in the longitudinal coordinate (Y-axis coordinate).
  • Fig. 4(k) is a comparison of the X-axis translational motion signal curve obtained by decomposing the original curve in Fig. 3(a) with the X-axis curve obtained by tracking the skeleton points in Fig. 2(a), wherein the solid line is decomposed The X-axis shifts the motion signal curve, and the dashed line is the X-axis curve obtained by tracking the bone points.
  • Fig. 4(l) is a comparison of the X-axis translational motion signal curve obtained by decomposing the original curve in Fig. 3(b) with the X-axis curve obtained by tracking the skeleton point in Fig. 2(b), wherein the solid line is decomposed The X-axis shifts the motion signal curve, and the dashed line is the X-axis curve obtained by tracking the bone points.
  • Figure 4 (m) is the Y-axis translational motion signal curve obtained by decomposing the original curve in Figure 3 (a) The line is compared with the Y-axis curve obtained by tracking the bone points in Fig. 2(a), where the solid line is the Y-axis translational motion signal curve obtained by the decomposition, and the broken line is the Y-axis curve obtained by tracking the skeleton point.
  • Figure 4(n) is a comparison of the Y-axis translational motion signal curve obtained by decomposing the original curve in Figure 3(b) with the Y-axis curve obtained by tracking the skeleton points in Figure 2(b), where the solid line is decomposed
  • the Y-axis translational motion signal curve, the dashed line is the Y-axis curve obtained by tracking the bone points.
  • Fig. 4(o) is a residual noise curve obtained by decomposing the original curve in Fig. 3(a), where the solid line is the change of the lateral coordinate (X-axis coordinate) of each feature point, and the broken line is the longitudinal direction of each feature point. The change in coordinates (Y-axis coordinates).
  • Fig. 4(p) is a residual noise curve obtained by decomposing the original curve in Fig. 3(b) to obtain a residual noise curve, wherein the solid line is the change of the lateral coordinate (X-axis coordinate) of each feature point, and the broken line is the characteristic point of each feature point. Change in longitudinal coordinate (Y-axis coordinate);
  • the invention particularly relates to a data separation method for selecting a feature point of a blood vessel structure and combining Fourier series expansion.
  • the method can solve the separation problem of the human body's translational motion signals, cardiac motion signals and respiratory motion signals.
  • the method is based on frequency domain filtering and multivariate optimization, and can accurately separate each motion parameter, which is simple, fast, and harmless, and can be widely applied regardless of single arm or double arm imaging.
  • X-ray angiography images not only record the projection of the heart's own motion on the coronary arteries, It also includes 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.
  • the steps include:
  • I is the number of feature points
  • N is the total number of frames of the image in the sequence of X-ray angiography images.
  • the vascular structural feature points need to be able to truly reflect the motion information of the whole blood vessel. Therefore, according to the knowledge of physiology and anatomy, the selected feature points are the bifurcation points of blood vessels.
  • the vascular structural feature points are obtained by first pre-processing the X-ray angiography image, obtaining the vascular skeleton by segmentation and refinement, and then selecting the vascular structural feature points through the vascular skeleton.
  • the tracking method adopted by the present invention is to obtain the coordinates of each feature point in each frame by establishing a correspondence relationship for each bifurcation point of different frames of the X-ray angiography image skeleton sequence.
  • Step 3 Determine the heart motion signal period N c according to the contrast sequence.
  • Step6 The comprehensive motion estimation curve is Take The closest to s i (n) is the optimal criterion, and the best curve translational motion signal curve, cardiac motion curve, respiratory motion signal curve and high frequency component motion signal are obtained.
  • the method of the invention takes into account that under the condition of actual X-ray angiography, suitable feature points (such as intersections of ribs, etc.) can not be found in all angiograms, so the selection of vascular structural feature points and Fourier frequency domain filtering are used.
  • the multivariate optimization approach combines multiple motion components to extract a wider range of applicability and flexibility than manual tracking to obtain respiratory or translational motion, and is applicable to almost all contrast sequence maps.
  • the method of setting the identification points in the near tissue and then tracking by the relevant imaging means has higher safety and operability, because the markers added in the tissue in the body are generally invasive and will generate more or less to the human body itself. The damage, and the addition, imaging, elimination, and respiratory motion extraction of the markers are complicated, which brings inevitable troubles and errors to the actual operation.
  • the invention provides a method for extracting multi-motion parameters such as translational motion signals, respiratory motion signals and cardiac motion signals directly from X-ray contrast sequence images, by selecting some feature points on the coronary vessels and tracking them in the sequence. Their motion curves over time are then separated by Fourier series expansion, frequency domain filtering and multivariate optimization to obtain separate translation signals, cardiac signals and respiratory motion signals. Has good clinical applicability.
  • the clinically measured normal cardiac motion is 60-100 beats/min with a period of 0.6-1.0 s, while the periodic motion with a frequency higher than 140 beats/min is diagnosed as pathological abnormal motion with a period of less than 0.42 s.
  • the coronary angiography can be used to analyze the exercise frequency higher than the normal cardiac motion, and the amplitude of the high-frequency signal exceeds 2-3 pixels, it can be reasonably considered that the high-frequency signal is a useful signal (or diagnostic lesion signal). This article is called high frequency component). This is because, in the process of processing the coronary angiography image, the inaccuracy of segmentation or skeleton extraction may cause a certain error, but the error generally does not exceed 3 pixels.
  • the frequency of the respiratory motion signal is much smaller than the two motions, and the period of the respiratory motion signal is much longer, usually 3-6 s, and may be longer when quiet.
  • each motion component including translational signal, cardiac motion signal, respiratory motion signal and high-frequency component motion signal
  • these components are separated by Fourier series expansion combined with frequency domain filtering.
  • continuous time period signals can be expressed by Fourier series.
  • the periodic sequence can also be represented by discrete Fourier series. If the period of the periodic sequence x p (n) is N p , then
  • the motion curve of the x-axis coordinate of the i-th feature point P i (x, y) on the coronary blood vessel in the contrast image sequence is x i (n) (n is the number of frames of the contrast frame, representing the time of discretization, sequence The total length is N)
  • the motion curve of the Y-axis coordinate is y i (n)
  • s i (n) (x i (n), y i (n))
  • s i (n) can be expressed as:
  • I denotes the number of feature points
  • r i (n) (x ri (n)
  • y ri (n)) represents the motion caused by respiratory motion
  • c i (n) (x ci (n), y ci (n)) represents the motion of the points on the blood vessel caused by the movement of the heart
  • s i (n) is used to represent the coordinate curve of the extracted vascular feature points along the X-axis and the Y-axis
  • t i (n) respectively represent coordinate curves of the translational motion signal, the respiratory motion signal, the cardiac motion signal, the high-frequency component motion signal, and the noise component along the X-axis and the Y-axis.
  • Is the harmonic component of the heart motion signal Is the harmonic component of the respiratory motion signal
  • Is the harmonic component of the respiratory motion signal It is the harmonic component of the high-frequency component motion signal
  • C(k), R(k), and H(k) represent the harmonic coefficients of c(n), r(n), and h(n), respectively.
  • their harmonic components cannot coincide (assuming that the translational motion d(n) has been eliminated), ie
  • N c , N r and N h are mutually homogeneous, the second condition is almost impossible to satisfy. Because under normal circumstances, the contrast agent stays in the human body for not too long, a sequence of contrast images will only last for about 6s, no more than 10s, assuming a frame rate of 80ms / frame (actually possibly smaller), then N c ⁇ 10, N r >30, N c ⁇ N r ⁇ N h >N c ⁇ N r >300>10/0.08, it is difficult to obtain a complete cycle sequence. Moreover, the patient's flat-moving signal curve, the respiratory motion signal cycle determination requirements are higher.
  • cardiac dynamic signals Although it is impossible to directly separate cardiac dynamic signals, respiratory signals, high-frequency component signals, and moving signal curves, they can be expanded by Fourier series, frequency domain filtering and multivariate optimization according to their respective motion characteristics. The method is used for approximate separation.
  • the separation algorithm is based on the following assumptions:
  • the translational motion signal estimation is assumed to be a polyline motion, as follows:
  • I is the number of feature points
  • N is the total number of frames of the image in the sequence of X-ray angiography images.
  • the vascular structural feature points need to be able to truly reflect the movement information of the whole blood vessel. Therefore, according to the knowledge of physiology and anatomy, the selected feature point is the bifurcation point of the blood vessel. As shown in Fig.
  • the vascular structural feature points are obtained by first pre-processing the X-ray angiography image, obtaining the vascular skeleton by segmentation and refinement, and then selecting the vascular structural feature points through the vascular skeleton.
  • the tracking method adopted by the present invention is to obtain the coordinates of each feature point in each frame by establishing a correspondence relationship for each bifurcation point of different frames of the X-ray angiography image skeleton sequence.
  • s i (n) can be expressed as:
  • Step2 estimation of the translational motion signal curve
  • the sequence of the X-axis movement change frame is determined by the contrast image Where T x is the number of times the X-axis movement direction is transformed, and the sequence of the Y-axis movement change frame Where T y is the number of times the Y-axis movement direction is changed.
  • Slope angle The range is (- ⁇ /2, ⁇ /2) to obtain the X-axis simulation translation curve And Y-axis simulation translation curve As in formula (7). Make Indicates the patient's translational motion signal.
  • the cardiac motion signal period N c is determined according to the sequence of the contrast map.
  • N c ⁇ 0.8/f, where N c represents the number of contrast frames in one cardiac cycle, and 0.8 represents the time of one cardiac cycle in seconds.
  • f denotes the frame rate, that is, the number of times the contrast image is taken every second.
  • Pass on Select a sequence of length n sc n sc is an integer multiple of N c .
  • the frequency domain response of cardiac motion is:
  • Step5 Estimation of high-frequency component motion signal curve
  • the optimal high frequency component motion signal period is obtained by the following equation (11) as the optimal criterion. And the best high frequency component motion signal curve
  • Step6 The motion estimation curve is The best translational motion signal curve is obtained by the optimal criterion of (12) And the corresponding heart motion curve Respiratory motor signal curve And high frequency component motion signal curve
  • the cardiac motion cycle was determined to be 11 frames by observing the sequence of the contrast image, and FIG. 4 is the result of motion separation of the motion of the marker points in FIG.
  • FIG. 4(j) show that the separated high frequency component is about 0.54 times of the heart motion, the period is about 6frames, and the amplitude is more than 3 pixels, if it is caused by segmentation or skeleton extraction.
  • the error should be in the range of 3 pixels, as shown in Figure 4(k), Figure 4(l), Figure 4(m) and Figure 4(n), so it can be assumed that it is some kind of regular abnormal motion of the heart (such as Heart tremors, etc.).
  • Figure 4(k), Figure 4(l), Figure 4(m) and Figure 4(n) show the contrast between the separated translation curve and the manually tracked translation curve. It can be seen that the two are in addition to the error caused by manual tracking. , basically the same.
  • Fig. 4(o) and Fig. 4(p) are the remaining curves separating the motions, all of which are less than 3 pixels in amplitude, which can be regarded as caused by algorithm errors such as segmentation and skeleton extraction, and the motion components are substantially completely separated.

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Multimedia (AREA)
  • Mathematical Physics (AREA)
  • Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Radiology & Medical Imaging (AREA)
  • Quality & Reliability (AREA)
  • General Health & Medical Sciences (AREA)
  • Apparatus For Radiation Diagnosis (AREA)
  • Image Analysis (AREA)

Abstract

一种X射线血管造影图像中多运动参数的分离估计方法，所述方法包括：首先根据所述的造影图像序列确定心脏运动信号周期以及平移运动变化帧的序列，然后通过对所述的造影图像序列中的血管结构特征点跟踪得到点的运动序列，再通过多变量参数寻优以及傅里叶频域滤波对所述的运动序列进行处理，根据所述的平移运动变化帧的序列，心脏运动信号的周期，呼吸运动信号周期的范围和高频运动信号周期的范围，可以分离出最佳的平移运动曲线、心脏运动曲线、呼吸运动曲线和高频运动曲线。所述方法由于考虑到X射线血管造影图像是反映了多种运动信号的综合，提出的多运动参数的分离估计，能够使分离出来的各个运动更加精确，为医生的诊断提供准确有效的帮助。

Description

一种X射线血管造影图像中多运动参数的分离估计方法 [技术领域]
本发明属于造影图像处理技术领域，更具体地，涉及一种X射线血管造影图像中多运动参数的分离估计方法。
[背景技术]
胸廓有节律的扩大和缩小，从而完成吸气与呼气，这就是呼吸运动。呼吸运动会造成心脏在三维空间中整体的平移，在X射线造影系统中，由于呼吸运动的影响冠状动脉血管在造影面上会发生二维运动。而且在造影过程中，医生为了能够使造影序列覆盖整个冠脉系统，可能会移动导管床，导致一个成像序列的不同帧之间整体存在二维平移。
因此，冠脉造影图像一方面记录有心脏自身运动在二维平面上的投影，同时也叠加有人体的呼吸运动引起的冠状动脉在造影面上的二维平移运动和病人的二维平移运动。那么，要得到更接近真实情况下的二维血管造影图，则需将人体的心脏运动，呼吸运动和平移运动等分离开来。
分离平移运动一般方法是进行帧间图像配准。许多文献报道的医学图像配准通常是分外部和内部。外部配准方法是通过在造影前设置一些明显的标记，在造影过程中跟踪，但是这种标记法通常是具有侵入性的；内部配准方法是在造影完成后又分为标记法和分割法等。标记法是在图中选取一些解剖结构点来进行配准，然而并不是每幅图都具有这样的解剖结构点，而且也需要对人体解剖结构比较了解。而分割法可以通过分割出来的解剖结构线来配准，可以适用于刚体模型也可适用于形变模型，但是它只能提取出解剖结构线的运动，却不能单独分离出整体的刚性平移运动。然而造影图中提取出来的血管运动不仅包含了病人平移运动，还有血管自身的生理运动和呼吸运动。虽然有些造影图中包含了不会有生理运动的脊椎等骨骼可以作为内部标记点来通过内部标记法配准，但不是所有的造影图都存 在这样的特征可以进行配准。当造影图不存在这样的特征时，平移运动是很难被提取出来的。而且这样的平移运动的提取在医学图像领域上少有提及。
文献报导在提取人体呼吸运动的时候都要通过预先设置标记点，然后对它们进行序列跟踪。根据呼吸运动的特点，人在进行呼吸的时候体内的一些器官如心脏，膈肌等会随着肺的运动一起同步进行三维空间运动。所以跟踪不会随着心脏一起运动的结构特征点得到其运动曲线，可近似认作呼吸运动。其中在X射线造影中用的比较多的是，直接手动跟踪造影图中的非心脏结构特征点或者在造影同时就记录这些点的运动。很显然，这两种方案都是有缺陷的。前者主要在于它的适用性很差，因为我们并不能保证每一帧造影图中都存在符合这种条件的标记点，而且找这种点也需要对人体解剖结构比较了解。所以，在造影图中不存在以上特征点时，呼吸运动是很难被提取出来的。而后者的实现则需要大量实验控制，对一般的临床应用不合适。此外，还有一种方法是在双臂X射线造影条件下实现的，其分离心脏运动与呼吸运动的思想是，取同一时刻不同投影角度的两幅造影图，对其中相对应的冠脉血管进行三维重建，获得该时刻的血管三维空间分布，那么，对一个呼吸周期中的所有造影图对进行匹配和重建后，得到的则是一组三维结构序列，而它们间的空间位移矢量便是呼吸运动。相对来说，通过该方法能比较可靠的呼吸运动结果。但是由于双臂x射线造影条件的约束，不能广泛的应用在医学实践中。
在专利一种从X射线造影图中提取呼吸运动参数的方法200910273528.5中，虽然也是使用的傅里叶频域分离，但是它只能通过提取心脏运动后再将剩余曲线视为呼吸运动，且不能分离出平移运动和高频成分，这样得到的呼吸运动不够准确，此方法的鲁棒性不高。
[发明内容]
针对现有技术的以上缺陷或改进需求，本发明提供了一种X射线血管 造影图像中多运动参数的分离估计方法，其目的在于基于频域滤波和多变量寻优，能够精确分离各个运动参数，简单快速、无伤害。
本发明提供了一种X射线血管造影图像中多运动参数的分离估计方法，包括下述步骤：
(1)对血管结构特征点在X射线血管造影图像序列中进行跟踪，得到点的跟踪曲线si(n)，i＝1,…,I，n＝1,...,N；I为特征点个数，N为X射线血管造影图像序列中图像的总帧数；
(2)根据平移运动变化帧的序列{Nt,t＝1,....,T}和平移运动的斜率角序列{αt,t＝1,....,T-1}得到模拟平移曲线dα(n)；其中T为运动方向变化的次数；
(3)根据所述X射线血管造影图像序列确定心脏运动信号周期Nc；根据所述跟踪曲线si(n)和所述模拟平移曲线dα(n)得到去除了平移运动后的多运动的综合运动曲线
Figure PCTCN2014085724-appb-000001
根据心脏运动周期Nc对所述综合运动曲线
Figure PCTCN2014085724-appb-000002
进行傅立叶频域滤波处理，得到心脏运动曲线
Figure PCTCN2014085724-appb-000003
(4)根据所述综合运动曲线
Figure PCTCN2014085724-appb-000004
和所述心脏运动曲线
Figure PCTCN2014085724-appb-000005
得到不包含平移运动信号和心脏运动信号的剩余曲线
Figure PCTCN2014085724-appb-000006
在呼吸运动信号周期[3Nc,10Nc]的范围内，根据每一个呼吸周期对所述剩余曲线
Figure PCTCN2014085724-appb-000007
进行傅立叶频域滤波，得到相应的呼吸运动曲线；以最接近所述剩余曲线
Figure PCTCN2014085724-appb-000008
的拟合曲线
Figure PCTCN2014085724-appb-000009
为最优准则，得到相对于当前模拟平移曲线的最佳呼吸运动曲线
Figure PCTCN2014085724-appb-000010
和最佳呼吸运动信号周期
Figure PCTCN2014085724-appb-000011
(5)根据所述剩余曲线
Figure PCTCN2014085724-appb-000012
和所述拟合曲线
Figure PCTCN2014085724-appb-000013
得到的曲线
Figure PCTCN2014085724-appb-000014
判断所述曲线
Figure PCTCN2014085724-appb-000015
的振幅是否小于3个像素，若是，则没有高频运动信号；若否，则在高频运动信号周期[1/7Nc,5/7Nc]的范围内，根据每一个高频运动周期对所述剩余曲线
Figure PCTCN2014085724-appb-000016
进行傅立叶频域滤波，得到相应的高频运动曲线；以最接近所述曲线
Figure PCTCN2014085724-appb-000017
为最优准则，得到相对于当 前模拟平移曲线的最佳高频运动曲线
Figure PCTCN2014085724-appb-000018
和最佳高频运动信号周期
Figure PCTCN2014085724-appb-000019
(6)根据模拟平移曲线dα(n)、呼吸运动曲线
Figure PCTCN2014085724-appb-000020
心脏运动曲线
Figure PCTCN2014085724-appb-000021
和高频运动曲线
Figure PCTCN2014085724-appb-000022
获得综合运动估计曲线
Figure PCTCN2014085724-appb-000023
以所述综合运动估计曲线
Figure PCTCN2014085724-appb-000024
与所述跟踪曲线si(n)最接近为最优准则，得到最佳的平移运动曲线，心脏运动曲线，呼吸运动曲线和高频运动曲线。
总体而言，通过本发明所构思的以上技术方案与现有技术相比，由于考虑到X射线血管造影图像是反映了多种运动信号如平移运动信号、心脏运动信号、呼吸运动信号和高频运动信号等综合，提出了多运动参数的分离估计，能够使分离出来的各个运动更加精确，为医生的诊断提供准确有效的帮助。
[附图说明]
图1是本发明实施例提供的X射线血管造影图像中多运动参数的分离估计方法的实现流程框图；
图2(a)是傅立叶变换提取中造影图及骨骼点的标记,对应的投影角度为(-26.5°，-20.9°)；
图2(b)是傅立叶变换提取中造影图及骨骼点的标记,对应的投影角度为(42.3°，26.8°)；
图3(a)是图2(a)的骨架图及特征点标记。
图3(b)是图2(b)的骨架图及特征点标记。
图4(a)是对图3(a)中所有标记点在造影图序列中跟踪得到的原始曲线，其中，实线为各个特征点的横向坐标(X轴坐标)的变化，虚线为各个特征点的纵向坐标(Y轴坐标)的变化。
图4(b)是对图3(b)中所有标记点在造影图序列中跟踪得到的原始曲线，其中，实线为各个特征点的横向坐标(X轴坐标)的变化，虚线为 各个特征点的纵向坐标(Y轴坐标)的变化。
图4(e)是对图3(a)中原始曲线分解得到的心脏运动曲线，其中，实线为各个特征点的横向坐标(X轴坐标)的变化，虚线为各个特征点的纵向坐标(Y轴坐标)的变化。
图4(f)是对图3(b)中原始曲线分解得到的心脏运动曲线，其中，实线为各个特征点的横向坐标(X轴坐标)的变化，虚线为各个特征点的纵向坐标(Y轴坐标)的变化。
图4(g)是对图3(a)中原始曲线分解得到的呼吸运动信号曲线，其中，实线为各个特征点的横向坐标(X轴坐标)的变化，虚线为各个特征点的纵向坐标(Y轴坐标)的变化。
图4(h)是对图3(b)中原始曲线分解得到的呼吸运动信号曲线，其中，实线为各个特征点的横向坐标(X轴坐标)的变化，虚线为各个特征点的纵向坐标(Y轴坐标)的变化。
图4(i)是对图3(a)中原始曲线分解得到的高频成分运动信号曲线，其中，实线为各个特征点的横向坐标(X轴坐标)的变化，虚线为各个特征点的纵向坐标(Y轴坐标)的变化。
图4(j)是对图3(b)中原始曲线分解得到的高频成分运动信号曲线，其中，实线为各个特征点的横向坐标(X轴坐标)的变化，虚线为各个特征点的纵向坐标(Y轴坐标)的变化。
图4(k)是对图3(a)中原始曲线分解得到的X轴平移运动信号曲线与跟踪图2(a)中的骨骼点得到的X轴曲线的对比，其中，实线为分解得到的X轴平移运动信号曲线，虚线为跟踪骨骼点得到的X轴曲线。
图4(l)是对图3(b)中原始曲线分解得到的X轴平移运动信号曲线与跟踪图2(b)中的骨骼点得到的X轴曲线的对比，其中，实线为分解得到的X轴平移运动信号曲线，虚线为跟踪骨骼点得到的X轴曲线。
图4(m)是对图3(a)中原始曲线分解得到的Y轴平移运动信号曲 线与跟踪图2(a)中的骨骼点得到的Y轴曲线的对比，其中，实线为分解得到的Y轴平移运动信号曲线，虚线为跟踪骨骼点得到的Y轴曲线。
图4(n)是对图3(b)中原始曲线分解得到的Y轴平移运动信号曲线与跟踪图2(b)中的骨骼点得到的Y轴曲线的对比，其中，实线为分解得到的Y轴平移运动信号曲线，虚线为跟踪骨骼点得到的Y轴曲线。
图4(o)是对图3(a)中原始曲线分解得到各个运动得到剩余噪声曲线，其中，实线为各个特征点的横向坐标(X轴坐标)的变化，虚线为各个特征点的纵向坐标(Y轴坐标)的变化。
图4(p)是对图3(b)中原始曲线分解得到各个运动曲线得到剩余噪声曲线，其中，实线为各个特征点的横向坐标(X轴坐标)的变化，虚线为各个特征点的纵向坐标(Y轴坐标)的变化；
图4中横坐标“Image Frame”表示“图像帧数”；纵坐标'displacement/mm'表示“运动位移/mm”。
[具体实施方式]
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。
本发明具体涉及一种血管结构特征点的选取和傅里叶级数展开相结合的数据分离方法。该方法可以解决人体的平移运动信号，心脏运动信号和呼吸运动信号等多种运动的分离难题。
本发明的目的在于提供一种X射线血管造影图像中多运动参数的分离估计方法。本方法基于频域滤波和多变量寻优，能够精确分离各个运动参数，简单快速、无伤害，而且无论单臂还是双臂造影都可以广泛应用。
X射线血管造影图像不仅记录了心脏自身运动在冠状动脉上的投影， 同时还包含了人体的呼吸运动引起的冠状动脉在造影面上的二维运动和病人的二维平移运动。
其步骤包括：
傅里叶级数展开和多变量寻优的方法分离平移信号，心脏信号与呼吸运动信号。具体算法如下：
Step1：对血管结构特征点在X射线血管造影图像序列中进行跟踪，得到点的跟踪曲线si(n)，i＝1,…,I，n＝1,...,N。其中I为特征点个数，N为X射线血管造影图像序列中图像的总帧数。其中，血管结构特征点需要能够真实综合反映出血管整体的运动信息，因此，根据生理学和解剖学等知识的指导，所选的特征点为血管的分叉点。血管结构特征点的得到是首先对X射线血管造影图像进行预处理，通过分割和细化得到血管骨架，然后通过血管骨架选取血管结构特征点。本发明采用的跟踪方法是通过X射线血管造影图像骨架序列不同帧的每个分叉点建立对应关系得到每个特征点在每一帧的坐标。
Step2：设平移运动变化帧的序列{Nt,t＝1,....,T}，平移运动的斜率角序列{αt,t＝1,....,T-1}，其中T为运动方向变化的次数；根据下式可以得到模拟平移曲线dα(n)。
Figure PCTCN2014085724-appb-000025
通过造影图确定X轴平移运动变化帧的序列
Figure PCTCN2014085724-appb-000026
其中Tx为X轴移动方向变换的次数，以斜率角
Figure PCTCN2014085724-appb-000027
范围为(-π/2,π/2)可以通过上述方法得到X轴模拟平移曲线
Figure PCTCN2014085724-appb-000028
同样，通过造影图确定Y轴移 动变化帧的序列
Figure PCTCN2014085724-appb-000029
其中Ty为Y轴移动方向变换的次数，以斜率角
Figure PCTCN2014085724-appb-000030
范围为(-π/2,π/2)可以通过上述方法得到Y轴模拟平移曲线
Figure PCTCN2014085724-appb-000031
使
Figure PCTCN2014085724-appb-000032
表示病人的平移运动。
Step3：根据造影图序列确定心脏运动信号周期Nc。对去除了平移运动后的多运动的综合信号
Figure PCTCN2014085724-appb-000033
即
Figure PCTCN2014085724-appb-000034
根据心脏运动周期Nc进行傅立叶频域滤波处理，得到心脏运动信号
Figure PCTCN2014085724-appb-000035
Step4：对不包含平移运动信号和心脏运动信号的剩余分量
Figure PCTCN2014085724-appb-000036
以呼吸运动信号周期Nαr＝mcNc,mc∈[3,10]的范围，对每一个呼吸周期进行傅立叶频域滤波，以最接近
Figure PCTCN2014085724-appb-000037
的拟合曲线
Figure PCTCN2014085724-appb-000038
为最优准则搜索得到相对于当前平移曲线的最佳呼吸运动信号曲线
Figure PCTCN2014085724-appb-000039
和最佳呼吸运动信号周期
Figure PCTCN2014085724-appb-000040
Step5：对不包含平移运动信号和心脏运动信号的剩余分量
Figure PCTCN2014085724-appb-000041
减去其拟合曲线
Figure PCTCN2014085724-appb-000042
后的曲线
Figure PCTCN2014085724-appb-000043
的振幅若小于3个像素，则可以认为没有高频成分运动信号；否则，以高频成分运动信号周期Nαh＝mhNc,mh∈[1/7,5/7]的范围，对每一个高频成分运动周期进行傅立叶频域滤波，以最接近
Figure PCTCN2014085724-appb-000044
为最优准则搜索得到相对于当前平移曲线的最佳高频成分运动信号曲线
Figure PCTCN2014085724-appb-000045
和最佳高频成分运动信号周期
Figure PCTCN2014085724-appb-000046
Step6：综合运动估计曲线为
Figure PCTCN2014085724-appb-000047
以
Figure PCTCN2014085724-appb-000048
与si(n)最接近为最优准则，得到最佳的曲线平移运动信号曲线，心脏运动曲线，呼吸运动信号曲线和高频成分运动信号。
本发明方法考虑到在实际X射线造影的条件下，并非在所有造影图中都能找到合适的特征点(像肋骨的交点等)，所以采用血管结构特征点的选取和傅立叶频域滤波相以及多变量寻优结合的途径来提取多个运动分量，比单纯的用手动跟踪来得到呼吸运动或者平移运动具有更广泛的适用性和灵活性，几乎可适用于所有造影序列图，同时此方法相较于直接在心脏附 近组织设置标识点再通过相关成像手段进行跟踪的方法拥有更高安全性和可操作性，这是因为在体内组织添加的标记物一般是可侵入性的，会对人体自身产生或多或少的损害，并且其标记物的添加、成像、排除、呼吸运动提取整个过程都是繁杂的，为实际操作中带来不可避免的麻烦与误差。
以下结合附图和实例对本发明进一步说明：
本发明提出了一种直接从X射线造影序列图像中提取平移运动信号，呼吸运动信号和心脏运动信号等多运动参数的方法，通过在冠脉血管上选取一些特征点并在序列中进行跟踪得到它们随时间的运动曲线，然后通过傅里叶级数展开，频域滤波和多变量寻优的方法分离得到单独的平移信号，心脏信号与呼吸运动信号等信号。具有很好的临床适用性。
临床上测得的正常的心脏运动为60～100次/分钟，其周期为0.6-1.0s，而将频率高于140次/分钟的周期运动则诊断为病理性异常运动，周期小于0.42s。如果能够用冠脉造影图分析出运动频率高于正常心脏运动，而且该高频信号的幅度超过2-3像素时，便可以合理的认为该高频信号为某种有用信号(或诊断病变信号，本文称为高频成分)。这是因为，在对冠脉造影图像进行处理的过程中，因分割或骨架提取的不精确性会导致一定的误差，但是这种误差一般不会超过3个像素。所以，当我们分析出存在这种高频运动而且幅度大于3个像素时，就可以为医生提供一种新的有用诊断信号。而呼吸运动信号的频率相比这两种运动则要小的多，呼吸运动信号的周期也就长得多，一般为3-6s，安静的时候可能更长。
根据造影图中各运动成分(包括平移信号、心脏运动信号、呼吸运动信号及高频成分运动信号等)的这些特点，本文通过傅立叶级数展开结合频域滤波的方法对这些成分进行分离。
下面首先介绍一下离散傅立叶级数变换。
1)离散傅立叶级数变换基础
由傅立叶级数展开的原理，连续时间周期信号可以用傅立叶级数来表 达，与此相应，周期序列也可用离散傅立叶级数来表示。若周期序列xp(n)的周期为Np，那么
xp(n)＝xp(n+rNp)(r为任意整数)           (1)
xp(n)可以通过公式(2)和(3)进行离散的傅立叶级数变换
Figure PCTCN2014085724-appb-000049
Figure PCTCN2014085724-appb-000050
对于式(3)作如下解释：式中
Figure PCTCN2014085724-appb-000051
是周期序列的基频分量，
Figure PCTCN2014085724-appb-000052
就是k次谐波分量，各次谐波的系数为Xp(k)；全部谐波分量只有Np个是独立的，因为
Figure PCTCN2014085724-appb-000053
因此，级数取和的项数是k＝0,...,Np-1，共Np个独立谐波分量。而式(2)正是由xp(n)决定系数Xp(k)的求和公式。
2)分离算法
假设造影图图像序列中冠脉血管上第i个特征点Pi(x,y)的X轴坐标的运动曲线为xi(n)(n为造影帧帧数，代表离散化的时间，序列总长度为N)，Y轴坐标的运动曲线为yi(n)，令si(n)＝(xi(n),yi(n))，可将si(n)表示为:
si(n)＝di(n)+ri(n)+ci(n)+hi(n)+ti(n)i∈[1,I]       (4)
其中I表示特征点的个数，d(n)＝(xd(n),yd(n))表示随时间变化的病人的刚性平移运动，ri(n)＝(xri(n),yri(n))表示呼吸运动引起的运动，ci(n)＝(xci(n),yci(n))表示心脏的运动引起的血管上各点的运动，hi(n)＝(xhi(n),yhi(n))表示病人可能具有的高频运动引起的运动，ti(n)＝(xti(n),yti(n))表示其他噪声分量(如因为血管分割、细化等原因造成的)。为方便表示，后面都用si(n)表示提取的血管特征点沿X轴，Y轴的坐标变化曲线，d(n),ri(n),ci(n),hi(n)和ti(n)分别表示由平移运动信号、呼吸运动信号、心脏运动信号、高频成分运动信号以及 噪声分量沿X轴、Y轴的坐标变化曲线。
假设心脏运动信号周期为Nc，呼吸运动信号周期为Nr，高频运动信号周期为Nh，则c(n)、r(n)和h(n)根据式(3)可变为：
Figure PCTCN2014085724-appb-000054
其中，
Figure PCTCN2014085724-appb-000055
是心脏运动信号的谐波分量，
Figure PCTCN2014085724-appb-000056
是呼吸运动信号的谐波分量，
Figure PCTCN2014085724-appb-000057
是高频成分运动信号的谐波分量,C(k)、R(k)和H(k)分别表示c(n)、r(n)和h(n)的各次谐波系数。要使c(n)、r(n)和h(n)完全分离，则其谐波分量不能发生重合(假设此时平移运动d(n)已消除)，即
Figure PCTCN2014085724-appb-000058
要满足式(6)，必须使Nc、Nr和Nh两两互质。同时要分离c(n)、r(n)和h(n)，还必须满足s(n)的序列长度大于或等于s(n)的一个周期，即大于Nc·Nr·Nh。
然而上述要求很难达到，就算Nc、Nr和Nh两两互质，第二个条件也几乎无法满足。因为一般情况下，造影剂在人体内停留的时间不会太长，一个造影图序列只会持续6s左右，不会超过10s，假设帧速为80ms/frame(实际可能更小)，则Nc≈10,Nr>30，Nc·Nr·Nh>Nc·Nr>300>10/0.08，很难获得一个完整周期的序列。而且病人的平移动信号曲线，呼吸运动信号周期确 定的要求较高。
虽然无法直接对心脏动信号、呼吸动信号、高频成分动信号和平移动信号曲线进行理想的分离，但是根据它们各自运动的特点，可以通过傅立叶级数展开，频域滤波和多变量寻优的方法来进行近似的分离。
分离算法基于如下假设：
(1)假设周期为Na,当序列长度为周期的整数倍时且不是其他运动周期的整数倍时，该运动只在频率
Figure PCTCN2014085724-appb-000059
上存在分量，这一点已经通过式(3)证明；而其他运动频率分量在
Figure PCTCN2014085724-appb-000060
处平滑。由于无论是呼吸运动信号，心脏运动信号还是高频成分运动信号都具有平滑性，且互相为对方的基频分量的可能性较小，完全有理由这样假设。
(2)由于平移运动信号是随时间变化的非周期性的刚体运动，所以平移运动信号估计假设为折线运动，如下式：
Figure PCTCN2014085724-appb-000061
其中，{Nt,t＝1,....,T}为平移运动信号变化帧，{αt,t＝1,....,T-1}为平移运动信号的斜率角，其中T为运动方向变化的次数。
如图1所示，为整体流程图。
具体算法如下：
Step1：对血管结构特征点在X射线血管造影图像序列中进行跟踪，得到点的跟踪曲线si(n)，i＝1,…,I，n＝1,...,N。其中I为特征点个数，N为X射线血管造影图像序列中图像的总帧数。其中，血管结构特征点需要能够真实综合反映出血管整体的运动信息，因此，根据生理学和解剖学等知识 的指导，所选的特征点为血管的分叉点。如图2所示，在投影角分别为(-26.5°，-20.9°)和(42.3°，26.8°)的一对造影图中，分别存在着5个编号后的特征点(如图3(a)、图3(b)中的标记点)。血管结构特征点的得到是首先对X射线血管造影图像进行预处理，通过分割和细化得到血管骨架，然后通过血管骨架选取血管结构特征点。本发明采用的跟踪方法是通过X射线血管造影图像骨架序列不同帧的每个分叉点建立对应关系得到每个特征点在每一帧的坐标。
在理想情况下,可将si(n)表示为:
si(n)≈d(n)+ri(n)+ci(n)+hi(n)i∈[1,I]        (8)
Step2:平移运动信号曲线估计；
若存在病人平移，则通过造影图确定X轴移动变化帧的序列
Figure PCTCN2014085724-appb-000062
其中Tx为X轴移动方向变换的次数，以及Y轴移动变化帧的序列
Figure PCTCN2014085724-appb-000063
其中Ty为Y轴移动方向变换的次数。并以斜率角
Figure PCTCN2014085724-appb-000064
范围为(-π/2,π/2)分别得到X轴模拟平移曲线
Figure PCTCN2014085724-appb-000065
和Y轴模拟平移曲线
Figure PCTCN2014085724-appb-000066
如公式(7)。使
Figure PCTCN2014085724-appb-000067
表示病人的平移运动信号。
去除了平移运动后的多运动的综合信号为
Figure PCTCN2014085724-appb-000068
Step3:心脏运动曲线估计
根据造影图序列确定心脏运动信号周期Nc，一般情况下，Nc≈0.8/f，其中，Nc表示在一个心动周期内的造影图帧数，0.8表示一个心动周期的时间，单位为秒，f表示帧频，即每秒钟造影图被拍摄的数量。通过对
Figure PCTCN2014085724-appb-000069
选取长度为nsc的序列
Figure PCTCN2014085724-appb-000070
nsc为Nc的整数倍。将
Figure PCTCN2014085724-appb-000071
分别分解为X方向的运动曲线xic(n)与Y方向的运动曲线yic(n)，再分别对它们进行离散傅里叶展开，得到各次谐波系数Xic(k)和Yic(k)，k＝0,...,nsc-1；
不包含平移运动信号和心脏运动信号的剩余运动信号曲线的频域响应为
Figure PCTCN2014085724-appb-000072
心脏运动的频域响应为：
Figure PCTCN2014085724-appb-000073
对
Figure PCTCN2014085724-appb-000074
进行离散傅里叶逆变换得到X轴心脏运动曲线
Figure PCTCN2014085724-appb-000075
同理可得
Figure PCTCN2014085724-appb-000076
即心脏运动曲线
Figure PCTCN2014085724-appb-000077
对不包含平移运动信号和心脏运动信号的剩余运动曲线进行傅里叶逆变换可得
Figure PCTCN2014085724-appb-000078
Step4:呼吸运动信号曲线估计
对呼吸运动信号周期Nαr＝mcNc,mc∈[3,10]，当原始运动信号序列长于呼吸运动信号周期Nαr，则通过对
Figure PCTCN2014085724-appb-000079
进行与心脏运动信号分离相同的方法进行傅立叶频域滤波得到的呼吸运动信号曲线
Figure PCTCN2014085724-appb-000080
当呼吸运动信号周期长于原始运动信号序列，则不包含呼吸运动信号的一个周期，则用以下步骤分离呼吸运动信号：
a.从序列
Figure PCTCN2014085724-appb-000081
第j点开始取长度为Nαr/2的序列
Figure PCTCN2014085724-appb-000082
对其进行傅里叶离散展开得到各次谐波系数
Figure PCTCN2014085724-appb-000083
b.对于长度为呼吸周期信号一半的曲线将认为其包含了周期为Nαr/2的运动曲线。该运动只在频率
Figure PCTCN2014085724-appb-000084
上存在分量。通过傅立叶频域滤波得到该运动曲线
Figure PCTCN2014085724-appb-000085
c.将
Figure PCTCN2014085724-appb-000086
取均值得到近似呼吸运动信号曲线
Figure PCTCN2014085724-appb-000087
然后，对
Figure PCTCN2014085724-appb-000088
进行曲线拟合得到呼吸运动信号估计曲线
Figure PCTCN2014085724-appb-000089
(因为呼 吸运动信号的周期最长，同时比高频成分运动信号的值要大得多，所以该曲线是呼吸运动信号的大致形状)。通过下式(10)为最优准则得到最佳呼吸运动信号周期
Figure PCTCN2014085724-appb-000090
以及最佳呼吸运动信号曲线
Figure PCTCN2014085724-appb-000091
Figure PCTCN2014085724-appb-000092
Step5:高频成分运动信号曲线估计
对不包含平移运动信号和心脏运动信号的剩余分量
Figure PCTCN2014085724-appb-000093
减去其拟合曲线
Figure PCTCN2014085724-appb-000094
后的曲线
Figure PCTCN2014085724-appb-000095
即
Figure PCTCN2014085724-appb-000096
若
Figure PCTCN2014085724-appb-000097
振幅不超过3个像素，则认为不具有高频成分(因为在对冠脉造影图进行处理的过程中，因分割或骨架提取的不精确性导致的误差一般不会超过2-3个像素)。否则，对于高频成分运动信号曲线周期Nαh＝mhNc,mh∈[1/7,5/7]，通过对
Figure PCTCN2014085724-appb-000098
与心脏运动信号分离相同的方法进行傅立叶频域滤波得到高频成分运动信号曲线
Figure PCTCN2014085724-appb-000099
通过下式(11)为最优准则得到最佳高频成分运动信号周期
Figure PCTCN2014085724-appb-000100
以及最佳高频成分运动信号曲线
Figure PCTCN2014085724-appb-000101
Figure PCTCN2014085724-appb-000102
Step6:运动估计曲线为
Figure PCTCN2014085724-appb-000103
通过(12)式为最优准则，得到最佳平移运动信号曲线
Figure PCTCN2014085724-appb-000104
以及相应的心脏运动曲线
Figure PCTCN2014085724-appb-000105
呼吸运动信号曲线
Figure PCTCN2014085724-appb-000106
和高频成分运动信号曲线
Figure PCTCN2014085724-appb-000107
Figure PCTCN2014085724-appb-000108
通过对造影图序列的观察确定心脏运动周期为11frames，图4是对图3中标记点的运动进行运动分离的结果。
从图4中可以看到，特征点的原始运动曲线受呼吸运动影响较大，比 较凌乱，而且不同特征点的运动幅度不同，这是因为不同特征点位于心脏的不同部位，而心脏不同部位的运动情况是不一样的。分离后的心脏运动曲线显示了良好的周期性，将所有标记点的呼吸运动曲线在一起进行比较，可以发现同一个造影面内所有特征点的呼吸运动曲线相差不大(如图4(g)和图4(h))。图4(i)和图4(j)为分离得到的高频成分大约是心脏运动的0.54倍，周期约为6frames，且其幅度基本都超过3个像素，若是因为分割或骨架提取等导致的误差应该是在3个像素范围内，如图4(k)，图4(l)，图4(m)和图4(n)，因此可以假定它是某种有规律的心脏异常运动(如心脏震颤等)。图4(k)，图4(l)，图4(m)和图4(n)为分离的平移曲线与手动跟踪的平移曲线的对比，可以看出两者除了因为手动跟踪引起的误差外，基本一致。图4(o)和图4(p)为分离各个运动的剩余曲线，其幅度都少于3个像素，可以看作是因为分割，骨架提取等算法误差引起的，运动分量基本完全分离。
本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (1)

  1. 一种X射线血管造影图像中多运动参数的分离估计方法，其特征在于，包括下述步骤：
    (1)对血管结构特征点在X射线血管造影图像序列中进行跟踪，得到点的跟踪曲线si(n)，i＝1,…,I，n＝1,...,N；I为特征点个数，N为X射线血管造影图像序列中图像的总帧数；
    (2)根据平移运动变化帧的序列{Nt,t＝1,....,T}和平移运动的斜率角序列{αt,t＝1,....,T-1}得到模拟平移曲线dα(n)；其中T为运动方向变化的次数；
    (3)根据所述X射线血管造影图像序列确定心脏运动信号周期Nc；根据所述跟踪曲线si(n)和所述模拟平移曲线dα(n)得到去除了平移运动后的多运动的综合运动曲线
    Figure PCTCN2014085724-appb-100001
    根据心脏运动周期Nc对所述综合运动曲线
    Figure PCTCN2014085724-appb-100002
    进行傅立叶频域滤波处理，得到心脏运动曲线
    Figure PCTCN2014085724-appb-100003
    (4)根据所述综合运动曲线
    Figure PCTCN2014085724-appb-100004
    和所述心脏运动曲线
    Figure PCTCN2014085724-appb-100005
    得到不包含平移运动信号和心脏运动信号的剩余曲线
    Figure PCTCN2014085724-appb-100006
    在呼吸运动信号周期[3Nc,10Nc]的范围内，根据每一个呼吸周期对所述剩余曲线
    Figure PCTCN2014085724-appb-100007
    进行傅立叶频域滤波处理，得到相应的呼吸运动曲线；以最接近所述剩余曲线
    Figure PCTCN2014085724-appb-100008
    的拟合曲线
    Figure PCTCN2014085724-appb-100009
    为最优准则，得到相对于当前模拟平移曲线的最佳呼吸运动曲线
    Figure PCTCN2014085724-appb-100010
    和最佳呼吸运动信号周期
    Figure PCTCN2014085724-appb-100011
    (5)根据所述剩余曲线
    Figure PCTCN2014085724-appb-100012
    和所述拟合曲线
    Figure PCTCN2014085724-appb-100013
    得到的曲线
    Figure PCTCN2014085724-appb-100014
    判断所述曲线
    Figure PCTCN2014085724-appb-100015
    的振幅是否小于3个像素，若是，则没有高频运动信号；若否，则在高频运动信号周期[1/7Nc,5/7Nc]的范围内，根据每一个高频运动周期对所述剩余曲线
    Figure PCTCN2014085724-appb-100016
    进行傅立叶频域滤波处理，得到相应的高频运动曲线；以最接近所述曲线
    Figure PCTCN2014085724-appb-100017
    为最优准则，得到相对于当前模拟平移曲线的最佳高频运动曲线
    Figure PCTCN2014085724-appb-100018
    和最佳高频运动信号周期
    Figure PCTCN2014085724-appb-100019
    (6)根据模拟平移曲线dα(n)、呼吸运动曲线
    Figure PCTCN2014085724-appb-100020
    心脏运动曲线
    Figure PCTCN2014085724-appb-100021
    和高频运动曲线
    Figure PCTCN2014085724-appb-100022
    获得综合运动估计曲线
    Figure PCTCN2014085724-appb-100023
    以所述综合运动估计曲线
    Figure PCTCN2014085724-appb-100024
    与所述跟踪曲线si(n)最接近为最优准则，得到最佳的平移运动曲线，心脏运动曲线，呼吸运动曲线和高频运动曲线。
PCT/CN2014/085724 2013-12-31 2014-09-02 一种x射线血管造影图像中多运动参数的分离估计方法 Ceased WO2015101059A1 (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
US14/698,894 US9286676B2 (en) 2013-12-31 2015-04-29 Method for separating and estimating multiple motion parameters in X-ray angiogram image

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN201310752751.4A CN103886615B (zh) 2013-12-31 2013-12-31 一种x射线血管造影图像中多运动参数的分离估计方法
CN201310752751.4 2013-12-31

Related Child Applications (1)

Application Number Title Priority Date Filing Date
US14/698,894 Continuation-In-Part US9286676B2 (en) 2013-12-31 2015-04-29 Method for separating and estimating multiple motion parameters in X-ray angiogram image

Publications (1)

Publication Number Publication Date
WO2015101059A1 true WO2015101059A1 (zh) 2015-07-09

Family

ID=50955488

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2014/085724 Ceased WO2015101059A1 (zh) 2013-12-31 2014-09-02 一种x射线血管造影图像中多运动参数的分离估计方法

Country Status (3)

Country Link
US (1) US9286676B2 (zh)
CN (1) CN103886615B (zh)
WO (1) WO2015101059A1 (zh)

Families Citing this family (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103886615B (zh) * 2013-12-31 2016-04-20 华中科技大学 一种x射线血管造影图像中多运动参数的分离估计方法
CN104517301B (zh) * 2014-12-30 2017-07-07 华中科技大学 多参数模型指导的迭代提取血管造影图像运动参数的方法
GB2545641B (en) * 2015-12-15 2020-05-27 Siemens Medical Solutions Usa Inc A method for detecting motion in a series of image data frames, and providing a corresponding warning to a user
CN105631864B (zh) * 2015-12-23 2018-06-12 电子科技大学 一种基于时空相关性的心脏表面目标点运动预测方法
US10447972B2 (en) * 2016-07-28 2019-10-15 Chigru Innovations (OPC) Private Limited Infant monitoring system
DE102017211677A1 (de) * 2017-07-07 2019-01-10 Siemens Healthcare Gmbh Bewegungsabhängige Rekonstruktion von Magnetresonanzabbildungen
US10499867B2 (en) * 2018-01-08 2019-12-10 Shenzhen Keya Medical Technology Corporation Method, storage medium, and system for analyzing image sequences of periodic physiological activities
CN108245178A (zh) * 2018-01-11 2018-07-06 苏州润迈德医疗科技有限公司 一种基于x射线冠脉造影图像的血液流动速度计算方法
CN109557509B (zh) * 2018-11-23 2020-08-11 安徽四创电子股份有限公司 一种用于改善脉间干扰的双脉冲信号合成器
JP2020141841A (ja) * 2019-03-06 2020-09-10 コニカミノルタ株式会社 動態解析装置及びプログラム
CN109978915A (zh) * 2019-03-11 2019-07-05 北京理工大学 X光造影图像序列中管状结构的跟踪方法及装置
US11357573B2 (en) * 2019-04-25 2022-06-14 International Business Machines Corporation Optimum treatment planning during coronary intervention by simultaneous simulation of a continuum of outcomes
TWI864040B (zh) * 2019-07-04 2024-12-01 新加坡商派拉梅維爾私人有限公司 診斷支援系統

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2001045374A (ja) * 1999-07-29 2001-02-16 Shimadzu Corp デジタルサブトラクション装置
US20080240536A1 (en) * 2007-03-27 2008-10-02 Elisabeth Soubelet Method of detection and compensation for respiratory motion in radiography cardiac images synchronized with an electrocardiogram signal
CN101283910A (zh) * 2008-06-05 2008-10-15 华北电力大学 一种获取冠状动脉血管运动信息的方法
CN101773395A (zh) * 2009-12-31 2010-07-14 华中科技大学 一种从单臂x射线造影图中提取呼吸运动参数的方法
US20120189183A1 (en) * 2010-09-28 2012-07-26 Siemens Corporation Automatic Registration of Image Series with Varying Contrast Based on Synthetic Images Derived from Intensity Behavior Model
WO2012167616A1 (en) * 2011-06-09 2012-12-13 The Hong Kong University Of Science And Technology Image based tracking
CN103886615A (zh) * 2013-12-31 2014-06-25 华中科技大学 一种x射线血管造影图像中多运动参数的分离估计方法

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO1998056294A1 (en) * 1997-06-12 1998-12-17 The Johns Hopkins University School Of Medicine Method for measuring myocardial motion and the like
PT1485007E (pt) * 2002-03-20 2007-10-11 Novadaq Technologies Inc ''sistema e método para visualização do fluxo de fluido através de vasos''
DE102004059182A1 (de) * 2004-12-08 2006-06-14 Siemens Ag Betriebsverfahren für einen Rechner und hiermit korrespondierende Einrichtungen
US9412044B2 (en) * 2009-06-09 2016-08-09 Siemens Aktiengesellschaft Method of compensation of respiratory motion in cardiac imaging
US8224056B2 (en) * 2009-12-15 2012-07-17 General Electronic Company Method for computed tomography motion estimation and compensation
US9129424B2 (en) * 2012-04-17 2015-09-08 Siemens Aktiengesellschaft Phase sensitive T1 mapping in magnetic resonance imaging

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2001045374A (ja) * 1999-07-29 2001-02-16 Shimadzu Corp デジタルサブトラクション装置
US20080240536A1 (en) * 2007-03-27 2008-10-02 Elisabeth Soubelet Method of detection and compensation for respiratory motion in radiography cardiac images synchronized with an electrocardiogram signal
CN101283910A (zh) * 2008-06-05 2008-10-15 华北电力大学 一种获取冠状动脉血管运动信息的方法
CN101773395A (zh) * 2009-12-31 2010-07-14 华中科技大学 一种从单臂x射线造影图中提取呼吸运动参数的方法
US20120189183A1 (en) * 2010-09-28 2012-07-26 Siemens Corporation Automatic Registration of Image Series with Varying Contrast Based on Synthetic Images Derived from Intensity Behavior Model
WO2012167616A1 (en) * 2011-06-09 2012-12-13 The Hong Kong University Of Science And Technology Image based tracking
CN103886615A (zh) * 2013-12-31 2014-06-25 华中科技大学 一种x射线血管造影图像中多运动参数的分离估计方法

Also Published As

Publication number Publication date
US9286676B2 (en) 2016-03-15
US20150317793A1 (en) 2015-11-05
CN103886615A (zh) 2014-06-25
CN103886615B (zh) 2016-04-20

Similar Documents

Publication Publication Date Title
CN103886615B (zh) 一种x射线血管造影图像中多运动参数的分离估计方法
CN101773395B (zh) 一种从单臂x射线造影图中提取呼吸运动参数的方法
JP6293495B2 (ja) 動いている器官の標的エリア内のオブジェクトを追跡するための方法および装置
CN101809618B (zh) 对介入工具的探测和跟踪
Shechter et al. Prospective motion correction of X-ray images for coronary interventions
Ferrario et al. Toward morphological thoracic EIT: major signal sources correspond to respective organ locations in CT
WO2018044549A1 (en) Respiration motion stabilization for lung magnetic navigation system
CN104519793A (zh) 用于内腔的线性映射的系统
US12458397B2 (en) Method and device for registration and tracking during a percutaneous procedure
JPWO2019135412A1 (ja) 診断支援プログラム
JP2012205899A (ja) ３次元的モデルを利用した身体臓器の映像生成方法及び装置並びにコンピュータ読み取り可能な記録媒体
JP2011516143A (ja) 呼吸測定器
CN108294768A (zh) 序列图像多参数配准的x射线心血管造影减影方法及系统
WO2015101060A1 (zh) 单臂x射线血管造影图像多运动参数分解估计方法
CN104517301B (zh) 多参数模型指导的迭代提取血管造影图像运动参数的方法
JP2008526420A (ja) 画像の位置合わせのための画像処理システム及び方法
CN111862312B (zh) 一种脑部血管显示装置及方法
CN114469146B (zh) X射线图像中的冠脉运动补偿方法和系统
CN112089438A (zh) 基于二维超声图像的四维重建方法及装置
Zhang et al. A novel structural features-based approach to automatically extract multiple motion parameters from single-arm X-ray angiography
WO2025088466A1 (en) Target movement modelling using electromagnetic navigation sensors
TWI828661B (zh) 診斷支援系統
CN117058191B (zh) 用于胸腹部肿瘤运动估计的肝脏超声血管跟踪方法及设备
JP4991181B2 (ja) ３次元画像処理装置、３次元画像処理方法及び３次元画像処理装置で使用される制御プログラム
CN115089294B (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: 14876077

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

Country of ref document: EP

Kind code of ref document: A1