U.S. patent application number 14/960461 was filed with the patent office on 2016-06-30 for method for iteratively extracting motion parameters from angiography images.
The applicant listed for this patent is Huazhong University of Science and Technology. Invention is credited to Lihua DENG, Zhenghua HUANG, Qiong SONG, Tianxu ZHANG, Yaozong ZHANG.
Application Number | 20160189394 14/960461 |
Document ID | / |
Family ID | 56164832 |
Filed Date | 2016-06-30 |
United States Patent
Application |
20160189394 |
Kind Code |
A1 |
ZHANG; Tianxu ; et
al. |
June 30, 2016 |
METHOD FOR ITERATIVELY EXTRACTING MOTION PARAMETERS FROM
ANGIOGRAPHY IMAGES
Abstract
A method for extracting motion parameters from angiography
images using a multi-parameter model. The method includes: 1)
extracting I vascular structural feature points automatically from
a medical image of an angiography image sequence, and auto-tracking
the feature points respectively in the angiography image sequence
to obtain a tracking sequence of each feature point; 2) performing
a discrete Fourier transformation on the tracking sequence of each
feature point to obtain a discrete Fourier transformation result;
initializing an iterative parameter, and obtaining amplitude range
and frequency range of each frequency point of the discrete Fourier
transformation result; 3) performing a Fourier transformation on a
tracking sequence of each frequency point in the amplitude range
and the frequency range thereof to obtain Fourier transformation
results; and 4) performing an inverse Fourier transformation on the
Fourier transformation results, and obtaining an estimated minimum
mean square error of each frequency point.
Inventors: |
ZHANG; Tianxu; (Wuhan,
CN) ; HUANG; Zhenghua; (Wuhan, CN) ; DENG;
Lihua; (Wuhan, CN) ; ZHANG; Yaozong; (Wuhan,
CN) ; SONG; Qiong; (Wuhan, CN) |
|
Applicant: |
Name |
City |
State |
Country |
Type |
Huazhong University of Science and Technology |
Wuhan |
|
CN |
|
|
Family ID: |
56164832 |
Appl. No.: |
14/960461 |
Filed: |
December 7, 2015 |
Related U.S. Patent Documents
|
|
|
|
|
|
Application
Number |
Filing Date |
Patent Number |
|
|
PCT/CN2015/072681 |
Feb 10, 2015 |
|
|
|
14960461 |
|
|
|
|
Current U.S.
Class: |
382/132 |
Current CPC
Class: |
G06T 2207/10016
20130101; G06T 7/246 20170101; A61B 6/486 20130101; A61B 6/5217
20130101; G06T 2207/10116 20130101; G06T 7/262 20170101; G06T
2207/30048 20130101; A61B 6/504 20130101; G06T 2207/20056
20130101 |
International
Class: |
G06T 7/20 20060101
G06T007/20; G06T 17/00 20060101 G06T017/00; A61B 6/00 20060101
A61B006/00; G06T 7/00 20060101 G06T007/00 |
Foreign Application Data
Date |
Code |
Application Number |
Dec 30, 2014 |
CN |
201410844139.4 |
Claims
1. A method for extracting motion parameters from angiography
images, the method comprising: (1) extracting I vascular structural
feature points from a medical image of an angiography image
sequence, and tracking the feature points respectively in the
angiography image sequence to obtain a tracking sequence
{s.sub.i(n), i=1, . . . , I} of each feature point, where n is
frame number of the medical image in the angiography image
sequence; (2) performing a discrete Fourier transformation on the
tracking sequence {s.sub.i(n), i=1, . . . , I} of each feature
point in (1) to obtain a discrete Fourier transformation result
S.sub.i(k); (3) initializing an iterative parameter j=0, and
obtaining an amplitude range and a frequency range of each
frequency point of the discrete Fourier transformation result
S.sub.i (k) in (2); (4) performing a Fourier transformation on a
tracking sequence of each frequency point in the amplitude range
and the frequency range thereof to obtain Fourier transformation
results; (5) performing an inverse Fourier transformation on the
Fourier transformation results in (4), and obtaining an estimated
minimum mean square error of each frequency point; (6) determining
whether the estimated minimum mean square error is greater than a
predetermined threshold, and proceeding to (7) if yes, otherwise
ending the process; (7) processing spectrums of each frequency
point by a multi-parameter iterative optimizing algori.sup.thm to
obtain (j+1).sup.th iterated time-domain signals; (8) processing a
residual signal by a translation model to obtain a (j+1).sup.th
iterated translation signal; (9) adding the (j+1).sup.th iterated
time-domain signals to the (j+1).sup.th iterated translation signal
to obtain an (j+1).sup.th iterated estimated mixed signal, and
calculating a (j+1).sup.th iterated minimum mean square error; and
(10) determining whether the (j+1).sup.th iterated minimum mean
square error is greater than the threshold in (6), and returning to
(7) if yes, otherwise ending the process.
2. The method of claim 1, wherein in (1), s.sub.i(n) is expressed
by the following equation:
s.sub.i(n)=L(n)+r.sub.i(n)+c.sub.i(n)+h.sub.i(n)+t.sub.i(n),i.di-elect
cons.[1, I], where L(n) represents translational movement,
r.sub.i(n) represents breathing movement of an i.sup.th feature
point, c.sub.i(n) represents cardiac movement of the i.sup.th
feature point, h.sub.i(n) represents tremor movement of the
i.sup.th feature point, and t.sub.i(n) represents other movements
of the i.sup.th feature point.
3. The method of claim 2, wherein in (2), .sup.Si (.sup.k) is
expressed by the following equation:
S.sub.i(k)=L(k)+R.sub.i(k)+C.sub.i(k)+H.sub.i(k), where k
represents a frequency point, and L(k), C(k), R(k) and H(k)
represent harmonic coefficients of L(n), c(n), r(n) and h(n)
correspondingly and respectively.
4. The method of claim 3, wherein in (5), the estimated minimum
mean square error {circumflex over (.epsilon.)}.sub.i.sup.j of the
frequency point is expressed by the following equation: ^ i j = min
( 1 N n ( s i ( n ) - s ^ i j ( n ) ) 2 ) , ##EQU00017## where
s.sub.i.sup.j(n)=L.sup.j(n)+r.sub.i.sup.j(n)+c.sub.i.sup.j(n)+h.sub.i.sup-
.j(n).
5. The method of claim 4, wherein (7) further comprises the
following sub-steps of: (7.1) calculating values L.sup.j(k.sub.ic),
R.sub.i.sup.j(k.sub.ic) and H.sub.i.sup.j(k.sub.icIC) of
L.sup.j(k), R.sub.i.sup.j(k) and H.sub.i.sup.j(k) near a frequency
point k.sub.ic in the frequency range respectively by the following
equation while keeping L.sup.j(k), R.sub.i.sup.j(k) and
H.sub.i.sup.j(k) constant: X p ( k ) = n = 0 N - 1 x p ( n ) - j (
2 .pi. N ) nk , ##EQU00018## calculating a (j+1).sup.th iterated
cardiac signal spectrum by an equation
C.sub.i.sup.j+1(k.sub.ic)=C.sub.i.sup.0(k.sub.ic)-R.sub.i.sup.j(k.sub.ic)-
-H.sub.i.sup.j(k.sub.ic), performing a discrete Fourier
transformation thereon in the frequency range, and obtaining an
optimized (j+1).sup.th cardiac time-domain signal
c.sub.i.sup.j+1(n) by the following equation: e ^ i j = min ( k =
.omega. - M .omega. + M ( S i ( k ) - S ^ i j ( k ) ) 2 ) ,
##EQU00019## where .omega. is a frequency point after the discrete
Fourier transformation, M represents a window size which is set to
3 normally, and
S.sub.i.sup.j(k)==L.sup.j(k)+R.sub.i.sup.j(k)+C.sub.i.sup.j+1(k)+H.s-
ub.i.sup.j(k); (7.2) calculating values L.sup.j(k.sub.ih),
R.sub.i.sup.j(k.sub.ih) and C.sub.i.sup.j+1(k.sub.ih) of
L.sup.j(k), R.sub.i.sup.j(k) and C.sub.i.sup.j+1(k) near a
frequency point k.sub.ih in the frequency range respectively by the
following equation while keeping L.sup.j(k), R.sub.i.sup.j(k) and
C.sub.i.sup.j+1(k) constant: X p ( k ) = n = 0 N - 1 x p ( n ) - j
( 2 .pi. N ) nk , ##EQU00020## calculating a (j+1).sup.th iterated
high-frequency signal spectrum by an equation
H.sub.i.sup.j+1(k.sub.ih)=H.sub.i.sup.0(k.sub.ih)-L.sup.j(k.sub.ih)-R.sub-
.i.sup.j(k.sub.ih)-C.sub.i.sup.j+1(k.sub.ih), performing a discrete
Fourier transformation thereon in the frequency range, and
obtaining an optimized (j+1).sup.th high-frequency time-domain
signal h.sub.i.sup.j-1(n) by the following equation: e ^ i j = min
( k = .omega. - M .omega. + M ( S i ( k ) - S ^ i j ( k ) ) 2 ) ,
##EQU00021## where
S.sub.i.sup.j(k)==L.sup.j(k)+R.sub.i.sup.j(k)+C.sub.i.sup.j+1(k)+H.sub.i.-
sup.j-1(k); and (7.3) calculating values L.sup.j(k.sub.ir),
H.sub.i.sup.j+1(k.sub.ir) of L.sup.j(k), H.sub.i.sup.j+1(k) and
C.sub.i.sup.j+1(k) near a frequency point k.sub.ir in the frequency
range respectively by the following equation while keeping
L.sup.j(k), H.sub.i.sup.j+1(k) and C.sub.i.sup.j+1(k) constant: X p
( k ) = n = 0 N - 1 x p ( n ) - j ( 2 .pi. N ) nk , ##EQU00022##
calculating a (j+1).sup.th iterated breathing signal spectrum by an
equation
R.sub.i.sup.j+1(k.sub.ir)=R.sub.i.sup.0(k.sub.ir)-L.sup.j(k.sub.ir)-C.sub-
.i.sup.j-1(k.sub.ir)-H.sub.i.sup.j+1(k.sub.ir), performing a
discrete Fourier transformation thereon in the frequency range, and
obtaining an optimized (j+1).sup.th breathing time-domain signal
r.sub.i.sup.j+1(n) by the following equation: e ^ i j = min ( k =
.omega. - M .omega. + M ( S i ( k ) - S ^ i j ( k ) ) 2 ) ,
##EQU00023## where
S.sub.i.sup.j(k)==L.sup.j(k)+R.sub.i.sup.j+1(k)+C.sub.i.sup.j+1(k)+H.sub.-
i.sup.j+1(k).
6. The method of claim 5, wherein in (9), the (j+1).sup.th iterated
minimum mean square error {circumflex over
(.epsilon.)}.sub.i.sup.j+1 is expressed by the following equation:
^ i j + 1 = min ( 1 N n ( s i ( n ) - s ^ i j + 1 ( n ) ) 2 ) .
##EQU00024##
7. A system for extracting motion parameters from angiography
images, the system comprising: a) a first module, operable for
extracting I vascular structural feature points from a medical
image of an angiography image sequence, and tracking the feature
points respectively in the angiography image sequence to obtain a
tracking sequence {s.sub.i(n), i=1, . . . , I} of each feature
point, where n is frame number of the medical image in the
angiography image sequence; b) a second module, operable for
performing a discrete Fourier transformation on the tracking
sequence {s.sub.i(n), i=1, . . . , I} of each feature point derived
by the first module to obtain a discrete Fourier transformation
result S.sub.i(k); c) a third module, operable for initializing an
iterative parameter j=0, and obtaining amplitude range and
frequency range of each frequency point of the discrete Fourier
transformation result S.sub.i(k) derived by the second module; d) a
fourth module, operable for performing a Fourier transformation on
a tracking sequence of each frequency point in the amplitude range
and the frequency range thereof to obtain Fourier transformation
results; e) a fifth module, operable for performing an inverse
Fourier transformation on the Fourier transformation results
derived by the fourth module, and obtaining an estimated minimum
mean square error of each frequency point; f) a sixth module,
operable for determining whether the estimated minimum mean square
error is greater than a predetermined threshold, and proceeding to
a seventh module if yes, otherwise ending the process; g) a seventh
module, operable for processing spectrums of each frequency point
by a multi-parameter iterative optimizing algori.sup.thm to obtain
(j+1).sup.th iterated time-domain signals; h) an eighth module,
operable for processing a residual signal by a translation model to
obtain a (j+1).sup.th iterated translation signal; i) a ninth
module, operable for adding the (j+1).sup.th iterated time-domain
signals to the (j+1).sup.th iterated translation signal to obtain
an (j+1).sup.th iterated estimated mixed signal, and calculating a
(j+1).sup.th iterated minimum mean square error; and j) a tenth
module, operable for determining whether the (j+1).sup.th iterated
minimum mean square error is greater than the threshold in the
sixth module, and returning to the seventh module if yes, otherwise
ending the process.
Description
CROSS-REFERENCE TO RELATED APPLICATIONS
[0001] This application is a continuation-in-part of International
Patent Application No. PCT/CN2015/072681 with an international
filing date of Feb. 10, 2015, designating the United States, now
pending, and further claims priority benefits to Chinese Patent
Application No. 201410844139.4 filed Dec. 30, 2014. The contents of
all of the aforementioned applications, including any intervening
amendments thereto, are incorporated herein by reference. Inquiries
from the public to applicants or assignees concerning this document
or the related applications should be directed to: Matthias Scholl
P. C., Attn.: Dr. Matthias Scholl Esq., 245 First Street, 18th
Floor, Cambridge, Mass. 02142.
BACKGROUND OF THE INVENTION
[0002] 1. Field of the Invention
[0003] The invention relates to a method for iteratively extracting
motion parameters from angiography images using a multi-parameter
model.
[0004] 2. Description of the Related Art
[0005] Typically, coronary angiography images record not only
projection of cardiac movement on a 2D surface, but also 2D
movement of coronary arteries on an angiograph surface. Thus, to
reconstruct a 3D cardiovascular tree from dynamic 2D angiograms,
the cardiac movement and breathing movement must be separated.
[0006] A typical method for separation of the cardiac movement and
breathing movement includes presetting labels in human organs such
as heart and diaphragm, and tracking the labels to obtain a motion
curve as approximate breathing movement. The method requires many
trials and, thus, is inapplicable for clinical applications.
[0007] Another method is to obtain parameters of the breathing
movement and cardiac movement by double-armed X-ray angiography.
The method is relatively complex, requires harmful contrast agents,
and the motion parameters obtained thereby are not accurate
enough.
SUMMARY OF THE INVENTION
[0008] In view of the above-mentioned problems, it is an objective
of the invention to provide a method for iteratively extracting
motion parameters from angiography images using a multi-parameter
model, so as to solve technical problems of failing to extract
translational movement automatically and separate breathing
movement and cardiac movement accurately in prior art.
[0009] The invention provides a method for iteratively extracting
parameters of multiple movements (including translational movement,
breathing movement and cardiac movement) from an X-ray angiography
image sequence, selecting a few feature points on coronary
arteries, obtaining motion curves thereof by tracking them in an
image sequence, optimizing components of the movements under the
condition of a minimum mean square error between a reconstructed
signal and an original signal at frequency points in frequency
domain through a discrete Fourier transformation and an inverse
discrete Fourier transformation and using a multi-parameter model,
and obtaining optimized 2D motion curves for cardiac movement,
tremor movement, breathing movement and translational movement,
which is applicable for clinical applications.
[0010] To achieve the above objective, according to one embodiment
of the present invention, there is provided a method for
iteratively extracting motion parameters from angiography images
using a multi-parameter model, the method comprising: [0011] (1)
extracting/vascular structural feature points automatically from a
medical image of an angiography image sequence, and auto-tracking
the feature points respectively in the angiography image sequence
to obtain a tracking sequence {s.sub.i(n), i=1, . . . , I} of each
feature point, where n is frame number of the medical image in the
angiography image sequence; [0012] (2) performing a discrete
Fourier transformation on the tracking sequence {s.sub.i(n), i=1, .
. . , I} of each feature point in (1) to obtain a discrete Fourier
transformation result S.sub.i(k); [0013] (3) initializing an
iterative parameter j=0, and obtaining an amplitude range and a
frequency range of each frequency point of the discrete Fourier
transformation result S.sub.i(k) in (2); [0014] (4) performing a
Fourier transformation on a tracking sequence of each frequency
point in the amplitude range and the frequency range thereof to
obtain Fourier transformation results; [0015] (5) performing an
inverse Fourier transformation on the Fourier transformation
results in (4), and obtaining an estimated minimum mean square
error of each frequency point; [0016] (6) determining whether the
estimated minimum mean square error is greater than a predetermined
threshold, and proceeding to (7) if yes, otherwise ending the
process; [0017] (7) processing spectrums of each frequency point by
a multi-parameter iterative optimizing algori.sup.thm to obtain
(j+1).sup.th iterated time-domain signals; [0018] (8) processing a
residual signal by a translation model to obtain a (j+1).sup.th
iterated translation signal; [0019] (9) adding the (j+1).sup.th
iterated time-domain signals to the (j+1).sup.th iterated
translation signal to obtain an (j+1).sup.th iterated estimated
mixed signal, and calculating a (j+1).sup.th iterated minimum mean
square error; and [0020] (10) determining whether the (j+1).sup.th
iterated minimum mean square error is greater than the threshold in
(6), and returning to (7) if yes, otherwise ending the process.
[0021] In a class of this embodiment, in (1), s.sub.i(n) is
expressed by the following equation:
s.sub.i(n)=L(n)+r.sub.i(n)+c.sub.i(n)+h.sub.i(n)+t.sub.i(n),i.di-elect
cons.[1, I],
[0022] where L(n) represents translational movement, r.sub.i(n)
represents breathing movement of an i.sup.th feature point,
c.sub.i(n) represents cardiac movement of the i.sup.th feature
point, h.sub.i(n) represents tremor movement of the i.sup.th
feature point, and t.sub.i(h) represents other movements of the
i.sup.th feature point.
[0023] In a class of this embodiment, in (2), S.sub.i(k) is
expressed by the following equation:
S.sub.i(k)=L(k)+R.sub.i(k)+C.sub.i(k)+H.sub.i(k),
[0024] where k represents a frequency point, and L(k), C(k), R(k)
and H(k) represent harmonic coefficients of L(n), c(n), r(n) and
h(n) correspondingly and respectively.
[0025] In a class of this embodiment, in (5), the estimated minimum
mean square error {circumflex over (.epsilon.)}.sub.i.sup.j of the
frequency point is expressed by the following equation:
^ i j = min ( 1 N n ( s i ( n ) - s ^ i j ( n ) ) 2 ) ,
##EQU00001##
[0026] where
s.sub.i.sup.j(n)=L.sup.j(n)+r.sub.i.sup.j(n)+c.sub.i.sup.j(n)+h.sub.i.sup-
.j(n).
[0027] In a class of this embodiment, (7) further comprises:
[0028] (7.1) calculating values L.sup.j(k.sub.ic),
R.sub.i.sup.j(k.sub.ic) and H.sub.i.sup.j(k.sub.ic) of L.sup.j(k),
R.sub.i.sup.j(k) and H.sub.i.sup.j(k) near a frequency point
k.sub.ic in the frequency range respectively by the following
equation while keeping L.sup.j(k), R.sub.i.sup.j(k) and
H.sub.i.sup.j(k) constant:
X p ( k ) = n = 0 N - 1 x p ( n ) - j ( 2 .pi. N ) nk ,
##EQU00002##
[0029] calculating a (j+1).sup.th iterated cardiac signal spectrum
by an equation
C.sub.i.sup.j-1(k.sub.ic)=C.sub.i.sup.0(k.sub.ic)-R.sub.i.sup.j(-
k.sub.ic)-H.sub.i.sup.j(k.sub.ic), performing a discrete Fourier
transformation thereon in the frequency range, and obtaining an
optimized (j+1).sup.th cardiac time-domain signal
c.sub.i.sup.j+1(n) by the following equation:
e ^ i j = min ( k = .omega. - M .omega. + M ( S i ( k ) - S ^ i j (
k ) ) 2 ) , ##EQU00003##
[0030] where .omega. is a frequency point after the discrete
Fourier transformation, M represents a window size which is set to
3 normally, and
S.sub.i.sup.j(k)==L.sup.j(k)+R.sub.i.sup.j(k)+C.sub.i.sup.j+1(k)+H.su-
b.i.sup.j(k);
[0031] (7.2) calculating values L.sup.i(k.sub.ih),
R.sub.i.sup.j(k.sub.ih) and C.sub.i.sup.j+1(k.sub.ih) of
L.sup.j(k), R.sub.i.sup.j(k) and C.sub.i.sup.j+1(k) near a
frequency point k.sub.ih in the frequency range respectively by the
following equation while keeping L.sup.j(k), R.sub.i.sup.j(k) and
C.sub.i.sup.j+1(k) constant:
X p ( k ) = n = 0 N - 1 x p ( n ) - j ( 2 .pi. N ) nk ,
##EQU00004##
[0032] calculating a (j+1).sup.th iterated high-frequency signal
spectrum by an equation
H.sub.i.sup.j+1(k.sub.ih)=H.sub.i.sup.0(k.sub.ih)-L.sup.j(k.sub.ih)-R.sub-
.i.sup.j(k.sub.ih)-C.sub.i.sup.j+1(k.sub.ih), performing a discrete
Fourier transformation thereon in the frequency range, and
obtaining an optimized (j+1).sup.th high-frequency time-domain
signal h.sup.j+1 (n) by the following equation:
e ^ i j = min ( k = .omega. - M .omega. + M ( S i ( k ) - S ^ i j (
k ) ) 2 ) , ##EQU00005##
[0033] where
S.sub.i.sup.j(k)==L.sup.j(k)+R.sub.i.sup.j(k)+C.sub.i.sup.j+1(k)+H.sub.i.-
sup.j(k); and
[0034] (7.3) calculating values L.sup.j(k.sub.ir),
H.sub.i.sup.j+1(k.sub.ir) and C.sub.i.sup.j+1(k.sub.ir) of
L.sup.j(k), H.sub.i.sup.j+1(k) and C.sub.i.sup.j+1(k) near a
frequency point k.sub.ir in the frequency range respectively by the
following equation while keeping L.sup.j(k), H.sub.i.sup.j+1(k) and
C.sub.i.sup.j+1(k) constant:
X p ( k ) = n = 0 N - 1 x p ( n ) - j ( 2 .pi. N ) nk ,
##EQU00006##
[0035] calculating a (j+1).sup.th iterated breathing signal
spectrum by an equation
R.sub.i.sup.j+1(k.sub.ir)=R.sub.i.sup.0(k.sub.ir)-L.sup.j(k.sub.-
ir)-C.sub.i.sup.j+1(k.sub.ir)-H.sub.i.sup.j+1(k.sub.ir), performing
a discrete Fourier transformation thereon in the frequency range,
and obtaining an optimized (j+1).sup.th breathing time-domain
signal r.sub.i.sup.j+1(n) by the following equation:
e ^ i j = min ( k = .omega. - M .omega. + M ( S i ( k ) - S ^ i j (
k ) ) 2 ) , ##EQU00007##
[0036] where
S.sub.i.sup.j(k)==L.sup.j(k)+R.sub.i.sup.j+1(k)+C.sub.i.sup.j+1(k)+H.sub.-
i.sup.j+1(k).
[0037] In a class of this embodiment, in (9), the (j+1).sup.th
iterated minimum mean square error {circumflex over
(.epsilon.)}.sub.i.sup.j+1 is expressed by the following
equation:
^ i j + 1 = min ( 1 N n ( s i ( n ) - s ^ i j + 1 ( n ) ) 2 ) .
##EQU00008##
[0038] According to another embodiment of the present invention,
there is provided a system for iteratively extracting motion
parameters from angiography images using a multi-parameter model,
the system comprising: [0039] a first module, operable for
extracting I vascular structural feature points automatically from
a medical image of an angiography image sequence, and auto-tracking
the feature points respectively in the angiography image sequence
to obtain a tracking sequence {s.sub.i(n), i=1, . . . , I} of each
feature point, where n is frame number of the medical image in the
angiography image sequence; [0040] a second module, operable for
performing a discrete Fourier transformation on the tracking
sequence {s.sub.i(n), i=1, . . . , I} of each feature point derived
by the first module to obtain a discrete Fourier transformation
result S.sub.i(k); [0041] a third module, operable for initializing
an iterative parameter j=0, and obtaining amplitude range and
frequency range of each frequency point of the discrete Fourier
transformation result S.sub.i(k) derived by the second module;
[0042] a fourth module, operable for performing a Fourier
transformation on a tracking sequence of each frequency point in
the amplitude range and the frequency range thereof to obtain
Fourier transformation results; [0043] a fifth module, operable for
performing an inverse Fourier transformation on the Fourier
transformation results derived by the fourth module, and obtaining
an estimated minimum mean square error of each frequency point;
[0044] a sixth module, operable for determining whether the
estimated minimum mean square error is greater than a predetermined
threshold, and proceeding to a seventh module if yes, otherwise
ending the process; [0045] a seventh module, operable for
processing spectrums of each frequency point by a multi-parameter
iterative optimizing algori.sup.thm to obtain (j+1).sup.th iterated
time-domain signals; [0046] an eighth module, operable for
processing a residual signal by a translation model to obtain a
(j+1).sup.th iterated translation signal; [0047] a ninth module,
operable for adding the (j+1).sup.th iterated time-domain signals
to the (j+1).sup.th iterated translation signal to obtain an
(j+1).sup.th iterated estimated mixed signal, and calculating a
(j+1).sup.th iterated minimum mean square error; and [0048] a tenth
module, operable for determining whether the (j+1).sup.th iterated
minimum mean square error is greater than the threshold in the
sixth module, and returning to the seventh module if yes, otherwise
ending the process.
[0049] Advantages of the method for iteratively extracting motion
parameters from angiography images using a multi-parameter model
according to embodiments of the invention are summarized as
follows:
[0050] 1. with (2), the invention obtains structural feature points
automatically for motion curves, avoiding setting labels on surface
of the heart;
[0051] 2. with the optimizing method by multiple iterations in (7),
the invention solves problems of failing to extract translational
movement automatically and separate other periodic movements (such
as breathing movement and cardiac movement); and
[0052] 3. with the optimizing method by multiple iterations in (7),
the invention solves problems of low robustness with inaccurate
extracted motion parameters.
BRIEF DESCRIPTION OF THE DRAWINGS
[0053] FIG. 1 is a flow chart of a method for iteratively
extracting motion parameters from angiography images using a
multi-parameter model of the invention;
[0054] FIG. 2 is an original angiography image obtained from an
angle;
[0055] FIG. 3 illustrates feature points selected from the vascular
skeleton in FIG. 2;
[0056] FIGS. 4A and 4B are original motion curves of feature points
A1-A4 in FIG. 2 in an X direction and a Y direction;
[0057] FIGS. 5A and 5B are spectrum components of the fourth
feature point in FIG. 2 in the X direction and the Y direction;
[0058] FIGS. 6A and 6B are cardiac motion curves of feature points
in FIG. 2 in the X direction and the Y direction;
[0059] FIGS. 7A and 7B are breathing motion curves of feature
points in FIG. 2 in the X direction and the Y direction;
[0060] FIGS. 8A and 8B are high-frequency motion curves of feature
points in FIG. 2 in the X direction and the Y direction;
[0061] FIGS. 9A and 9B illustrate translational motion curves of
feature points in FIG. 2 in the X direction and the Y direction
compared with those obtained by tracking skeleton points
manually;
[0062] FIGS. 10A and 10B are residual curves of feature points in
FIG. 2 in the X direction and the Y direction;
[0063] FIG. 11 is an original angiography image obtained from
another angle;
[0064] FIG. 12 illustrates feature points selected from the
vascular skeleton in FIG. 11;
[0065] FIGS. 13A and 13B are original motion curves of feature
points A1-A5 in FIG. 11 in an X direction and a Y direction;
[0066] FIGS. 14A and 14B are spectrum components of the fourth
feature point in FIG. 11 in the X direction and the Y
direction;
[0067] FIGS. 15A and 15B are cardiac motion curves of feature
points in FIG. 11 in the X direction and the Y direction;
[0068] FIGS. 16A and 16B are cardiac motion curves of feature
points in FIG. 11 in the X direction and the Y direction;
[0069] FIGS. 17A and 17B are breathing motion curves of feature
points in FIG. 11 in the X direction and the Y direction; and
[0070] FIGS. 18A and 18B are residual curves of feature points in
FIG. 11 in the X direction and the Y direction.
DETAILED DESCRIPTION OF THE EMBODIMENTS
[0071] For clear understanding of the objectives, features and
advantages of the invention, detailed description of the invention
will be given below in conjunction with accompanying drawings and
specific embodiments. It should be noted that the embodiments are
only meant to explain the invention, and not to limit the scope of
the invention.
[0072] Technical terms of the invention are illustrated as
follows.
[0073] Minimum mean square error: a minimum value of a mean-squared
superimposition of differences between estimated signals and
original signals. Multi-parameter model of the invention obtains
optimized extracted signals by a method of minimum mean square
error performing iterations in time domain and frequency domain to
achieve a minimum residual signal in terms of the property of a
discrete Fourier transformation.
[0074] To solve the problems in prior art, the invention provides a
multi-parameter model of structural feature points, selecting a few
vascular structural feature points from angiography images,
obtaining motion curves thereof by tracking them in an image
sequence, optimizing components of the movements under the
condition of a minimum mean square error between a reconstructed
signal and an original signal at frequency points in frequency
domain through a discrete Fourier transformation and an inverse
discrete Fourier transformation combined with a translation model,
and obtaining optimized 2D motion curves for cardiac movement,
tremor movement, breathing movement and translational movement. The
method is flexible and can be widely used in almost all angiography
images (with a clear vascular distribution and covering two or more
heart beat periods, which is easy to be satisfied) including
double-armed X-ray angiography images.
[0075] As in FIG. 1, a method for iteratively extracting motion
parameters from angiography images using a multi-parameter model
comprises steps of:
[0076] (1) extracting I vascular structural feature points
automatically from a medical image of an angiography image sequence
(where I is a natural number), and auto-tracking the feature points
respectively in the angiography image sequence to obtain a tracking
{s.sub.i(n), i=1, . . . , I} sequence of each feature point, where
n is frame number of the medical image in the angiography image
sequence;
[0077] s.sub.i(n) is expressed by the following equation:
s.sub.i(n)=L(n)+r.sub.i(n)+c.sub.i(n)+h.sub.i(n)+t.sub.i(n),i.di-elect
cons.[1, I],
[0078] where L(n) represents translational movement, r.sub.i(n)
represents breathing movement of an i.sup.th feature point,
c.sub.i(n) represents cardiac movement of the i.sup.th feature
point, h.sub.i(n) represents tremor movement of the i.sup.th
feature point, and t.sub.i(n) represents other movements of the
i.sup.th feature point.
[0079] (2) performing a discrete Fourier transformation on the
tracking sequence {s.sub.i(n), i=1, . . . , I} of each feature
point in (1) to obtain a discrete Fourier transformation result
S.sub.i(k);
[0080] S.sub.i(k) is expressed by the following equation:
S.sub.i(k)=L(k)+R.sub.i(k)+C.sub.i(k)+H.sub.i(k),
[0081] where k represents a frequency point, and L(k), C(k), R(k)
and H(k) represent harmonic coefficients of L(n), c(n) , r(n) and
h(n) correspondingly and respectively.
[0082] (3) initializing an iterative parameter j=0, and obtaining
an amplitude range and a frequency range of each frequency point of
the discrete Fourier transformation result S.sub.i(k) in (2);
[0083] (4) performing a Fourier transformation on a tracking
sequence of each frequency point in the amplitude range and the
frequency range thereof to obtain Fourier transformation results
L.sup.j(k), R.sub.i.sup.j(k), C.sub.i.sup.j(k) and
H.sub.i.sup.j(k);
[0084] (5) performing an inverse Fourier transformation on the
Fourier transformation results L.sup.j(k), R.sub.i.sup.j(k),
C.sub.i.sup.j(k) and H.sub.i.sup.j(k) in (4), and obtaining an
estimated minimum mean square error {circumflex over
(.epsilon.)}.sub.i.sup.j of each frequency point:
^ i j = min ( 1 N n ( s i ( n ) - s ^ i j ( n ) ) 2 ) ,
##EQU00009##
[0085] where
s.sub.i.sup.j(n)=L.sup.j(n)+r.sub.i.sup.j(n)+c.sub.i.sup.j(n)+h.sub.i.sup-
.j(n);
[0086] (6) determining whether the estimated minimum mean square
error {circumflex over (.epsilon.)}.sub.i.sup.j is greater than a
predetermined threshold (a positive integer not greater than 2),
and proceeding to (7) if yes, otherwise ending the process;
[0087] (7) processing spectrums L.sup.j(k), R.sub.i.sup.j(k) and
H.sub.i.sup.j(k) of each frequency point by a multi-parameter
iterative optimizing algori.sup.thm to obtain (j+1).sup.th iterated
time-domain signals c.sub.i.sup.j+1(n), h.sub.i.sup.j-1(n),
r.sub.i.sup.j+1(n), etc.;
[0088] Step (7) further comprises the following sub-steps of:
[0089] (7.1) calculating values L.sup.j(k.sub.ic),
R.sub.i.sup.j(k.sub.ic) and H.sub.i.sup.j(k.sub.ic) of L.sup.j(k),
R.sub.i.sup.j(k) and H.sub.i.sup.j(k) near a frequency point
k.sub.ic in the frequency range respectively by the following
equation while keeping L.sup.j(k), R.sub.i.sup.j(k) and
H.sub.i.sup.j(k) constant:
X p ( k ) = n = 0 N - 1 x p ( n ) - j ( 2 .pi. N ) nk ,
##EQU00010##
[0090] calculating a (j+1).sup.th iterated cardiac signal spectrum
by an equation
C.sub.i.sup.j-1(k.sub.ic)=C.sub.i.sup.0(k.sub.ic)-L.sup.j(k.sub.-
ic)-R.sub.i.sup.j(k.sub.ic)-H.sub.i.sup.j(k.sub.ic), performing a
discrete Fourier transformation thereon in the frequency range, and
obtaining an optimized (j+1).sup.th cardiac time-domain signal
c.sub.i.sup.j+1(n) by the following equation:
e ^ i j = min ( k = .omega. - M .omega. + M ( S i ( k ) - S ^ i j (
k ) ) 2 ) , ##EQU00011##
[0091] where .omega. is a frequency point after the discrete
Fourier transformation, M represents a window size which is set to
3 normally, and
S.sub.i.sup.j(k)==L.sup.j(k)+R.sub.i.sup.j+1(k)+C.sub.i.sup.j+1(k)+H.-
sub.i.sup.j(k);
[0092] (7.2) calculating values L.sup.j(k.sub.ih),
R.sub.i.sup.j(k.sub.ih) and C.sub.i.sup.j+1(k.sub.ih) of
L.sup.j(k), R.sub.i.sup.j(k) and C.sub.i.sup.j+1(k) near a
frequency point k.sub.ih in the frequency range respectively by the
following equation while keeping L.sup.j(k), R.sub.i.sup.j(k) and
C.sub.i.sup.j+1(k) constant:
X p ( k ) = n = 0 N - 1 x p ( n ) - j ( 2 .pi. N ) nk ,
##EQU00012##
[0093] calculating a (j+1).sup.th iterated high-frequency signal
spectrum by an equation
H.sub.i.sup.j+1(k.sub.ih)=H.sub.i.sup.o(k.sub.ih)-L.sup.j(k.sub.ih)-R.sub-
.i.sup.j(k.sub.ih)-C.sub.i.sup.j+1(k.sub.ih), performing a discrete
Fourier transformation thereon in the frequency range, and
obtaining an optimized (j+1).sup.th high-frequency time-domain
signal h.sub.i.sup.j+1(n) by the following equation:
e ^ i j = min ( k = .omega. - M .omega. + M ( S i ( k ) - S ^ i j (
k ) ) 2 ) , ##EQU00013##
[0094] where
S.sub.i.sup.j(k)==L.sup.j(k)+R.sub.i.sup.j(k)+C.sub.i.sup.j+1(k)+H.sub.i.-
sup.j-1(k); and
[0095] (7.3) calculating values L.sup.j(k.sub.ir),
H.sub.i.sup.j+1(k.sub.ir) and C.sub.i.sup.j+1(k.sub.ir) of
L.sup.j(k), H.sub.i.sup.j+1(k) and C.sub.i.sup.j+1(k) near a
frequency point k.sub.ir in the frequency range respectively by the
following equation while keeping L.sup.j(k), H.sub.i.sup.j+1(k) and
C.sub.i.sup.j+1(k) constant:
X p ( k ) = n = 0 N - 1 x p ( n ) - j ( 2 .pi. N ) nk ,
##EQU00014##
[0096] calculating a (j+1).sup.th iterated breathing signal
spectrum by an equation
R.sub.i.sup.j+1(k.sub.ir)=R.sub.i.sup.0(k.sub.ir)-I.sub.i.sup.j(-
k.sub.ir)-C.sub.i.sup.j+1(k.sub.ir)-H.sub.i.sup.j+1(k.sub.ir),
performing a discrete Fourier transformation thereon in the
frequency range, and obtaining an optimized (j+1)t.sup.h breathing
time-domain signal r.sub.i.sup.j-1(n) by the following
equation:
e ^ i j = min ( k = .omega. - M .omega. + M ( S i ( k ) - S ^ i j (
k ) ) 2 ) , ##EQU00015##
[0097] where
S.sub.i.sup.j(k)==L.sup.j(k)+R.sub.i.sup.j+1(k)+C.sub.i.sup.j+1(k)+H.sub.-
i.sup.j+1(k).
[0098] (8) processing a residual signal
v.sub.i.sup.j+1(n)=s.sub.i(n)-c.sub.i.sup.j+1(n)-h.sub.i.sup.j+1(n)-r.sub-
.i.sup.j+1(n) by a translation model to obtain a (j+1).sup.th
iterated translation signal L.sup.j+1(n);
[0099] (9) adding the (j+1).sup.th iterated time-domain signals
c.sub.i.sup.j+1(n), h.sub.i.sup.j+1(n) and r.sub.i.sup.j+1(n) to
the (j+1).sup.th iterated translation signal (n) to obtain an (j+1)
iterated estimated mixed signal s.sub.i.sup.j+1(n), and calculating
a (j+1).sup.th iterated minimum mean square error {circumflex over
(.epsilon.)}.sub.i.sup.j+1:
^ i j + 1 = min ( 1 N n ( s i ( n ) - s ^ i j + 1 ( n ) ) 2 ) ;
##EQU00016## and ##EQU00016.2##
[0100] (10) determining whether the (j+1).sup.th iterated minimum
mean square error {circumflex over (.epsilon.)}.sub.i.sup.j+1 is
greater than the threshold in (6), and returning to (7) if yes,
otherwise ending the process.
[0101] An angiography image of a patient 1 obtained from an angle
of (50.8.degree. , 30.2.degree.) is shown in FIG. 2, and bifurcate
points are extracted as structural feature points automatically, as
shown in FIG. 3. A cardiac movement period equals 10 frames, and a
sampling rate of an angiography image sequence is 12.5 frames per
second. For angiography time is quite short, structural feature
points A1-A4 with comparatively long residence times are selected
for experiment.
[0102] An angiography image of a patient 4 obtained from an angle
of (30.8.degree. , 15.3.degree.) is shown in FIG. 11, and bifurcate
points are extracted as structural feature points automatically, as
shown in FIG. 12. The patient suffers from arrhythmia, two cardiac
movement periods in a range of the patient's cardiac period are
extracted by the method of the invention, equaling 9 frames and 12
frames respectively, and a sampling rate of an angiography image
sequence is 12.5 frames per second.
[0103] It can be seen from FIGS. 4A, 4B, 13A and 13B that original
motion curves of the feature points are irregular and significantly
influenced by breathing movement, and motion amplitudes are
different for different feature points for the feature points are
located on different parts of the heart with different movements.
FIGS. 5A, 5B, 14A and 14B show spectrums obtained by performing a
discrete Fourier transformation on mixed signals. It can be seen
that each frequency point features an apparent peak, the number of
which illustrates the number of frequency points contained in the
mixed signal. In particular, as a peak appears at a frequency point
0, a translation signal must exist in the mixed signal. FIGS. 6A,
6B, 15A, 15B, 16A and 16B are extracted cardiac motion curves which
show excellent periodicity. FIGS. 7A, 7B, 17A and 17B are extracted
breathing motion curves which show excellent periodicity. It can be
seen that peaks of the feature points occur almost simultaneously
for angiography time is quite short and a corresponding phase
change of each feature point is extremely small. Amplitudes of the
peaks illustrate distribution of the feature points on a vascular
surface, that is, a greater peak corresponds to a location closer
to the lung. FIGS. 8A and 8B are extracted tremor signals which
show that the patient trembled significantly during the treatment
and amplitudes thereof are different with distribution of the
feature points. FIGS. 9A and 9B illustrate extracted translational
motion curves of feature points in FIG. 2 compared with those
obtained by tracking skeleton points manually. It can be seen that
results of the two methods are almost consistent in addition to
errors caused by manual tracking. FIGS. 10A, 10B, 18A and 18B are
residual curves of feature points in FIG. 11 after all movements
being extracted. It can be seen that amplitudes thereof are smaller
than 2 pixels which may be caused by segmentation and skeleton
extraction, and motion components are extracted and separated
completely.
[0104] Different from the patient 1 of FIG. 2, for the patient 4 of
FIG. 11, there exist two different frequencies in a frequency range
of a cardiac signal instead of a high frequency of tremor, which
provides an important reference for diagnosis of a patient.
[0105] Overall, considering qualified feature points (such as
intersections of ribs) do not exist in every angiography image in
single-armed X-ray angiography, the invention extracts multiple
motion components by selecting vascular structural feature points
combined with filtering in frequency domain by Fourier
transformations and optimizing multiple parameters, which is more
flexible and can be widely used in almost all angiography images
compared with obtaining breathing movement or translational
movement merely by manual tracking. Besides, the method of the
invention features a higher level of security and operability
compared with setting labels on tissues adjacent to the heart and
tracking them by related imaging methods, for the labels are
usually intrusive and may damage a human body more or less, and the
whole process of label setting, imaging and eliminating and
breathing movement extraction is complicated which brings troubles
and errors in operation inevitably.
[0106] While particular embodiments of the invention have been
shown and described, it will be obvious to those skilled in the art
that changes and modifications may be made without departing from
the invention in its broader aspects, and therefore, the aim in the
appended claims is to cover all such changes and modifications as
fall within the true spirit and scope of the invention.
* * * * *