CN101773395A - Method for extracting respiratory movement parameter from one-arm X-ray radiography picture - Google Patents
Method for extracting respiratory movement parameter from one-arm X-ray radiography picture Download PDFInfo
- Publication number
- CN101773395A CN101773395A CN200910273528A CN200910273528A CN101773395A CN 101773395 A CN101773395 A CN 101773395A CN 200910273528 A CN200910273528 A CN 200910273528A CN 200910273528 A CN200910273528 A CN 200910273528A CN 101773395 A CN101773395 A CN 101773395A
- Authority
- CN
- China
- Prior art keywords
- mrow
- msub
- msubsup
- curve
- prime
- 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.)
- Granted
Links
- 230000033001 locomotion Effects 0.000 title claims abstract description 188
- 230000000241 respiratory effect Effects 0.000 title claims abstract description 125
- 238000000034 method Methods 0.000 title claims abstract description 23
- 238000002601 radiography Methods 0.000 title claims abstract description 15
- 238000001914 filtration Methods 0.000 claims abstract description 8
- 230000000747 cardiac effect Effects 0.000 claims description 19
- 210000004351 coronary vessel Anatomy 0.000 claims description 15
- 238000000605 extraction Methods 0.000 claims description 10
- 230000003434 inspiratory effect Effects 0.000 claims description 6
- 238000002059 diagnostic imaging Methods 0.000 abstract description 2
- 210000001367 artery Anatomy 0.000 abstract 1
- 230000000007 visual effect Effects 0.000 abstract 1
- 210000004204 blood vessel Anatomy 0.000 description 12
- 238000003384 imaging method Methods 0.000 description 9
- 239000003550 marker Substances 0.000 description 7
- 230000029058 respiratory gaseous exchange Effects 0.000 description 7
- 230000000737 periodic effect Effects 0.000 description 6
- 238000000926 separation method Methods 0.000 description 5
- 230000002792 vascular Effects 0.000 description 4
- 210000000056 organ Anatomy 0.000 description 3
- 210000003484 anatomy Anatomy 0.000 description 2
- 238000006073 displacement reaction Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 230000009466 transformation Effects 0.000 description 2
- 210000000038 chest Anatomy 0.000 description 1
- 239000002872 contrast media Substances 0.000 description 1
- 230000002950 deficient Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 210000004072 lung Anatomy 0.000 description 1
- 230000035479 physiological effects, processes and functions Effects 0.000 description 1
- 230000001360 synchronised effect Effects 0.000 description 1
- 238000012800 visualization Methods 0.000 description 1
Images
Landscapes
- Apparatus For Radiation Diagnosis (AREA)
Abstract
The invention relates to a method for extracting a respiratory movement parameter from a one-arm X-ray radiography picture, which belongs to the crossed field of digital signal processing and medical imaging, and aims at separating the respiratory movement from the complex movement of an actual human body heart. The method comprises the following steps: selecting a group of architectural feature points (including the head and tail points of a vessel section) on a coronal artery vessel and tracking the motion curve of the architectural feature point along with the time in a sequence; and separating to obtain a bi-dimensional respiratory movement curve by the Fourier expansion and frequency-domain filtering method. In addition, in the invention, in order to describe the respiratory movement more intuitively and effectively, the bi-dimensional respiratory movement obtained from two different radiography visual angles is rebuilt according to the three-dimensional reconstruction theory of the mid point of a radiography system, therefore the three-dimensional respiratory movement is obtained. The invention has the advantages of better extracting result of the respiratory movement and wider applicability and flexibility, and meets clinical requirements.
Description
Technical Field
The invention belongs to the field of digital signal processing and medical imaging intersection, and particularly relates to a method for extracting respiratory motion parameters from a single-arm x-ray radiography image.
Background
The thorax rhythmically expands and contracts, thereby completing inspiration and expiration, which is the breathing movement. Respiratory motion causes overall translational motion of the human heart in three dimensions. In an X-ray imaging system, the coronary vessels are subject to two-dimensional translational motion on the imaging plane due to the influence of respiratory motion.
The coronary images thus record, on the one hand, a projection of the movement of the heart onto a two-dimensional plane, and at the same time superimpose a two-dimensional translational movement of the coronary arteries on the image plane, which is caused by the respiratory movement of the body. Then, to obtain a two-dimensional angiogram closer to the real situation and use it for three-dimensional reconstruction of blood vessels, it is necessary to separate the two kinds of motion information so as to analyze the cardiac motion and respiratory motion of the human body separately.
Many foreign documents set marking points in advance when extracting human respiratory motion and then carry out sequence tracking on the marking points. In x-ray radiography, it is often used to directly and manually track non-cardiac structural feature points in the angiogram. According to the characteristics of breathing movement, people can drive other organs in the body to move together when breathing. It is believed that these organs translate in three dimensions as the lungs move, and that their movements are synchronized. It can be assumed that the motion of the heart caused by the respiratory motion and the motion of the organs adjacent to it are also identical in the contrast image plane, and some characteristic points on other tissues outside the heart can be found in the contrast image as marker points. The marked points are tracked in the whole sequence, the motion situation of the marked points is obtained, and then the motion of the marked points is approximated to the breathing motion on the two-dimensional projection surface. Another approach also makes use of the structural feature points which do not move with the heart, except that the movement of the marker points is recorded simultaneously with the imaging. This therefore requires that the individual feature points be selected and marked before the visualization. Both of these solutions are clearly deficient. The former is mainly due to its poor applicability, because we cannot guarantee that there are some marker points (some feature points on other tissues outside the heart) meeting this condition in each frame of the angiogram, and it also needs experience to find such points (it needs to know the human anatomy relatively). Therefore, when the above feature points do not exist in the histogram, the respiratory motion is difficult to extract. The latter implementation requires a lot of experimental control and is not suitable for general clinical application.
At present, a method for obtaining respiratory motion parameters under the condition of double-arm x-ray radiography also appears, and the idea of separating cardiac motion and respiratory motion is to take two radiography images at the same moment and different projection angles, and perform three-dimensional reconstruction on corresponding coronary vessels to obtain the three-dimensional spatial distribution of the vessels at the moment. Then, after all pairs of the contrast maps in a respiratory cycle are matched and reconstructed, a set of three-dimensional structural sequences is obtained, and the spatial displacement vector between the three-dimensional structural sequences is respiratory motion. Relatively reliable respiratory motion estimation results can be obtained by this method, but cannot be widely applied in medical practice due to the constraints of two-arm x-ray imaging conditions.
Disclosure of Invention
The invention aims to provide a method for extracting respiratory motion parameters from a single-arm X-ray contrast image sequence, which comprises the steps of selecting a group of characteristic points on coronary vessels, tracking the characteristic points in the sequence to obtain motion curves of the characteristic points along with time, and separating by Fourier series expansion and frequency domain filtering to obtain independent heart and respiratory motion information. Has good clinical applicability.
A human body respiratory motion parameter extraction method based on frequency domain filtering comprises the following steps:
(1) acquiring a single-arm X-ray radiography image sequence of coronary vessels, and determining a heart motion period N1;
(2) selecting coronary vessel structure characteristic points;
(3) extracting a respiratory motion curve in a specific mode:
(3.1) automatically tracking the coronary vessel structure characteristic points in an angiogram image sequence to obtain a characteristic point tracking sequence s (n), wherein n is the length of the characteristic point tracking sequence;
(3.2) selecting one section in the characteristic point tracking sequence s (n) as a target sequence
ns-N% N1, where% is the remainder symbol;
(3.3) to the target sequence
Performing discrete Fourier transform to obtain a target frequency domain response S (k), wherein k is 0.
(3.4) frequency-Domain response to respiratory motion
By the respiratory motion extraction method, a two-dimensional respiratory motion curve under a certain contrast angle is finally obtained. However, in order to describe the respiratory motion more intuitively and more effectively, the two-dimensional motion obtained under two different contrast view angles can be reconstructed according to the three-dimensional reconstruction principle of the midpoint of the contrast system, so that the three-dimensional respiratory motion is obtained.
The method comprises the following specific steps:
step 1: determining two contrast image sequences under different contrast angles, recording as a left contrast image sequence and a right contrast image sequence, determining any corresponding reference point in the left contrast image sequence and the right contrast image sequence, respectively recording as pl(x, y, t) and pr(x, y, t), where (x, y) is the coordinate of any one of the reference points in the image sequence, t is the time, i.e. the frame number of the contrast image sequence, and t is the frame number of the contrast image sequenceIntegers between 0 and 70, respectively, at time tl，trCoordinates of lower (x)l，yl)、(xr，yr) As reference point pl(x, y, t) and prInitialization value of (x, y, t), i.e. pl(xl，yl，tl) And pr(xr，yr，tr) Wherein t isl，trSelecting end diastole corresponding to the same moment of the cardiac cycle;
step 2: extracting respiratory motion curves of left and right contrast image sequences according to the respiratory motion parameter extraction method, and respectively recording the respiratory motion curves as curvel(x, y, t) and curer(x, y, t), then correspondingly selecting an extreme point of an inspiratory end or an expiratory end in the respiratory motion period as a respiratory reference point on the two curves respectively, namely selecting an extreme point of the inspiratory end or selecting an extreme point of the expiratory end on the two curves respectively, and marking as curvel(xcl，ycl，tcl) And curver(xcr，ycr，tcr) Wherein t iscl，tcrThe time of the end of inspiration or expiration corresponding to the left and right contrast images in the sequence, (x)cl，ycr)、(xcr，ycr) Are respectively curvedl(x, y, t) and curerAt (x, y, t) corresponding to time tcl，tcrThe coordinates of (a);
step 3: for initialized point pairs p in left and right contrast image sequencel(xl，yl，tl) And pr(xr，yr，tr) Performing respiratory motion compensation, i.e. tl，trThe respiratory motion at the moment is respectively compensated to the corresponding respiratory reference moment tcl，tcrThe following steps:
in the formula, curvel(xil，ytl，tl)、curver(xtr，ytr，tr) Respectively represent tl，trCurve at timel(x, y, t) and curerPoint on (x, y, t), p'l(x′l，y′l，t′l) And p'r(x′r，y′r，tr) Is a compensated reference point. Then, three-dimensional reconstruction is carried out on the two compensated reference points to obtain a three-dimensional point P (x)clr，yclr，zclr，tclr) Wherein t isclrRepresents the sum of tcl(or t)cr) Corresponding to a consistent end of inspiration or end of expiration;
step 4: it is reasonable to assume that the heart itself is stationary and that it is only the respiration that causes the left coronary vascular tree to move. Then p 'obtained in Step3 can be obtained'l(x′l，y′l，tl) And p'r(x′r，y′r，tr) Further compensation to other moments of the respiratory motion cycle without the influence of cardiac motion:
similarly, point pairs p 'for other times of the respiratory motion cycle'l(x′l+i，y′l+i，tl+i) And p'r(x′r+i，y′r+i，tr+i) Performing three-dimensional reconstruction to obtain a three-dimensional point P (x)clr+i，yclr+i，zclr+i，tclr+i) Where i denotes the time t relative to the reference instant of respirationlAnd trThe number of frames before (negative integer) or after (positive integer);
step 5: and combining the results obtained by the Step3 and the Step4 to obtain the three-dimensional respiratory motion.
The technical effects of the invention are as follows: considering that under the condition of single-arm X-ray radiography, proper characteristic points (such as rib intersection points and the like) can not be found in all radiography images, a new way of combining the selected vascular structure characteristic points and Fourier frequency domain filtering is provided to extract respiratory motion, the method has wider applicability and flexibility compared with the method of pure manual tracking, and can be almost suitable for all radiography sequence images (the requirement of clear blood vessel distribution and two or more cardiac cycles). Meanwhile, compared with a method for directly arranging the identification points on the tissues near the heart and then tracking the tissues through a related imaging means, the method has higher safety and operability. This is because the markers added to the body tissue are generally invasive, causing more or less damage to the body itself, and the whole process of adding, imaging, removing, and extracting respiratory motion is complicated, which brings inevitable troubles and errors to the actual operation. In addition, the characteristic points selected by the method relate to all levels of blood vessels of the left coronary artery, and the motion information of the left coronary artery is comprehensively considered, so that the method has better reliability and accuracy.
Drawings
FIG. 1 is a block flow diagram of the present invention;
FIG. 2(a) is a contrast image and its labeled points selected in the respiratory motion extraction method, corresponding to projection angles (-26.8 °, -27.2 °);
FIG. 2(b) is a contrast image and its labeled points selected in the respiratory motion extraction method, corresponding to projection angles (50.8 degrees, 30.2 degrees);
FIG. 3(a) is a raw curve traced in a sequence of angiograms for marker point 1 in FIG. 2(a), where the solid line is the change in the lateral coordinate (X-axis coordinate) and the dashed line is the change in the longitudinal coordinate (Y-axis coordinate); FIG. 3(b) is a raw curve traced in a sequence of angiograms for marker point 1 in FIG. 2(b), where the solid line is the change in the lateral coordinate (X-axis coordinate) and the dashed line is the change in the longitudinal coordinate (Y-axis coordinate);
FIG. 3(c) is a heart motion curve decomposed from the original curve in FIG. 3(a), in which the solid line represents the variation of the horizontal coordinate (X-axis coordinate) and the dotted line represents the variation of the vertical coordinate (Y-axis coordinate);
FIG. 3(d) is a heart motion curve decomposed from the original curve in FIG. 3(b), in which the solid line represents the variation of the horizontal coordinate (X-axis coordinate) and the dotted line represents the variation of the vertical coordinate (Y-axis coordinate);
FIG. 3(e) is a respiratory motion curve decomposed from the original curve in FIG. 3(a), wherein the solid line represents the variation of the horizontal coordinate (X-axis coordinate) and the dotted line represents the variation of the vertical coordinate (Y-axis coordinate);
FIG. 3(f) is a respiratory motion curve decomposed from the original curve in FIG. 3(b), wherein the solid line represents the variation of the horizontal coordinate (X-axis coordinate) and the dotted line represents the variation of the vertical coordinate (Y-axis coordinate);
FIG. 3(g) is the result of 6 curve fits to the respiratory motion curve of FIG. 3(e), where the solid line is the change in the lateral coordinate (X-axis coordinate) and the dashed line is the change in the longitudinal coordinate (Y-axis coordinate);
FIG. 3(h) is the result of 6 curve fits to the respiratory motion curve of FIG. 3(f), where the solid line is the change in the lateral coordinate (X-axis coordinate) and the dashed line is the change in the longitudinal coordinate (Y-axis coordinate);
FIG. 4(a) is a plot of respiratory motion for all marker points at projection angles (-26.8 deg., -27.2 deg.), with the solid line being the change in lateral coordinate (X-axis coordinate) and the dashed line being the change in longitudinal coordinate (Y-axis coordinate);
fig. 4(b) is a respiratory motion curve for all marked points at a projection angle (50.8 °, 30.2 °), where the solid line is the change in the lateral coordinate (X-axis coordinate) and the dotted line is the change in the longitudinal coordinate (Y-axis coordinate);
FIG. 4(c) is a result of fitting a motion curve for all points in FIG. 4(a), where the solid line is the change in the lateral coordinate (X-axis coordinate) and the dashed line is the change in the longitudinal coordinate (Y-axis coordinate);
FIG. 4(d) is the result of fitting a motion curve for all points in FIG. 4(b), where the solid line is the change in the lateral coordinate (X-axis coordinate) and the dashed line is the change in the longitudinal coordinate (Y-axis coordinate);
fig. 5 is a three-dimensional result obtained after reconstructing two-dimensional respiratory motion at contrast angles (-26.8 deg., -27.2 deg.) and (50.8 deg., 30.2 deg.).
Detailed Description
The invention is further described below with reference to the accompanying drawings:
the invention provides a human body respiratory motion parameter extraction method based on frequency domain filtering, which comprises the following steps as shown in figure 1:
(1) acquiring a single-arm X-ray radiography image sequence of coronary vessels, and determining a heart motion period N1;
(2) and selecting the characteristic points of the vascular structure.
According to the guidance of knowledge of physiology, anatomy and the like, the marked feature points need to comprehensively reflect the motion information of the whole blood vessel, so that the selected special points (i.e. structural feature points) generally comprise the starting point and the ending point of each blood vessel segment and each inflection point (the point with the maximum curvature) between the blood vessel segments. In the sequence of the cartographic images under two different projection angles, all the feature points are numbered, and the corresponding feature points have the same name. As shown in fig. 2, in a pair of the contrast maps having projection angles of (-26.8 °, -27.2 °) and (50.8 °, 30.2 °), respectively, there are 15 numbered distinctive points, and the relationship thereof is one-to-one according to the numerical designation.
(3) Fourier series expansion method for separating heart and respiratory motion
The heart and the human respiratory motion are periodic motion, but the frequency of the respiratory motion is much smaller than that of the heart motion, generally speaking, the frequency of the normal motion of the heart is 60-100 times/min, the period of the normal motion is 0.6-1.0s, and the period of the respiratory motion is much longer, generally 3-6s, and the rest time can be longer. On the other hand, the heart motion is relatively severe, and the amplitude of the respiratory motion is small, i.e. the resulting displacement is small, which is a relatively smooth process.
According to the characteristics of the heart motion and the respiratory motion in the contrast image, the invention separates the heart motion from the respiratory motion by combining Fourier series expansion with frequency domain filtering.
The following first describes the discrete fourier series expansion.
1) Discrete Fourier series expansion
By the principle of fourier series expansion, a continuous-time periodic signal can be expressed by a fourier series, and correspondingly, a periodic sequence can also be expressed by a discrete fourier series. If periodic sequence xp(N) has a period of N, then
xp(n)＝xp(n + rN) (r is an arbitrary integer) (1)
xp(n) discrete Fourier series transformation by equations (2) and (3)
The following explanation is made for formula (3): in the formula
Is the fundamental frequency component of the periodic sequence,is the k harmonic component, Xp(k) Is the coefficient of each harmonic, which can also be considered as xp(n) a frequency domain response; only N of all harmonic components are independent, because
Thus, the number of terms for the sum of the series is N, from k-0 to N-1Independent harmonic components. And the formula (2) is just formed by xp(n) determining the coefficient Xp(k) The summation formula of (1).
2) Separation algorithm
Assuming that a motion curve of x-axis coordinates of a certain point P (x, y) on a coronary vessel in an angiogram image sequence is x (n) (n is the number of frames of an angiogram frame), and a motion curve of y-axis coordinates is y (n), let s (n) be (x (n), y (n)), s (n) be decomposed into the following equation:
s(n)＝c(n)+r(n)+t(n)
wherein c (n) ═ xc(n)，yc(n)) represents the motion of a blood vessel point caused by the motion of the heart, and r (n) ═ xr(n)，yr(n)) represents the motion caused by respiratory motion, and t (n) ═ xt(n)，yt(n)) represents some other motion, including the movement of the body of the person during the imaging process, and the movement of the imaging equipment, etc. For convenience, s (n) is used to represent the coordinate variation curves of the extracted blood vessel points along the x-axis and the y-axis, c (n) is used to represent the coordinate variation curves of the blood vessel points along the x-axis and the y-axis caused by heart motion, r (n) is used to represent the coordinate variation curves of the blood vessel points along the x-axis and the y-axis caused by respiratory motion, so that the operation on s (n) is the operation on x (n) and y (n), and the operation on c (n) is the operation on x (n), respectivelyc(n) and yc(n) the operations on r (n) are respectively on xr(n) and yr(n) operation. Since t (n) is not measurable, neglecting it from consideration here, there are
s(n)≈c(n)+r(n)
Assuming that the period of cardiac motion is N1 and the period of respiratory motion is N2, c (N) and r (N) may become according to equation (3):
wherein,
are the harmonic components of the motion of the heart,harmonic components of respiratory motion, c (k), r (k) represent the respective harmonic coefficients of c (n) and r (n), respectively. To completely separate c (n) and r (n), the harmonic components cannot be coincident (zero frequency components can be disregarded because although both the frequency of the respiratory motion and the frequency domain components of the cardiac motion are present above zero frequency (dc components), there is no way to completely separate them, but since the zero frequency component varies to a constant in the time domain, it does not change the shape of the periodic sequence curve, nor does it affect our analysis of the cardiac and respiratory motion), i.e. the zero frequency component is a constant in the time domain
(8)
To satisfy equation (8), N1 and N2 must be made relatively prime. To separate c (N) and r (N), it is necessary to satisfy that the sequence length of s (N) is greater than or equal to one period of s (N), i.e. greater than N1 × N2.
However, the above requirements are difficult to achieve, and even if N1 and N2 are relatively prime, the second condition cannot be satisfied. Since the time during which the contrast agent stays in the body is generally not too long, a sequence of contrast images lasts only about 6s, not more than 10s, and assuming a frame rate of 80ms/frame (which may actually be smaller), N1 ≈ 10, N2 > 30, N1 × N2 > 300 > 10/0.08, a sequence of one cycle is not obtained at all.
Although the ideal separation of cardiac and respiratory motion is not possible, the approximate separation can be performed by fourier series transformation and frequency domain filtering, depending on the characteristics of their respective motions.
The separation step is based on the following assumptions:
(1) heart movement only at frequency
There is a component. This is demonstrated by formula (3);
(2) frequency domain component of respiratory motion in frequency
Is smoothed. Due to the smooth breathing motion, an
This is entirely reasonable for the low probability of the fundamental frequency component of respiratory motion.
The specific algorithm is as follows:
step 1: determining a cardiac motion cycle N1 from a sequence of human eye-observed angiograms, typically N1 ≈ 0.8/f, where N1 denotes the number of angiogram frames in one cardiac cycle, 0.8 denotes the time of one cardiac cycle in seconds, and f denotes the frame rate, i.e. the number of angiograms captured per second;
step 2: tracking the blood vessel structure characteristic points selected in the step (1) in the whole contrast map sequence to obtain a tracking sequence s (n) of the points;
step 3: selecting a sequence with the length of ns from the tracking sequence s (n) of the point
ns is an integer multiple of N1. If the original sequence length is N, selecting N1 with the sequence length ns being N-N%, wherein,% is a remainder symbol;
step 4: will be provided withIs decomposed intoMotion x (n) in the x direction and motion y (n) in the y direction, and then performing discrete fourier transform on the motion x (n) and the motion y (n) in the y direction to obtain subharmonic coefficients x (k) and y (k), wherein k is 0.
Step 5: frequency domain response of respiratory motion
By analyzing the sequence of contrast maps, the cardiac motion period N1 can be determined to be 10 frames, and fig. 3 is the result of separating the respiratory motion from the cardiac motion for the motion of marker point 1 in fig. 2.
As can be seen from fig. 3, the original motion curve of the marked points is greatly influenced by the respiratory motion and is relatively disordered, while the separated heart motion curve shows good periodicity, but the respiratory motion curve still has some burrs, which is also the reason for incomplete separation. In order to remove the residual heart motion component in the respiratory motion, curve fitting is performed on the extracted original respiratory motion curve to remove the burrs on the respiratory motion curve, and as shown in fig. 3(g) and 3(h), the fitted respiratory motion curve shows a distinct periodicity.
Comparing the respiratory motion curves of all the marked points together, it can be found that the respiratory motion curves of all the points in the same contrast surface are not very different (as shown in fig. 4(a) and (b)), so we fit the respiratory motion curves of all the marked points into a curve as the respiratory motion curve of the coronary vessel in the contrast surface (as shown in fig. 4(c) and (d)).
(3) Three-dimensional reconstruction of respiratory motion
By the above-described respiratory motion extraction method, a two-dimensional respiratory motion curve at a certain contrast angle is finally obtained, as shown in fig. 4(c) and 4 (d). However, it can be seen from the figure that the two-dimensional representation method obviously cannot clearly express respiratory information, so that in order to describe respiratory motion more intuitively and more effectively, two-dimensional motion obtained under two different contrast angles can be reconstructed according to a three-dimensional reconstruction principle of a midpoint of an imaging system, so as to obtain three-dimensional respiratory motion. The method comprises the following specific steps:
step 1: determining two contrast image sequences under different contrast angles, recording as a left contrast image sequence and a right contrast image sequence, determining any corresponding reference point in the left contrast image sequence and the right contrast image sequence, respectively recording as pl(x, y, t) and pr(x, y, t), where (x, y) is the coordinate of any reference point in the image sequence, t is time, i.e. the frame number of the contrast image sequence, and t is an integer between 0 and 70, the time t in t is selected respectivelyl，trCoordinates of lower (x)l，yl)、(xr，yr) As reference point pl(x, y, t) and prInitialization value of (x, y, t), i.e. pl(xl，yl，tl) And pr(xr，yr，tr) Wherein t isl，trSelecting end diastole corresponding to the same moment of the cardiac cycle;
step 2: extracting respiratory motion curves of left and right contrast image sequences according to the respiratory motion parameter extraction method, and respectively recording the respiratory motion curves as curvel(x, y, t) and curer(x, y, t), then correspondingly selecting an extreme point of an inspiratory end or an expiratory end in the respiratory motion period as a respiratory reference point on the two curves respectively, namely selecting an extreme point of the inspiratory end or selecting an extreme point of the expiratory end on the two curves respectively, and marking as curvel(xcl，ycl，tcl) And curver(xcr，ycr，tcr) Wherein t iscl，tcrThe time of the end of inspiration or expiration corresponding to the left and right contrast images in the sequence, (x)cl，ycl)、(xcr，ycr) Are respectively curvedl(x, y, t) and curerAt (x, y, t) corresponding to time tcl，tcrThe coordinates of (a);
step 3: for initialized point pairs p in left and right contrast image sequencel(xl，yl，tl) And pr(xr，yr，tr) Performing respiratory motion compensation, i.e. tl，trThe respiratory motion at the moment is respectively compensated to the corresponding respiratory reference moment tcl，tcrThe following steps:
in the formula, curvel(xil，ytl，tl)、curver(xtr，ytr，tr) Respectively represent tl，trCurve at timel(x, y, t) and curerPoint on (x, y, t), p'l(x′l，y′l，tl) And p'r(x′r，y′r，tr) Is a compensated reference point. Then, three-dimensional reconstruction is carried out on the two compensated reference points to obtain a three-dimensional point P (x)clr，yclr，zclr，tclr) Wherein t isclrRepresents the sum of tcl(or t)cr) Corresponding to a consistent end of inspiration or end of expiration;
step 4: it is reasonable to assume that the heart itself is stationary and that it is only the respiration that causes the left coronary vascular tree to move. Then p 'obtained in Step3 can be obtained'l(x′l，y′l，tl) And p'r(x′r，y′r，tr) Further compensation to other moments of the respiratory motion cycle without the influence of cardiac motion:
similarly, point pairs p 'for other times of the respiratory motion cycle'l(x′l+i，y′l+i，tl+i) And p'r(x′r+i，y′r+i，tr+i) Performing three-dimensional reconstruction to obtain a three-dimensional point P (x)clr+i，yclr+i，zclr+i，tclr+i) Where i denotes the time t relative to the reference instant of respirationlAnd trThe number of frames before (negative integer) or after (positive integer);
step 5: and combining the results obtained by the Step3 and the Step4 to obtain the three-dimensional respiratory motion.
By analyzing the sequence of the contrast images under the angles of (-26.8 degrees, -27.2 degrees) (left) and (50.8 degrees, 30.2 degrees (right), respectively, the 22 nd frame and the 43 th frame which are at the end diastole are taken as reference frames, respectively, the blood vessel starting points thereon are taken as reference points and are initialized to (35, 62) and (303, 94), and then the respiratory motion curves extracted from the sequence of the left contrast image and the right contrast image are analyzed, and respectively the respiratory reference time is determined to be the 38 th frame and the 26 th frame. Then, from the parameters selected above, a three-dimensional respiratory motion can be found, as shown in fig. 5.
As can be seen from the figure, the overall trend of the three-dimensional respiratory motion is to make a reciprocating translational motion, but a motion bevel still exists locally, and the motion bevel is not simple linear motion. It has been analyzed that the incidence of this dog-ear phenomenon should be caused by the heart motion effects not being completely separated from the respiratory motion. However, the motion deviation caused by the folding angle has little influence on the whole, so that the three-dimensional respiratory motion can be reasonably assumed to be the piecewise linear motion under the real condition.
Claims (2)
1. A human body respiratory motion parameter extraction method based on frequency domain filtering comprises the following steps:
(1) acquiring a single-arm X-ray radiography image sequence of coronary vessels, and determining a heart motion period N1;
(2) selecting coronary vessel structure characteristic points;
(3) extracting a respiratory motion curve in a specific mode:
(3.1) automatically tracking the coronary vessel structure characteristic points in the image sequence of the contrast map to obtain a characteristic point tracking sequence s (n), wherein n is characteristic point trackingThe length of the sequence; (3.2) selecting one section in the characteristic point tracking sequence s (n) as a target sequence
ns-N% N1, where% is the remainder symbol;
(3.3) to the target sequence
Performing discrete Fourier transform to obtain a target frequency domain response S (k), wherein k is 0.
(3.4) frequency-Domain response to respiratory motion
And performing inverse discrete Fourier transform on the R (k) to obtain a respiratory motion curve r (ns), namely, extracting respiratory motion parameters.
2. A method for three-dimensional reconstruction of the respiratory motion curve of claim 1 to obtain three-dimensional respiratory motion, comprising the steps of:
(I) determining two contrast image sequences under different contrast angles, recording as a left contrast image sequence and a right contrast image sequence, determining any corresponding reference point in the left contrast image sequence and the right contrast image sequence, respectively recording as pl(x, y, t) and pr(x, y, t), where (x, y) is the coordinate of the reference point in the image sequence, t is time, i.e. the frame number of the contrast image sequence, and t is an integer between 0 and 70, the time t in t is selected respectivelyl，trCoordinates of lower (x)l，yl)、(xr，yr) As reference point pl(x, y, t) and prInitialization value of (x, y, t), i.e. pl(xl，yl，tl) And pr(xr，yr，tr) Wherein t isl，trThe same time corresponding to the cardiac cycle in the left and right angiogram sequences;
(II) extracting respiratory motion curves of left and right contrast image sequences according to the respiratory motion parameter extraction method of claim 1, and respectively recording the curves as curvel(x, y, t) and curer(x, y, t), then correspondingly selecting an extreme point of an inspiratory end or an expiratory end in the respiratory motion period as a respiratory reference point on the two curves respectively, namely selecting an extreme point of the inspiratory end or selecting an extreme point of the expiratory end on the two curves respectively, and marking as curvel(xcl，ycl，tcl) And curver(xcr，ycr，tcr) Wherein t iscl，tcrThe time of the end of inspiration or expiration corresponding to the left and right contrast images in the sequence, (x)cl，ycl)、(xcr，ycr) Are respectively curvedl(x, y, t) and curerAt (x, y, t) corresponding to time tcl、tcrThe coordinates of (a);
(III) for initialized point pairs p in left and right contrast image sequencesl(xl，yl，tl) And pr(xr，yr，tr) Performing respiratory motion compensation, i.e. tl，trThe respiratory motion at the moment is respectively compensated to the corresponding respiratory reference moment tcl，tcrThe following steps:
in the formula, curvel(xtl，ytl，tl)、curver(xtr，ytr，tr) Respectively represent tl，trCurve at timel(x, y, t) and curerPoint on (x, y, t), p'l(x′l，y′l，tl) And p'r(x′r，y′r，tr) Is a compensated reference point. Then, three-dimensional reconstruction is carried out on the two compensated reference points to obtain a three-dimensional point P (x)clr，yclr，zclr，tclr) Wherein t isclrRepresents the sum of tclOr tcrCorresponding to a consistent end of inspiration or end of expiration;
(IV) reacting p 'obtained in step (III)'l(x′l，y′l，tl) And p'r(x′r，y′r，tr) Further compensation to other moments of the respiratory motion cycle without the influence of cardiac motion:
similarly, point pairs p 'for other times of the respiratory motion cycle'l(x′l+i，y′l+i，tl+i) And p'r(x′r+i，y′r+i，tr+i) Performing three-dimensional reconstruction to obtain a three-dimensional point P (x)clr+i，yclr+i，zclr+i，tclr+i) Where i denotes at the time tlAnd trThe number of previous or subsequent frames, i is a negative integer when the number of previous frames is the number of previous frames, and i is a positive integer when the number of subsequent frames is the number of subsequent frames;
and (V) integrating the results obtained in the step (III) and the step (IV) to obtain the three-dimensional respiratory motion.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN2009102735285A CN101773395B (en) | 2009-12-31 | 2009-12-31 | Method for extracting respiratory movement parameter from one-arm X-ray radiography picture |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN2009102735285A CN101773395B (en) | 2009-12-31 | 2009-12-31 | Method for extracting respiratory movement parameter from one-arm X-ray radiography picture |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN101773395A true CN101773395A (en) | 2010-07-14 |
| CN101773395B CN101773395B (en) | 2011-05-18 |
Family
ID=42510080
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN2009102735285A Expired - Fee Related CN101773395B (en) | 2009-12-31 | 2009-12-31 | Method for extracting respiratory movement parameter from one-arm X-ray radiography picture |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN101773395B (en) |
Cited By (13)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103083023A (en) * | 2011-10-31 | 2013-05-08 | Ge医疗系统环球技术有限公司 | Method and device for measuring respiratory cycle through computerized tomography (CT) scanner |
| CN103083030A (en) * | 2011-10-31 | 2013-05-08 | Ge医疗系统环球技术有限公司 | Device-less 4 dimensional-computed tomography (D4D-CT) imaging method, device and system |
| CN103340619A (en) * | 2013-06-18 | 2013-10-09 | 浙江大学 | Electromagnetic wave physiological movement imaging system |
| CN103810721A (en) * | 2013-12-30 | 2014-05-21 | 华中科技大学 | Single-arm x-ray angiography image multiple motion parameter decomposition and estimation method |
| CN103830848A (en) * | 2012-11-27 | 2014-06-04 | Ge医疗系统环球技术有限公司 | Method and system for gating radiotherapy |
| CN104517301A (en) * | 2014-12-30 | 2015-04-15 | 华中科技大学 | Method for iteratively extracting movement parameters of angiography image guided by multi-parameter model |
| WO2015101059A1 (en) * | 2013-12-31 | 2015-07-09 | 华中科技大学 | Separation and estimation method for multiple motion parameters in x-ray angiographic image |
| CN106056589A (en) * | 2016-05-24 | 2016-10-26 | 西安交通大学 | Ultrasound contrast perfusion parameter imaging method based on respiratory motion compensation |
| CN107505584A (en) * | 2016-06-14 | 2017-12-22 | 西门子(深圳)磁共振有限公司 | A kind of magnetic resonance data acquisition triggering method and device |
| CN107569251A (en) * | 2017-08-29 | 2018-01-12 | 上海联影医疗科技有限公司 | Medical imaging procedure and system and non-transient computer readable storage medium storing program for executing |
| CN108245178A (en) * | 2018-01-11 | 2018-07-06 | 苏州润迈德医疗科技有限公司 | A kind of blood flowing speed computational methods based on X ray coronary angiography image |
| CN108364290A (en) * | 2018-01-08 | 2018-08-03 | 深圳科亚医疗科技有限公司 | Method, medium and the system that the image sequence of cyclical physiological activity is analyzed |
| CN108470357A (en) * | 2018-03-27 | 2018-08-31 | 中科超精（安徽）科技有限公司 | Couple the elastic registrating method of respiratory phase |
-
2009
- 2009-12-31 CN CN2009102735285A patent/CN101773395B/en not_active Expired - Fee Related
Cited By (21)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103083030B (en) * | 2011-10-31 | 2017-04-12 | Ge医疗系统环球技术有限公司 | Device-less 4 dimensional-computed tomography (D4D-CT) imaging method, device and system |
| CN103083030A (en) * | 2011-10-31 | 2013-05-08 | Ge医疗系统环球技术有限公司 | Device-less 4 dimensional-computed tomography (D4D-CT) imaging method, device and system |
| CN103083023A (en) * | 2011-10-31 | 2013-05-08 | Ge医疗系统环球技术有限公司 | Method and device for measuring respiratory cycle through computerized tomography (CT) scanner |
| CN103830848A (en) * | 2012-11-27 | 2014-06-04 | Ge医疗系统环球技术有限公司 | Method and system for gating radiotherapy |
| CN103830848B (en) * | 2012-11-27 | 2018-09-14 | Ge医疗系统环球技术有限公司 | Method and system for gating radiotherapy |
| CN103340619A (en) * | 2013-06-18 | 2013-10-09 | 浙江大学 | Electromagnetic wave physiological movement imaging system |
| CN103810721A (en) * | 2013-12-30 | 2014-05-21 | 华中科技大学 | Single-arm x-ray angiography image multiple motion parameter decomposition and estimation method |
| WO2015101060A1 (en) * | 2013-12-30 | 2015-07-09 | 华中科技大学 | Decomposition and estimation method for multiple motion parameters in single-arm x-ray angiographic image |
| WO2015101059A1 (en) * | 2013-12-31 | 2015-07-09 | 华中科技大学 | Separation and estimation method for multiple motion parameters in x-ray angiographic image |
| CN104517301A (en) * | 2014-12-30 | 2015-04-15 | 华中科技大学 | Method for iteratively extracting movement parameters of angiography image guided by multi-parameter model |
| CN104517301B (en) * | 2014-12-30 | 2017-07-07 | 华中科技大学 | The method of the iterative extraction angiographic image kinematic parameter that multi-parameters model is instructed |
| WO2016106959A1 (en) * | 2014-12-30 | 2016-07-07 | 华中科技大学 | Multi-parameter model guided-method for iteratively extracting movement parameter of angiography image of blood vessel |
| CN106056589A (en) * | 2016-05-24 | 2016-10-26 | 西安交通大学 | Ultrasound contrast perfusion parameter imaging method based on respiratory motion compensation |
| CN106056589B (en) * | 2016-05-24 | 2018-12-07 | 西安交通大学 | A kind of ultrasonic contrast perfusion parametric imaging method of respiration motion compensation |
| CN107505584A (en) * | 2016-06-14 | 2017-12-22 | 西门子(深圳)磁共振有限公司 | A kind of magnetic resonance data acquisition triggering method and device |
| CN107569251A (en) * | 2017-08-29 | 2018-01-12 | 上海联影医疗科技有限公司 | Medical imaging procedure and system and non-transient computer readable storage medium storing program for executing |
| CN108364290A (en) * | 2018-01-08 | 2018-08-03 | 深圳科亚医疗科技有限公司 | Method, medium and the system that the image sequence of cyclical physiological activity is analyzed |
| CN108364290B (en) * | 2018-01-08 | 2020-10-09 | 深圳科亚医疗科技有限公司 | Method, medium, and system for analyzing a sequence of images of periodic physiological activity |
| CN108245178A (en) * | 2018-01-11 | 2018-07-06 | 苏州润迈德医疗科技有限公司 | A kind of blood flowing speed computational methods based on X ray coronary angiography image |
| CN108470357A (en) * | 2018-03-27 | 2018-08-31 | 中科超精（安徽）科技有限公司 | Couple the elastic registrating method of respiratory phase |
| CN108470357B (en) * | 2018-03-27 | 2022-04-01 | 中科超精(南京)科技有限公司 | Elastic registration method for coupling respiratory phases |
Also Published As
| Publication number | Publication date |
|---|---|
| CN101773395B (en) | 2011-05-18 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN101773395B (en) | Method for extracting respiratory movement parameter from one-arm X-ray radiography picture | |
| CN103886615B (en) | The separating estimation method of parameter of doing more physical exercises in a kind of X ray angiographic image | |
| JP5067398B2 (en) | Ultrasound image and CT image matching system | |
| US20100061611A1 (en) | Co-registration of coronary artery computed tomography and fluoroscopic sequence | |
| US9173626B2 (en) | Method for performing dynamic registration, overlays, and 3D views with fluoroscopic images | |
| JP2015515913A (en) | Method for linear mapping of lumens | |
| JP2016507304A (en) | System for detecting and tracking and superimposing objects | |
| US20120087563A1 (en) | Method and System for Intraoperative Guidance Using Physiological Image Fusion | |
| JP2014140742A (en) | Method and apparatus for tracking object in target area of moving organ | |
| US20110245651A1 (en) | Medical image playback device and method, as well as program | |
| JP2018508276A (en) | Tracking-based 3D model enhancements | |
| KR101768526B1 (en) | Method and apparatus for creating model of patient specified target organ based on blood vessel structure | |
| CN101107628A (en) | Image processing system and method for alignment of images | |
| Lee et al. | Synthesis of electrocardiogram V-lead signals from limb-lead measurement using R-peak aligned generative adversarial network | |
| CN108294768B (en) | X-ray angiocardiography subtraction method and system based on sequence image multi-parameter registration | |
| WO2015101060A1 (en) | Decomposition and estimation method for multiple motion parameters in single-arm x-ray angiographic image | |
| JP6472606B2 (en) | X-ray diagnostic equipment | |
| WO2016106959A1 (en) | Multi-parameter model guided-method for iteratively extracting movement parameter of angiography image of blood vessel | |
| Gao et al. | Intraoperative registration of preoperative 4D cardiac anatomy with real-time MR images | |
| Fischer et al. | An MR-based model for cardio-respiratory motion compensation of overlays in X-ray fluoroscopy | |
| CN112089438B (en) | Four-dimensional reconstruction method and device based on two-dimensional ultrasonic image | |
| Metz et al. | Alignment of 4D coronary CTA with monoplane X-ray angiography | |
| US20100030572A1 (en) | Temporal registration of medical data | |
| Ma et al. | Real-time registration of 3D echo to x-ray fluoroscopy based on cascading classifiers and image registration | |
| Jeong et al. | Deep-learning-based registration of diagnostic angiogram and live fluoroscopy for percutaneous coronary intervention |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| C14 | Grant of patent or utility model | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20110518 |
|
| CF01 | Termination of patent right due to non-payment of annual fee |