CN101059869A - Method for separating blood vessel data in digital coronary angiography - Google Patents

Abstract

The invention discloses a method for dividing vessel data in digit coronary contrast image, comprising that (1), strengthening the frequency domain of the coronary contrast image sequence, to obtain strengthened new image sequence f(x, y, t), (2), extracting object single frame strengthened image f(x, y) from f(x, y, t), to process disc structure low-cap conversion on aorta and tiny vessel, to obtain gray images g1(x, y) and g2(x, y), (3), respectively processing relative threshold comparisons on g1(x, y) and g2(x, y) to extract the vessel in the gray images to obtain binary images B1(x, y) and B2(x, y), (4), combining the binary images to obtain divided image. The invention can completely utilize the motion, structure, size, shape and gray information of vessel, to effective divide contrast image, with simple operation and high efficiency.

Description

The method of in the digital coronary contrastographic picture, cutting apart blood vessel data

Technical field

The invention belongs to the medical imaging technology field, be specifically related to a kind of method of in the digital coronary contrastographic picture, cutting apart blood vessel data.

Background technology

The digital angiographic technology has been used more than 20 year clinical, is the important evidence of cardiovascular and cerebrovascular disease non-invasive diagnosis and interventional therapy surgical navigational.The coronary angiography art is to use the most general being used to diagnose coronary artery disease and assesses the image method of stenosis.A mission critical during contrastographic picture is handled carries out image segmentation exactly, so that the physiological characteristic of blood vessel can clearerly show, subsequent operations such as structure analysis, motion analysis, three-dimensional visualization, and applied researcies such as the image guiding is performed the operation, treatment assessment are all based on image segmentation.Yet, make that the signal to noise ratio (S/N ratio) of heart coronary artery contrastographic picture is very low because tissue thickness's difference, the contrast concentration in the blood of X ray process is inhomogeneous and different background tissue and target are mixed in together.In addition, the institutional framework of human body and shape are very complicated, and sizable difference is arranged between men.Therefore the blood vessel segmentation of contrastographic picture is a very task of difficulty.

C.Kirbas and F.Quek (" A review of vessel extraction techniques andalgorithms; " ACM Computing Surveys, vol.36 pp.81-121 June 2004.) proposes, existing blood vessel segmentation technology roughly can be divided into 6 big classes: mode identification method, based on the method for model, based on method, artificial intelligence approach, neural net method, the tubular articles detection method of following the tracks of, and these traditional blood vessel segmentation methods at contrastographic picture nearly all are to utilize the gray scale or the blood vessel structure information of single image.Because the otherness between the coronary angiography image is very big and the single image signal to noise ratio (S/N ratio) is very low, make the adaptability of these class methods have very big problem and can not extract all blood vessels, especially the low part of contrast.

Summary of the invention

The object of the present invention is to provide a kind of method of cutting apart blood vessel data in the digital coronary contrastographic picture, this method can effectively be cut apart contrastographic picture, and easy and simple to handle, high efficiency.

The method of in the digital coronary contrastographic picture, cutting apart blood vessel data provided by the invention, its step comprises:

(1) to image sequence f ₀(x, y t) carry out frequency domain and strengthen, obtain new sequence f (x, y, t), image sequence f wherein ₀(x, y, frame number t) they are T, the single image size is M * N, 0≤x＜M, 0≤y＜N, 0≤t＜T, x, y, t is integer;

(2) ((x, y), (x, y) carrying out diameter respectively is D to f t) to take out single frames enhancing image f to be split for x, y at new sequence f ₁And D ₂The low cap conversion of disc structure, obtain g ₁(x, y) and g ₂(x, y); Wherein, D ₁＞D ₂, D ₁D ₂Be positive integer, D _MAX＜D ₁≤ (D _MAX+ 5), D _MIN＜D ₂≤ (D _MIN+ (D _MAX-D _MIN)/2), D _MAXBe the wideest diameter coronarius, D _MINBe the narrowest diameter coronarius;

(3) according to following step respectively to g ₁(x, y) and g ₂(x y) carries out threshold ratio, obtains bianry image B ₁(x, y) and B ₂(x, y):

(3.1) (x, y) (x, (x, y), wherein, (x y) is g to g y) to obtain differential chart d with strengthening image f to utilize the relatively lower cap changing image g of formula (I) ₁(x, y) or g ₂(x, y);

(3.2) utilize formula (II) calculated threshold T ₁, T ₁Be d (x, y) mean value of the non-zero pixels value among the figure

$T_1=\frac{\underset{d(x,y)>0}{\mathrm{\Σ}}d(x,y)}{N_d}---\left(\mathrm{II}\right)$

N wherein _dBe d (x, y) number of pixels of non-zero pixels value among the figure;

(3.3) if d (x, y) 〉=T ₁, then this is the blood vessel picture element, if d (x, y)＜T ₁, then this is a background dot, obtains binary map B ₁(x, y) and B ₂(x, y);

(4) stack binary map B ₁(x, y) and B ₂(x y), removes area and is lower than threshold value T ₂The block distortion zone, T ₂Be threshold value big or small according to image and that the blood vessel thickness is chosen, 0＜T ₂≤ (D _MIN* D _MIN), resulting binary map is extraction result coronarius.

At traditional digital coronarogram nearly all is the half-tone information that adopts single image as dividing method, the limited defective of the information of utilizing, the present invention at first utilizes the different motion of coronary artery image sequence medium vessels and background to change, and adopts time domain Fourier high-pass transform to weaken the low frequency ground unrest and strengthen blood vessel.Simultaneously, also consider tubular structure, the diameter information of blood vessel, adopt the multiple dimensioned morphologic filtering that meets blood vessel structure further to strengthen single-frame images to be split.At last, the gray difference consistance between different images according to blood vessel and background is partitioned into bianry image.The inventive method has made full use of motion, structure, yardstick, the shape of blood vessel, the information of gray scale five aspects, good between the different images than big-difference adaptability, the coronary artery segmentation result is complete, for the part a little less than the contrasts such as minute blood vessel or blood vessel end fine performance is arranged also.

Description of drawings

Fig. 1 is the process flow diagram of the inventive method;

Fig. 2 is the contrastographic picture sequence (48*512*512 size, 8) of arteria coroaria sinistra (LCA) tree;

Fig. 3 is the enhancing image sequence of Fig. 2 after through step (1) (β=1/5);

Fig. 4 (a) is the 22nd two field picture that takes out from original sequence (Fig. 2), and Fig. 4 (b) is the 22nd two field picture that takes out from strengthen image sequence (Fig. 3);

Fig. 5 is to be the artery leaching process of example with Fig. 4 (b).Fig. 5 (a) illustrates diameter D ₁=35 disc structure hangs down cap (Bottom-Hat) conversion conversion figure; Fig. 5 (b) illustrates diameter D ₁=5 disc structure hangs down cap (Bottom-Hat) conversion conversion figure; Fig. 5 (c), Fig. 5 (d) is respectively Fig. 5 (a), and Fig. 5 (b) " relative threshold is relatively " is figure as a result;

Reflection shown in Figure 6 in exemplary embodiments of the present invention, use a comparative result of said method and single image dividing method, Fig. 6 (a) is the 22nd two field picture among the original series figure, 6 (b) represent the segmentation result of this method; Fig. 6 (c) expression segmentation result of S.Ehio based on morphological method; Fig. 6 (d) expression result of global threshold dividing method; Dark shape strip zone in the white rectangle frame among Fig. 6 (a) is used to inject the conduit of contrast preparation in order to mark.

Embodiment

The present invention is further detailed explanation below in conjunction with accompanying drawing and example.

As shown in Figure 1, the step of the inventive method comprises:

(1) to image sequence f ₀(x, y t) carry out frequency domain and strengthen this image sequence f ₀(x, y, frame number t) they are T, the single image size is M * N, 0≤x＜M, 0≤y＜N, 0≤t＜T, x, y, t is integer.

Image sequence is carried out the frequency domain enhancing can adopt several different methods, now illustrate wherein a kind of concrete steps of method:

(1.1) to image sequence f ₀(x, y t) carry out M * N discrete Fourier transform (DFT) along time shaft t, obtain frequency-region signal F ₀(k), 0≤k＜T, discrete Fourier transform (DFT) formula are suc as formula (1): wherein j is an imaginary unit for x, y

$F_0(x,y,k)=\frac{1}{T}\underset{t=0}{\overset{T-1}{\mathrm{\Σ}}}{f}_0(x,y,t)e^{-j2\mathrm{\πkt}/T}$

0≤k＜T，0≤x＜M，0≤y＜N (1)

(1.2) to frequency-region signal F ₀(k) employing formula (2) is carried out M * N high-pass filtering along the k axle for x, y, the high-pass filtering function be H (x, y, k)=(1-e ^{-β t}), wherein β is a filter constant given in advance, in order to suppress the ground unrest of some frequency, through repeatedly experiment, the preferable range of β value is 1/8--1/5.

F(x，y，k)＝F ₀(x，y，k)(1-e ^-βm)

0≤k＜T，0≤x＜M，0≤y＜N，

(2)

(1.3) utilize inverse discrete fourier transform (formula (3), wherein j is an imaginary unit) to act on filtered frequency-region signal, obtain new image sequence:

$f(x,y,t)=\underset{k=0}{\overset{T-1}{\mathrm{\Σ}}}\stackrel{\‾}{F}(x,y,k)e^{j2\mathrm{\πk}/T}$ (3)

0≤t＜T，0≤x＜M，0≤y＜N

(2) ((x, y), (x, y) carrying out diameter respectively is D to f t) to take out single frames enhancing image f to be split for x, y at new sequence f ₁And D ₂Low cap (Bottom-Hat) conversion of disc structure, obtain g ₁(x, y) and g ₂(x, y); Wherein, D ₁＞D ₂, D ₁D ₂Be positive integer, D _MAX＜D ₁≤ (D _MAX+ 5), D _MIN＜D ₂≤ (D _MIN+ (D _MAX-D _MIN)/2), D _MAXBe the wideest diameter coronarius, D _MINBe the narrowest diameter coronarius.

It is the morphological operator that the dark structure of image is extracted in a kind of common being used to that morphology hangs down cap (Bottom-Hat) conversion, and this operator is represented with g, is defined as g=I-(Ib), and wherein I is an input picture, and b is a structural element, and fb represents closed operation.

(3) respectively to g ₁(x, y) and g ₂(x y) carries out " relative threshold relatively ", obtains bianry image B ₁(x, y) and B ₂(x, y).Process is as follows:

(3.1) (x, y) (x, (x, y), wherein, (x y) is g to g y) to obtain differential chart d with strengthening image f to utilize the relatively lower cap changing image g of formula (4) ₁(x, y) or g ₂(x, y);

(4)

(3.2) utilize formula (5) calculated threshold T ₁, T ₁Be d (x, y) mean value of the non-zero pixels value among the figure

$T_1=\frac{\underset{d(x,y)>0}{\mathrm{\Σ}}d(x,y)}{N_d}$ (5)

N wherein _dBe d (x, y) number of pixels of non-zero pixels value among the figure.

(3.3) if d (x, y) 〉=T ₁, then this is the blood vessel picture element, if d (x, y)＜T ₁, then this is a background dot, obtains binary map B ₁(x, y) and B ₂(x, y).

(4) binary map merges: stack binary map B ₁(x, y) and B ₂(x y), removes area and is lower than threshold value T ₂The block distortion zone, T ₂Be threshold value big or small according to image and that the blood vessel thickness is chosen, 0＜T ₂≤ (D _MIN* D _MIN), resulting binary map is extraction result coronarius.

Example:

This example image contains dark target and bright background as shown in Figure 2 in the image, dark object representation arteries, the step of this example of detailed description below:

(1) to image sequence f shown in Figure 2 ₀(t) (0≤x＜512,0≤y＜512,0≤t＜48) are carried out frequency domain and are strengthened for x, y, this sequence comprise 48 width of cloth from paradoxical expansion to the cardiac enlargement phase again to the image of paradoxical expansion.Its process is:

(1.1) to image sequence f ₀(x, y t) carry out the discrete Fourier transform (DFT) along time shaft t 512 * 512 times, obtain 48 * 512 * 512 frequency-region signal F that count ₀(x, y, k), and 0≤k＜48, the discrete Fourier transform (DFT) formula is as follows:

$F_0(x,y,k)=\frac{1}{48}\underset{t=0}{\overset{47}{\mathrm{\Σ}}}{f}_0(x,y,t)e^{-j2\mathrm{\πkt}/48}$

0≤k＜48，0≤x＜512，0≤y＜512 (1)

(1.2) to frequency-region signal F ₀(x, y k) adopt (2) formula to carry out 512 * 512 high-pass filterings along the k axle, and the high-pass filtering function is $H(x,y,k)=(1-e^{-\frac{1}{5}t})$

$\stackrel{\‾}{F}(x,y,k)=F_0(x,y,k)(1-e^{-\frac{1}{5}m})$

(2)

(1.3) to filtered frequency-region signal F (x, y k) carry out 512 * 512 inverse discrete fourier transforms (formula (3)), obtain new time domain image sequence as shown in Figure 4:

$f(x,y,t)=\underset{k=0}{\overset{47}{\mathrm{\Σ}}}\stackrel{\‾}{F}(x,y,k)e^{-j2\mathrm{\πk}/48}$

0≤t＜48，0≤x＜M，0≤y＜N (3)

(2) from strengthen image sequence, (x, y), (x y) carries out diameter D respectively to f to take out the 22nd frame single-frame images f to be split ₁=35, D ₂=5 disc structure hangs down cap (Bottom-Hat) conversion, obtains g ₁(x, y) (Fig. 5 (a)) and g ₂(x, y) Fig. 5 (b).

(3) respectively to g ₁(x, y) and g ₂(x y) carries out " relative threshold relatively " and obtains bianry image B ₁(x, y) (Fig. 5 (c)) and B ₂(x, y) (Fig. 5 (d)).

(3.1) (x, y) (x, (x, y), wherein, (x y) is g to g y) to obtain differential chart d with strengthening image f to utilize the relatively lower cap changing image g of formula (4) ₁(x, y) or g ₂(x, y);

(4)

(3.2) utilize formula (5) calculated threshold T ₁, T ₁Be d (x, y) mean value of the non-zero pixels value among the figure

$T_1=\frac{\underset{d(x,y)>0}{\mathrm{\Σ}}d(x,y)}{N_d}$ (5)

N wherein _dBe d (x, y) number of pixels of non-zero pixels value among the figure.

(3.3) if d (x, y) 〉=T ₁, then this is the blood vessel picture element, if d (x, y)＜T ₁, then this is a background dot, obtains binary map B ₁(x, y) and B ₂(x, y).

(4) binary map merges: stack binary map B ₁(x, y) and B ₂(x y), removes area and is lower than threshold value T ₂=20 block distortion zone, resulting binary map (Fig. 6 (b)) is extraction result coronarius.

Found out that by Fig. 6 (b) segmentation result (Fig. 6 (b)) that adopts the present invention to obtain is cut apart complete, residual noise is few, and the static contrast preparation duct portion overwhelming majority among Fig. 6 (a) in the white rectangle frame is removed; But the segmentation result (Fig. 6 (c), Fig. 6 (d)) that obtains according to the half-tone information and the structural information method of single image only, blood vessel structure is imperfect: the weak little blood vessel of contrast is removed, cuts apart the very big noise of introducing, dwell catheters is present in the segmentation result.

Adopt the inventive method, effectively raise the contrast of image, can from the low signal-to-noise ratio contrastographic picture, automatically extract whole arterial tree, fine performance is also arranged for minute blood vessel part or blood vessel end part.

Claims (2)

1, a kind of method of in the digital coronary contrastographic picture, cutting apart blood vessel data, its step comprises:

(1) to image sequence f ₀(x, y t) carry out frequency domain and strengthen, obtain new sequence f (x, y, t), image sequence f wherein ₀(x, y, frame number t) they are T, the single image size is M * N, 0≤x＜M, 0≤y＜N, 0≤t＜T, x, y, t is integer;

(2) ((x, y), (x, y) carrying out diameter respectively is D to f t) to take out single frames enhancing image f to be split for x, y at new sequence f ₁And D ₂The low cap conversion of disc structure, obtain g ₁(x, y) and g ₂(x, y); Wherein, D ₁＞D ₂, D ₁D ₂Be positive integer, D _MAX＜D ₁≤ (D _MAX+ 5), D _MIN＜D ₂≤ (D _MIN+ (D _MAX-D _MIN)/2), D _MAXBe the wideest diameter coronarius, D _MINBe the narrowest diameter coronarius;

(3) according to following step respectively to g ₁(x, y) and g ₂(x y) carries out threshold ratio, obtains bianry image B ₁(x, y) and B ₂(x, y):

(3.1) (x, y) (x, (x, y), wherein, (x y) is g to g y) to obtain differential chart d with strengthening image f to utilize the relatively lower cap changing image g of formula (I) ₁(x, y) or g ₂(x, y);

(3.2) utilize formula (II) calculated threshold T ₁, T ₁Be d (x, y) mean value of the non-zero pixels value among the figure

$T_1=\frac{\underset{d(x,y)>0}{\mathrm{\Σ}}d(x,y)}{N_d}...\left(\mathrm{II}\right)$

N wherein _dBe d (x, y) number of pixels of non-zero pixels value among the figure;

(3.3) if d (x, y) 〉=T ₁, then this is the blood vessel picture element, if d (x, y)＜T ₁, then this is a background dot, obtains binary map B ₁(x, y) and B ₂(x, y);

(4) stack binary map B ₁(x, y) and B ₂(x y), removes area and is lower than threshold value T ₂The block distortion zone, T ₂Be threshold value big or small according to image and that the blood vessel thickness is chosen, 0＜T ₂≤ (D _MIN* D _MIN), resulting binary map is extraction result coronarius.

2, method according to claim 1 is characterized in that: step (1) is handled by following process:

(1.1) to image sequence f ₀(x, y t) carry out M * N discrete Fourier transform (DFT) along time shaft t, obtain frequency-region signal F ₀(k), 0≤k＜T, discrete Fourier transform (DFT) formula are suc as formula (III): wherein j is an imaginary unit for x, y

$F_0(x,y,k)=\frac{1}{T}\underset{t=0}{\overset{T-1}{\mathrm{\Σ}}}{f}_0(x,y,t)e^{-2\mathrm{\πkt}/T}$

0≤k＜T，0≤x＜M，0≤y＜N (III)

(1.2) to frequency-region signal F ₀(k) employing formula (IV) is carried out M * N high-pass filtering along the k axle for x, y, the high-pass filtering function be H (x, y, k)=(1-e ^{-β t}), wherein β is a filter constant given in advance,

F(x，y，k)＝F ₀(x，y，k)(1-e ^-βm)

(1.3) utilize formula (V) to carry out inverse discrete fourier transform, act on filtered frequency-region signal, acquisition new sequence f (x, y, t):

$f(x,y,t)=\underset{k=0}{\overset{T-1}{\mathrm{\Σ}}}\stackrel{\‾}{F}(x,y,k)e^{j2\mathrm{\πk}/T}...\left(V\right)$

0≤t＜T，0≤x＜M，0≤y＜N。

Images