CN101773395A - Method for extracting respiratory movement parameter from one-arm X-ray radiography picture - Google Patents

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

A kind of method of from one-arm X-ray radiography picture, extracting respiratory movement parameter

Technical field

The invention belongs to Digital Signal Processing and the field that medical imaging intersects, be specifically related to a kind of method of from single armed x ray contrast figure, extracting respiratory movement parameter.

Background technology

The rhythmic expansion of thorax and dwindling, thus finish air-breathing with exhale respiratory movement that Here it is.The respiratory movement meeting causes the translational motion of human heart integral body in three dimensions.In the x-ray imaging system, because respirometric influence, the two-dimension translational motion can take place in coronary artery on the radiography face.

Therefore, the coronary angiography image records the projection of motion on two dimensional surface of heart on the one hand, also is superimposed with simultaneously the two-dimension translational motion of coronary artery on the radiography face that the respiratory movement of human body causes.So, obtain more near the two-dimentional angiogram under the truth and be used for blood vessel 3 D reconstructing, the heart movement and the respiratory movement of human body be analyzed separately thereby then need these two kinds of movable informations are separated.

External a lot of document all will carry out sequential tracks to them by setting in advance gauge point then when extracting the human body respiration motion.Wherein in the x ray contrast, use many be directly manually to follow the tracks of the architectural feature point of the non-heart among the radiography figure.According to respirometric characteristics, the people can drive intravital other organs and move together when breathing.It is generally acknowledged that these organs can carry out three-dimensional translation along with the motion of lung, and their motion all is synchronous.So the motion that can suppose the heart that respiratory movement causes also is consistent on the radiography plan with the motion of the organ adjacent with it, can find other the more structural characteristic points outside the heart to serve as a mark a little in radiography figure.In whole sequence, follow the tracks of these gauge points, obtain the motion conditions of these gauge points, then with the respiratory movement on the approximate two-dimensional projection's face for this reason of the motion of these gauge points.What another kind of way was utilized equally is that these can not follow the architectural feature point that heart moves together, and different is that it is the motion of writing down these gauge points in radiography.Therefore, this just requires just each characteristic point to be chosen and labelling before radiography.Obviously, these two kinds of schemes all are defective.The former is that mainly its suitability is very poor, because we can not guarantee to exist among each frame radiography figure the gauge point (other more structural characteristic points that heart is outer) that meets this condition, and look for this point also to need experience (need know quite well human anatomic structure).So when not having above characteristic point in radiography figure, respiratory movement is difficult to be extracted out.The latter's realization then needs a large amount of experiment controls, and is improper to general clinical applications.

At present, a kind of employing realizes obtaining respiratory movement parameter under both arms x ray contrast condition method has also appearred, its isolating cardiac motion with respirometric thought is, get two width of cloth radiography figure of the different projection angles of synchronization, corresponding arteria coronaria blood vessel is wherein carried out three-dimensional reconstruction, obtain the blood vessel three-dimensional spatial distribution in this moment.So, to after mating and rebuilding, what obtain then is one group of three dimensional structure sequence to all the radiography figure in the breathing cycle, and the space displacement vector between them is respiratory movement.Comparatively speaking, can obtain reliable respiratory movement estimated result, still, can not be applied in the medical practice widely owing to the constraint of both arms x ray contrast condition by this method.

Summary of the invention

The objective of the invention is to propose a kind of method of from the one-arm X-ray radiography picture image sequence, extracting respiratory movement parameter, obtain their curve movements in time by on the arteria coronaria blood vessel, choosing a group of feature point and in sequence, following the tracks of, separate independent heart and the respiratory movement information of obtaining by the method for Fourier expansion and frequency domain filtering then.Has the good clinical suitability.

A kind of human body respiration kinematic parameter extracting method based on frequency domain filtering comprises the steps:

(1) obtain the one-arm X-ray radiography picture image sequence of arteria coronaria blood vessel, and definite heart movement cycle N1;

(2) choose arteria coronaria blood vessel structure characteristic point;

(3) extract the respiratory movement curve, concrete mode is:

(3.1) in radiography figure image sequence described arteria coronaria blood vessel structure characteristic point is carried out from motion tracking, obtain feature point tracking sequence s (n), n is the feature point tracking sequence length;

(3.2) in feature point tracking sequence s (n), choose one section as target sequence

Ns=n-n%N1, wherein, % is for getting surplus symbol;

(3.3) to target sequence

Carry out discrete Fourier transform (DFT), obtain target frequency domain response S (k), k=0 ..., ns-1;

(3.4) make respirometric frequency domain response

$R\left(k\right)=\left\{\begin{array}{cc}S\left(k\right),& (k\≠\frac{\mathrm{ns}}{N1}l,l=1,...,N1-1)\\ (S\left(k\right)+S(k+1))/2,& (k=\frac{\mathrm{ns}}{N1}l,l=1,...,N1-1)\end{array}\right.,$ R (k) is carried out inverse discrete Fourier transformer inverse-discrete obtain respiratory movement curve r (ns), promptly can extract respiratory movement parameter.

By respiratory movement extracting method recited above, what obtain at last is two-dimentional respiratory movement curve under certain radiography angle.But in order more intuitively, more effectively to describe respiratory movement, we can come the two dimensional motion that obtains under two different radiographies visual angles is rebuild according to the three-dimensional reconstruction principle of radiography system mid point, and obtain three-dimensional respiratory movement.

Concrete steps are as follows:

Step1: determine two radiography figure image sequences under the different radiography angles, be designated as left radiography figure image sequence and right radiography figure image sequence, in the radiography figure image sequence of the left and right sides, determine respectively more corresponding arbitrary reference point to be designated as p respectively _l(x, y, t) and p _r(x, y, t), wherein (x is the coordinate of this arbitrary reference point in image sequence y), and t is the time, i.e. the frame number of radiography figure image sequence, and t is the integer between 0 to 70, chooses the moment t among the t respectively _l, t _rUnder coordinate (x _l, y _l), (x _r, y _r) as reference point p _l(x, y, t) and p _r(x, y, initialization value t), i.e. p _l(x _l, y _l, t _l) and p _r(x _r, y _r, t _r), t wherein _l, t _rCorresponding to the synchronization of cardiac cycle, generally select diastole latter stage;

Step2: extract the respiratory movement curve of left and right sides radiography figure image sequence according to above-mentioned respiratory movement parameter extracting method, be designated as curve respectively _l(x, y, t) and curve _r(x, y, t), then respectively on two curves corresponding the extreme point of air-breathing latter stage in the respiration motion cycle or EEP of selecting as breathing reference point, promptly on two curves, all select extreme point or all select the EEP extreme point in air-breathing latter stage, be designated as curve _l(x _Cl, y _Cl, t _Cl) and curve _r(x _Cr, y _Cr, t _Cr), t wherein _Cl, t _CrBe respectively the air-breathing moment in latter stage or the EEP moment corresponding in the radiography figure image sequence of the left and right sides, (x _Cl, y _Cr), (x _Cr, y _Cr) be respectively curve curve _l(x, y, t) and curve _r(x, y t) go up corresponding t constantly _Cl, t _CrCoordinate;

Step3: to the point after the initialization in the radiography figure image sequence of the left and right sides to p _l(x _l, y _l, t _l) and p _r(x _r, y _r, t _r) carry out respiration motion compensation, be about to t _l, t _rThe time respiratory movement of inscribing compensate to corresponding breathing respectively with reference to moment t _Cl, t _CrDown:

$\left\{\begin{array}{c}{p}_l^{\′}(x_l^{\′},y_l^{\′},t_l)=p_l(x_l,y_l,t_l)-({\mathrm{curve}}_l(x_{\mathrm{tl}},y_{\mathrm{tl}},t_l)-{\mathrm{curve}}_l(x_{\mathrm{cl}},y_{\mathrm{cl}},t_{\mathrm{cl}}))\\ p_r^{\′}(x_r^{\′},y_r^{\′},t_r)=p_r(x_r,y_r,t_r)-({\mathrm{curve}}_r(x_{\mathrm{tr}},y_{\mathrm{tr}},t_r)-{\mathrm{curve}}_r(x_{\mathrm{cr}},y_{\mathrm{cr}},t_{\mathrm{cr}}))\end{array}\right.$

In the formula, curve _l(x _Il, y _Tl, t _l), curve _r(x _Tr, y _Tr, t _r) represent t respectively _l, t _rMoment curve curve _l(x, y, t) and curve _r(x, y, the t) point on, p ' _l(x ' _l, y ' _l, t ' _l) and p ' _r(x ' _r, y ' _r, t _r) for compensating the back reference point.Then the reference point after two compensation is carried out three-dimensional reconstruction, obtain three-dimensional point P (x _Clr, y _Clr, z _Clr, t _Clr), t wherein _ClrExpression and t _Cl(or t _Cr) corresponding consistent air-breathing latter stage or EEP;

Step4: the reasonable assumption heart itself is immobilized, and what cause the motion of left coronary artery vascular tree only is Repiration.So, can be with the p ' that obtains among the Step3 _l(x ' _l, y ' _l, t _l) and p ' _r(x ' _r, y ' _r, t _r) under the influence that does not have heart movement, further compensate to respiration motion cycle other constantly:

$\left\{\begin{array}{c}{p}_l^{\′}(x_{l+i}^{\′},y_{l+i}^{\′},t_{l+i})=p_l^{\′}(x_l^{\′},y_l^{\′},t_l)+({\mathrm{curve}}_l(x_{l+i},y_{l+i},t_{l+i})-{\mathrm{curve}}_l(x_{\mathrm{cl}},y_{\mathrm{cl}},t_{\mathrm{cl}}))\\ p_r^{\′}(x_{r+i}^{\′},y_{r+i}^{\′},t_{r+i})=p_r^{\′}(x_r^{\′},y_r^{\′},t_r)+({\mathrm{curve}}_r(x_{r+i},y_{r+i},t_{r+i})-{\mathrm{curve}}_r(x_{\mathrm{cr}},y_{\mathrm{cr}},t_{\mathrm{cr}}))\end{array}\right.$

In like manner, to other points constantly of respiration motion cycle to p ' _l(x ' _L+i, y ' _L+i, t _L+i) and p ' _r(x ' _R+i, y ' _R+i, t _R+i) carry out three-dimensional reconstruction, obtain three-dimensional point P (x _Clr+i, y _Clr+i, z _Clr+i, t _Clr+i), wherein i represents with respect to breathing with reference to moment t _lAnd t _rThe frame number of preceding (getting negative integer) or back (getting positive integer);

Step5: the result that comprehensive Step3 and Step4 ask for can obtain three-dimensional respiratory movement.

Technique effect of the present invention is embodied in: consider under the condition of one-arm X-ray radiography, be not in all radiography figure, can both find suitable feature point (riblike intersection point etc.), extract respiratory movement so proposed to choose the new way that blood vessel structure characteristic point and Fourier domain filtering combines, have the suitability widely and motility than simple manual tracking, almost applicable to all radiography sequence chart (need satisfying have more clearly vascularity and comprise two or more cardiac cycles).Simultaneously the method compared to directly near heart tissue the method that identification point follows the tracks of by the dependent imaging means again be set have greater security and operability.This is because the label of tissue interpolation generally is can be invasive in vivo, can be to human body self generation infringement more or less, and it all is numerous and diverse that whole process is extracted in the interpolation of its label, imaging, eliminating, respiratory movement, for bringing inevitable trouble and error in the practical operation.In addition, the characteristic point that the inventive method is chosen relates to the blood vessels at different levels of left coronary artery, has taken all factors into consideration the movable information of left coronary artery, thereby has better reliability and accuracy.

Description of drawings

Fig. 1 is a FB(flow block) of the present invention;

Fig. 2 (a) is radiography figure and the gauge point thereof that is selected in the respiratory movement extracting method, and corresponding projection angle is (26.8 ° ,-27.2 °);

Fig. 2 (b) is radiography figure and the gauge point thereof that is selected in the respiratory movement extracting method, and corresponding projection angle is (50.8 °, 30.2 °);

Fig. 3 (a) follows the tracks of the primitive curve that obtains to gauge point 1 among Fig. 2 (a) in the radiography graphic sequence, wherein solid line is the variation of lateral coordinates (X-axis coordinate), and dotted line is the variation of along slope coordinate (Y-axis coordinate); Fig. 3 (b) follows the tracks of the primitive curve that obtains to gauge point 1 among Fig. 2 (b) in the radiography graphic sequence, wherein solid line is the variation of lateral coordinates (X-axis coordinate), and dotted line is the variation of along slope coordinate (Y-axis coordinate);

Fig. 3 (c) decomposes the heart movement curve that obtains to primitive curve among Fig. 3 (a), and wherein solid line is the variation of lateral coordinates (X-axis coordinate), and dotted line is the variation of along slope coordinate (Y-axis coordinate);

Fig. 3 (d) decomposes the heart movement curve that obtains to primitive curve among Fig. 3 (b), and wherein solid line is the variation of lateral coordinates (X-axis coordinate), and dotted line is the variation of along slope coordinate (Y-axis coordinate);

Fig. 3 (e) decomposes the respiratory movement curve that obtains to primitive curve among Fig. 3 (a), and wherein solid line is the variation of lateral coordinates (X-axis coordinate), and dotted line is the variation of along slope coordinate (Y-axis coordinate);

Fig. 3 (f) decomposes the respiratory movement curve that obtains to primitive curve among Fig. 3 (b), and wherein solid line is the variation of lateral coordinates (X-axis coordinate), and dotted line is the variation of along slope coordinate (Y-axis coordinate);

Fig. 3 (g) is the result that respiratory movement curve among Fig. 3 (e) is carried out 6 curve fittings, and wherein solid line is the variation of lateral coordinates (X-axis coordinate), and dotted line is the variation of along slope coordinate (Y-axis coordinate);

Fig. 3 (h) is the result that respiratory movement curve among Fig. 3 (f) is carried out 6 curve fittings, and wherein solid line is the variation of lateral coordinates (X-axis coordinate), and dotted line is the variation of along slope coordinate (Y-axis coordinate);

Fig. 4 (a) is that wherein solid line is the variation of lateral coordinates (X-axis coordinate) at projection angle respiratory movement curve for (26.8 ° ,-27.2 °) following all gauge points, and dotted line is the variation of along slope coordinate (Y-axis coordinate);

Fig. 4 (b) is that wherein solid line is the variation of lateral coordinates (X-axis coordinate) at projection angle respiratory movement curve for (50.8 °, 30.2 °) following all gauge points, and dotted line is the variation of along slope coordinate (Y-axis coordinate);

Fig. 4 (c) is the result that the curve movement of being had a few among Fig. 4 (a) is carried out match, and wherein solid line is the variation of lateral coordinates (X-axis coordinate), and dotted line is the variation of along slope coordinate (Y-axis coordinate);

Fig. 4 (d) is the result that the curve movement of being had a few among Fig. 4 (b) is carried out match, and wherein solid line is the variation of lateral coordinates (X-axis coordinate), and dotted line is the variation of along slope coordinate (Y-axis coordinate);

Fig. 5 rebuilds the three-dimensional result that the back obtains to the radiography angle for the two-dimentional respiratory movement under (26.8 ° ,-27.2 °) and (50.8 °, 30.2 °).

The specific embodiment

The invention will be further described below in conjunction with accompanying drawing:

The present invention proposes a kind of human body respiration kinematic parameter extracting method based on frequency domain filtering, as shown in Figure 1, comprises the steps:

(1) obtain the one-arm X-ray radiography picture image sequence of arteria coronaria blood vessel, and definite heart movement cycle N1;

(2) choose the blood vessel structure characteristic point.

Guidance according to knowledge such as physiology and anatomy, the characteristic point of our labelling needs can concentrated expression to go out the movable information of blood vessel integral body, therefore, the different point of selected shape (being the architectural feature point) generally comprises the starting point and the terminating point of each vessel segment, and each flex point between vessel segment (point of curvature maximum).And in the radiography graphic sequence under two different projection angles,, corresponding mutually then have an identical name to all characteristic points numbering.As shown in Figure 2, be respectively at projected angle among a pair of radiography figure of (26.8 ° ,-27.2 °) and (50.8 °, 30.2 °), exist 15 different points of the shape after the numbering, their relation is according to the numeral name one to one.

(3) isolating cardiac and respirometric Fourier expansion method

Heart and human body respiration motion all are periodic motions, it is much smaller that but respirometric frequency is compared heart movement, in general, the frequency of heart proper motion is 60～100 times/minute, its cycle is 0.6-1.0s, the respirometric cycle is then much longer, is generally 3-6s, may be longer in the time of quietly.On the other hand, heart movement is more violent, and respirometric amplitude is less, i.e. the displacement of Chan Shenging is less, is a process relatively stably.

According to the motion of radiography figure cardiac and respirometric these characteristics, the present invention separates with respiratory movement heart in conjunction with the method for frequency domain filtering by fourier progression expanding method.

At first introduce the discrete Fourier series expansion below.

1) discrete Fourier series expansion

By the principle of fourier progression expanding method, continuous time, periodic signal can be expressed with fourier series, and is corresponding therewith, and the also available discrete Fourier progression of periodic sequence is represented.If periodic sequence x _p(n) cycle is N, so

x _p(n)=x _p(n+rN) (r is an arbitrary integer) (1)

x _p(n) the fourier series conversion that can disperse by formula (2) and (3)

$X_p\left(k\right)=\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}{x}_p\left(n\right)*e^{-j\left(\frac{2\mathrm{\π}}{N}\right)\mathrm{nk}}---\left(2\right)$

$x_p\left(n\right)=\frac{1}{N}\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}{X}_p\left(k\right)*e^{j\left(\frac{2\mathrm{\π}}{N}\right)\mathrm{nk}}---\left(3\right)$

Make description below for formula (3): in the formula

Be the fundamental component of periodic sequence, Be exactly the k order harmonic components, X _p(k) be the coefficient of each harmonic, also can think x _p(n) frequency domain response; It is independently that whole harmonic components have only N, because

$e^{j\left(\frac{2\mathrm{\π}}{N}\right)n(k+N)}=e^{j\left(\frac{2\mathrm{\π}}{N}\right)\mathrm{nk}},$

Therefore, progression get and item number be from k=0 to N-1, N independent harmonic component altogether.And formula (2) is just by x _p(n) coefficient of determination X _p(k) sum formula.

2) separation algorithm

Suppose that certain puts P (x on the arteria coronaria blood vessel in the radiography figure image sequence, the curve movement of x axial coordinate y) is x (n) (n is the frame number of angiogram frames), and the curve movement of y axial coordinate is y (n), makes s (n)=(x (n), y (n)), s (n) can be resolved into following formula:

s(n)＝c(n)+r(n)+t(n)

C (n)=(x wherein _c(n), y _c(n)) motion of the kinetic puncta vasculosa of expression heart, r (n)=(x _r(n), y _r(n)) represent the motion that respiratory movement causes, t (n)=(x _t(n), y _t(n)) represent some other motion (comprise people's the moving of health in angiographic procedure, and the moving etc. of radiography equipment).Be convenient expression, the back all uses puncta vasculosa that s (n) represents to extract along the x axle, the changes in coordinates curve of y axle, the puncta vasculosa that c (n) expression heart movement causes is along the x axle, the changes in coordinates curve of y axle, and the puncta vasculosa that r (n) expression respiratory movement causes is along the x axle, the changes in coordinates curve of y axle, therefore, be exactly respectively to the operation of x (n) and y (n) to the operation of s (n), be exactly respectively to x to the operation of c (n) _c(n) and y _c(n) operation is exactly respectively to x to the operation of r (n) _r(n) and y _r(n) operation.Because t (n) is immesurable, it is neglected do not consider here, then has

s(n)≈c(n)+r(n)

The cycle of supposing heart movement is N1, and respiration motion cycle is N2, and then c (n) and r (n) can be changed into according to formula (3):

$c\left(n\right)=\frac{1}{N1}\underset{k=0}{\overset{N1-1}{\mathrm{\Σ}}}C\left(k\right)*e^{j\left(\frac{2\mathrm{\π}}{N1}\right)\mathrm{nk}}$

$r\left(n\right)=\frac{1}{\mathrm{<figure-callout\; id="2"\; label="N"\; filenames="H2009102735285E00013.png"\; state="\{\{state\}\}">N</figure-callout>}2}\underset{k=0}{\overset{N2-1}{\mathrm{\&Sigma;}}}R\left(k\right)*e^{j\left(\frac{2\mathrm{\&pi;}}{N2}\right)\mathrm{nk}}$

Wherein,

Be the harmonic component of heart movement, Be respirometric harmonic component, C (k), R (k) represent the each harmonic coefficient of c (n) and r (n) respectively.C (n) is separated fully with r (n), then its harmonic component can not overlap that (the zero-frequency component can not considered, because though the frequency domain components of respirometric frequency and heart movement (DC component) on zero-frequency all exists, here also have no idea and they can be separated fully, but because the zero-frequency component variation is a constant in time domain, therefore it can not change the shape of periodic sequence curve, can not influence us to heart and respirometric analysis yet), promptly

$\&ForAll;k1\left(n1\right)=1,...,N1-1,k2\left(n2\right)=1,...,N2-1,$

(8)

$\&Exists;\left(\frac{2\mathrm{\&pi;}}{N1}\right)n1k1\&NotEqual;\left(\frac{2\mathrm{\&pi;}}{N2}\right)n2\mathrm{<figure-callout\; id="2"\; label="k"\; filenames="H2009102735285E00013.png"\; state="\{\{state\}\}">k</figure-callout>}2,$ (n1, n2 are the time domain independent variable, and k1, k2 are the frequency domain independent variable)

Satisfy formula (8), must make N1 and N2 relatively prime.To separate c (n) and r (n) simultaneously, also must satisfy the one-period of the sequence length of s (n), promptly greater than N1*N2 more than or equal to s (n).

Yet above-mentioned requirements is difficult to reach, even if N1 and N2 are relatively prime, second condition also can't satisfy.Because generally speaking, the time that contrast agent stops in human body can be not oversize, a radiography graphic sequence only can continue about 6s, can not surpass 10s, suppose that frame speed is 80ms/frame (actual capabilities are littler), then N1 ≈ 10, N2＞30, N1*N2＞300＞10/0.08, the therefore basic sequence that does not obtain one-period.

Though can't carry out ideal the separation with respiratory movement to heart,, can carry out proximate the separation by the method for fourier series conversion and frequency domain filtering according to the characteristics of their each autokinesis.

Separating step is based on following hypothesis:

(1) heart movement is only in frequency

On have component.This point through type (3) proves;

(2) respirometric frequency domain components is in frequency

The place is level and smooth.Because respiratory movement has flatness, and

Less for the probability of respiratory movement fundamental component, have reason fully to suppose like this.

Specific algorithm is as follows:

Step1: determine heart movement cycle N1 according to eye-observation radiography graphic sequence, generally speaking, N1 ≈ 0.8/f, wherein, N1 is illustrated in the radiography figure frame number in the cardiac cycle, and the time of a cardiac cycle of 0.8 expression, unit is second, f represents frame frequency, i.e. the quantity that each second, radiography figure was taken;

Step2: the blood vessel structure characteristic point of choosing in the step (1) is followed the tracks of in whole radiography graphic sequence, obtain tracking sequence s (n) a little;

Step3: in the tracking sequence s (n) of point, choose the sequence that length is ns

Ns is the integral multiple of N1.If original series length is n, then choosing sequence length is ns=n-n%N1, and wherein, % is for getting surplus symbol;

Step4: will Motion x (n) that is decomposed in the x direction and the motion y (n) on the y direction carry out discrete Fourier transform (DFT) to them respectively again, obtain each harmonic coefficient X (k) and Y (k), k=0 ..., ns-1;

Step5: make respirometric frequency domain response

$R_x\left(k\right)=\left\{\begin{array}{cc}X\left(k\right),& (k\&NotEqual;\frac{\mathrm{ns}}{N1}l,l=1,...,N1-1)\\ (X\left(k\right)+X(k+1))/2,& (k=\frac{\mathrm{ns}}{N1}l,l=1,...,N1-1)\end{array}\right.,$ To R _x(k) carry out inverse discrete Fourier transformer inverse-discrete and obtain respiratory movement r _x(n), in like manner can get r _y(n).

By analysis to the radiography graphic sequence, can determine that heart movement cycle N1 is 10 frames, Fig. 3 carries out respiratory movement and the isolating result of heart movement to the motion of gauge point among Fig. 21.

As can see from Figure 3, it is bigger, more in disorder that the original motion curve of gauge point is influenced by respiratory movement, and the heart movement curve display after separating good period, but then still there are some burrs in the respiratory movement curve, and this also is to separate inadequately reason completely.In order to remove heart movement component also residual in the respiratory movement, here the primary respiratory movement curve that extraction is obtained carries out curve fitting, remove the burr on the respiratory movement curve, shown in Fig. 3 (g) and 3 (h), the respiratory movement curve display after the match obvious periodic.

The respiratory movement curve of all gauge points is compared together, can find that the respiratory movement curve of being had a few in the same radiography face is more or less the same (as Fig. 4 (a) with (b)), so we become the respiratory movement curve (as Fig. 4 (c) with (d) shown in) of a curve as arteria coronaria blood vessel in this radiography face with the respiratory movement match of all gauge points.

(3) respirometric three-dimensional reconstruction

By respiratory movement extracting method recited above, what obtain at last is two-dimentional respiratory movement curve under certain radiography angle, as Fig. 4 (c) and Fig. 4 (d).But as can be seen from the figure, this two-dimensional representation method obviously can not be expressed respiration information clearly, so, in order more intuitively more effectively to describe respiratory movement, we can come the two dimensional motion that obtains under two different radiographies angles is rebuild according to the three-dimensional reconstruction principle of radiography system mid point, thereby obtain three-dimensional respiratory movement.Concrete steps are as follows:

Step1: determine two radiography figure image sequences under the different radiography angles, be designated as left radiography figure image sequence and right radiography figure image sequence, in the radiography figure image sequence of the left and right sides, determine respectively more corresponding arbitrary reference point to be designated as p respectively _l(x, y, t) and p _r(x, y, t), wherein (x is the coordinate of this arbitrary reference point in image sequence y), and t is the time, i.e. the frame number of radiography figure image sequence, and t is the integer between 0 to 70, chooses the moment t among the t respectively _l, t _rUnder coordinate (x _l, y _l), (x _r, y _r) as reference point p _l(x, y, t) and p _r(x, y, initialization value t), i.e. p _l(x _l, y _l, t _l) and p _r(x _r, y _r, t _r), t wherein _l, t _rCorresponding to the synchronization of cardiac cycle, generally select diastole latter stage;

Step2: extract the respiratory movement curve of left and right sides radiography figure image sequence according to above-mentioned respiratory movement parameter extracting method, be designated as curve respectively _l(x, y, t) and curve _r(x, y, t), then respectively on two curves corresponding the extreme point of air-breathing latter stage in the respiration motion cycle or EEP of selecting as breathing reference point, promptly on two curves, all select extreme point or all select the EEP extreme point in air-breathing latter stage, be designated as curve _l(x _Cl, y _Cl, t _Cl) and curve _r(x _Cr, y _Cr, t _Cr), t wherein _Cl, t _CrBe respectively the air-breathing moment in latter stage or the EEP moment corresponding in the radiography figure image sequence of the left and right sides, (x _Cl, y _Cl), (x _Cr, y _Cr) be respectively curve curve _l(x, y, t) and curve _r(x, y t) go up corresponding t constantly _Cl, t _CrCoordinate;

Step3: to the point after the initialization in the radiography figure image sequence of the left and right sides to p _l(x _l, y _l, t _l) and p _r(x _r, y _r, t _r) carry out respiration motion compensation, be about to t _l, t _rThe time respiratory movement of inscribing compensate to corresponding breathing respectively with reference to moment t _Cl, t _CrDown:

$\left\{\begin{array}{c}{p}_l^{\&prime;}(x_l^{\&prime;},y_l^{\&prime;},t_l)=p_l(x_l,y_l,t_l)-({\mathrm{curve}}_l(x_{\mathrm{tl}},y_{\mathrm{tl}},t_l)-{\mathrm{curve}}_l(x_{\mathrm{cl}},y_{\mathrm{cl}},t_{\mathrm{cl}}))\\ p_r^{\&prime;}(x_r^{\&prime;},y_r^{\&prime;},t_r)=p_r(x_r,y_r,t_r)-({\mathrm{curve}}_r(x_{\mathrm{tr}},y_{\mathrm{tr}},t_r)-{\mathrm{curve}}_r(x_{\mathrm{cr}},y_{\mathrm{cr}},t_{\mathrm{cr}}))\end{array}\right.$

In the formula, curve _l(x _Il, y _Tl, t _l), curve _r(x _Tr, y _Tr, t _r) represent t respectively _l, t _rMoment curve curve _l(x, y, t) and curve _r(x, y, the t) point on, p ' _l(x ' _l, y ' _l, t _l) and p ' _r(x ' _r, y ' _r, t _r) for compensating the back reference point.Then the reference point after two compensation is carried out three-dimensional reconstruction, obtain three-dimensional point P (x _Clr, y _Clr, z _Clr, t _Clr), t wherein _ClrExpression and t _Cl(or t _Cr) corresponding consistent air-breathing latter stage or EEP;

Step4: the reasonable assumption heart itself is immobilized, and what cause the motion of left coronary artery vascular tree only is Repiration.So, can be with the p ' that obtains among the Step3 _l(x ' _l, y ' _l, t _l) and p ' _r(x ' _r, y ' _r, t _r) under the influence that does not have heart movement, further compensate to respiration motion cycle other constantly:

$\left\{\begin{array}{c}{p}_l^{\&prime;}(x_{l+i}^{\&prime;},y_{l+i}^{\&prime;},t_{l+i})=p_l^{\&prime;}(x_l^{\&prime;},y_l^{\&prime;},t_l)+({\mathrm{curve}}_l(x_{l+i},y_{l+i},t_{l+i})-{\mathrm{curve}}_l(x_{\mathrm{cl}},y_{\mathrm{cl}},t_{\mathrm{cl}}))\\ p_r^{\&prime;}(x_{r+i}^{\&prime;},y_{r+i}^{\&prime;},t_{r+i})=p_r^{\&prime;}(x_r^{\&prime;},y_r^{\&prime;},t_r)+({\mathrm{curve}}_r(x_{r+i},y_{r+i},t_{r+i})-{\mathrm{curve}}_r(x_{\mathrm{cr}},y_{\mathrm{cr}},t_{\mathrm{cr}}))\end{array}\right.$

In like manner, to other points constantly of respiration motion cycle to p ' _l(x ' _L+i, y ' _L+i, t _L+i) and p ' _r(x ' _R+i, y ' _R+i, t _R+i) carry out three-dimensional reconstruction, obtain three-dimensional point P (x _Clr+i, y _Clr+i, z _Clr+i, t _Clr+i), wherein i represents with respect to breathing with reference to moment t _lAnd t _rThe frame number of preceding (getting negative integer) or back (getting positive integer);

Step5: the result that comprehensive Step3 and Step4 ask for can obtain three-dimensional respiratory movement.

By being (26.8 ° to the radiography angle,-27.2 °) (left side) and (50.8 °, 30.2 °) radiography graphic sequence under (right side) analyzes, choose the 22nd frame that is in diastole latter stage together and the 43rd frame respectively as the reference frame, choose blood vessel starting point on it again as the reference point, and be initialized as (35,62) and (303,94), analyze the respiratory movement curve that extracts in the radiography graphic sequence of the left and right sides then, determining to breathe reference respectively is the 38th frame and the 26th frame constantly.So, according to each parameter of above selection, can obtain three-dimensional respiratory movement, as shown in Figure 5.

As can be seen from the figure, three-dimensional respirometric overall trend is to come and go to do translational motion, but still there is the motion knuckle in the part, is not simple rectilinear motion.By analysis, the generation of this knuckle phenomenon should not separated from respiratory movement fully by heart movement influence and caused.But because the motion excursion that knuckle causes is little to entire effect, therefore, under current conditions, three-dimensional respiratory movement can reasonable assumption be the segmentation rectilinear motion.

Claims (2)

1. the human body respiration kinematic parameter extracting method based on frequency domain filtering comprises the steps:

(1) obtain the one-arm X-ray radiography picture image sequence of arteria coronaria blood vessel, and definite heart movement cycle N1;

(2) choose arteria coronaria blood vessel structure characteristic point;

(3) extract the respiratory movement curve, concrete mode is:

(3.1) in radiography figure image sequence described arteria coronaria blood vessel structure characteristic point is carried out from motion tracking, obtain feature point tracking sequence s (n), n is the feature point tracking sequence length; (3.2) in feature point tracking sequence s (n), choose one section as target sequence

Ns=n-n%N1, wherein, % is for getting surplus symbol;

(3.3) to target sequence

Carry out discrete Fourier transform (DFT), obtain target frequency domain response S (k), k=0 ..., ns-1;

(3.4) make respirometric frequency domain response

$R\left(k\right)=\left\{\begin{array}{cc}S\left(k\right),& (k\&NotEqual;\frac{\mathrm{ns}}{N1}l,l=1,...,N1-1)\\ (S\left(k\right)+S(k+1))/2,& (k=\frac{\mathrm{ns}}{N1}l,l=1,...,N1-1)\end{array}\right.,$

R (k) is carried out inverse discrete Fourier transformer inverse-discrete, just obtain respiratory movement curve r (ns), promptly can extract respiratory movement parameter.

2. one kind is carried out three-dimensional reconstruction to the respiratory movement curve described in the claim 1 and obtains three-dimensional respirometric method, and concrete steps comprise:

(I) determine two radiography figure image sequences under the different radiography angles, be designated as left radiography figure image sequence and right radiography figure image sequence, in the radiography figure image sequence of the left and right sides, determine respectively more corresponding arbitrary reference point to be designated as p respectively _l(x, y, t) and p _r(x, y, t), wherein (x y) is the coordinate of this reference point in image sequence, and t is the time, i.e. the frame number of radiography figure image sequence, and t is the integer between 0 to 70, chooses the moment t among the t respectively _l, t _rUnder coordinate (x _l, y _l), (x _r, y _r) as reference point p _l(x, y, t) and p _r(x, y, initialization value t), i.e. p _l(x _l, y _l, t _l) and p _r(x _r, y _r, t _r), t wherein _l, t _rIn the radiography graphic sequence of the left and right sides corresponding to the synchronization of cardiac cycle;

(II) extract the respiratory movement curve of left and right sides radiography figure image sequence according to the described respiratory movement parameter extracting method of claim 1, be designated as curve respectively _l(x, y, t) and curve _r(x, y, t), then respectively on two curves corresponding the extreme point of air-breathing latter stage in the respiration motion cycle or EEP of selecting as breathing reference point, promptly on two curves, all select extreme point or all select the EEP extreme point in air-breathing latter stage, be designated as curve _l(x _Cl, y _Cl, t _Cl) and curve _r(x _Cr, y _Cr, t _Cr), t wherein _Cl, t _CrBe respectively the air-breathing moment in latter stage or the EEP moment corresponding in the radiography figure image sequence of the left and right sides, (x _Cl, y _Cl), (x _Cr, y _Cr) be respectively curve curve _l(x, y, t) and curve _r(x, y t) go up corresponding t constantly _Cl, t _CrCoordinate;

(III) to the point after the initialization in the radiography figure image sequence of the left and right sides to p _l(x _l, y _l, t _l) and p _r(x _r, y _r, t _r) carry out respiration motion compensation, be about to t _l, t _rThe time respiratory movement of inscribing compensate to corresponding breathing respectively with reference to moment t _Cl, t _CrDown:

$\left\{\begin{array}{c}{p}_l^{\&prime;}(x_l^{\&prime;},y_l^{\&prime;},t_l)=p_l(x_l,y_l,t_l)-({\mathrm{curve}}_l(x_{\mathrm{tl}},y_{\mathrm{tl}},t_l)-{\mathrm{curve}}_l(x_{\mathrm{cl}},y_{\mathrm{cl}},t_{\mathrm{cl}}))\\ p_r^{\&prime;}(x_r^{\&prime;},y_r^{\&prime;},t_r)=p_r(x_r,y_r,t_r)-({\mathrm{curve}}_r(x_{\mathrm{tr}},y_{\mathrm{tr}},t_r)-{\mathrm{curve}}_r(x_{\mathrm{cr}},y_{\mathrm{cr}},t_{\mathrm{cr}}))\end{array}\right.$

In the formula, curve _l(x _Tl, y _Tl, t _l), curve _r(x _Tr, y _Tr, t _r) represent t respectively _l, t _rMoment curve curve _l(x, y, t) and curve _r(x, y, the t) point on, p ' _l(x ' _l, y ' _l, t _l) and p ' _r(x ' _r, y ' _r, t _r) for compensating the back reference point.Then the reference point after two compensation is carried out three-dimensional reconstruction, obtain three-dimensional point P (x _Clr, y _Clr, z _Clr, t _Clr), t wherein _ClrExpression and t _ClOr t _CrCorresponding consistent air-breathing latter stage or EEP;

(IV) with the p ' that obtains in the step (III) _l(x ' _l, y ' _l, t _l) and p ' _r(x ' _r, y ' _r, t _r) under the influence that does not have heart movement, further compensate to respiration motion cycle other constantly:

$\left\{\begin{array}{c}{p}_l^{\&prime;}(x_{l+i}^{\&prime;},y_{l+i}^{\&prime;},t_{l+i})=p_l^{\&prime;}(x_l^{\&prime;},y_l^{\&prime;},t_l)+({\mathrm{curve}}_l(x_{l+i},y_{l+i},t_{l+i})-{\mathrm{curve}}_l(x_{\mathrm{cl}},y_{\mathrm{cl}},t_{\mathrm{cl}}))\\ p_r^{\&prime;}(x_{r+i}^{\&prime;},y_{r+i}^{\&prime;},t_{r+i})=p_r^{\&prime;}(x_r^{\&prime;},y_r^{\&prime;},t_r)+({\mathrm{curve}}_r(x_{r+i},y_{r+i},t_{r+i})-{\mathrm{curve}}_r(x_{\mathrm{cr}},y_{\mathrm{cr}},t_{\mathrm{cr}}))\end{array}\right.$

In like manner, to other points constantly of respiration motion cycle to p ' _l(x ' _L+i, y ' _L+i, t _L+i) and p ' _r(x ' _R+i, y ' _R+i, t _R+i) carry out three-dimensional reconstruction, obtain three-dimensional point P (x _Clr+i, y _Clr+i, z _Clr+i, t _Clr+i), wherein i is illustrated in t constantly _lAnd t _rBefore or after frame number, when for frame number the preceding, i gets negative integer, for after frame number the time, i gets positive integer;

(V) result that asks for of comprehensive step (III) and step (IV) can obtain three-dimensional respiratory movement.

Images