CN105741299A - Coronary artery CT angiography image segmentation method - Google Patents

Abstract

The invention provides a coronary artery CT angiography image segmentation method. The method comprises the following steps of firstly preprocessing a CTA image to obtain a binary image, wherein the preprocessing comprises median filtering and threshold segmentation; secondly marking all connected regions in the binary image and calculating characteristic information of each connected region; thirdly judging whether the coronary artery is split or not by utilizing a characteristic matching mechanism; and finally tracking a target by utilizing the characteristic matching mechanism in more than ten frames of initial images and then tracking the target by utilizing Kalman. The target is tracked through combination of the characteristic matching mechanism and the Kalman filter, so that the tracking of a target region in part of special forms can be realized and the tracking accuracy of the target region is improved. Moreover, in the whole tracking process, the characteristic matching mechanism judges whether the coronary artery in the image is split or not, so that whether the coronary artery is split or not and which of split parts can be accurately judged and the tracking accuracy can be further ensured.

Description

A kind of coronary artery CT angiographic image dividing method

Technical field

The present invention relates to computer assisted image processing technical field, specifically the coronary artery CT angiographic image dividing method in a kind of medical image.

Background technology

Coronary heart disease is the disease that worldwide mortality rate is the highest, and its sickness rate presents lasting ascendant trend, and healthy to people constitutes serious threat.As effective non-invasive means of diagnosis of coronary heart disease, the development of multislice spiral CT angiography (ComputedTomographyAngiography, CTA) technology is comparatively rapid.The segmentation of CTA image is the important foundation to the reconstruct of arteria coronaria (i.e. coronary artery) 3-D geometric model, and the diagnostic mode that threedimensional model can provide more accurately quickly for Cardiologists, repeatability is higher.Therefore, arteria coronaria profile is extracted exactly based on arteria coronaria CTA sequence image, it it is the important clinical assistant analysis instrument of coronary stricture, and can provide calcification degree, speckle burden and the quantitative analysis of stenosis, therefore become the focus of field of medical image processing research based on the blood vessel segmentation of CTA image.

Currently, domestic and international researcher has done substantial amounts of research about blood vessel structure sequences segmentation especially coronarius.Basis thought is also all the seriality utilizing image interframe.But, most blood vessel segmentation method heavy dependence and user's is mutual, for instance manual definition starts and end point and vessel directions, or is manually inserted intermediate point and narrows the gap.These methods make the accuracy of image segmentation result be strongly depend on the interaction of user.In addition, also having certain methods is as new seed points to realize sequences segmentation using the barycenter of target area or outline projection, although these methods are capable of the CTA Segmentation of Image Sequences of certain organ, tissue, but they or be poorly suited for use in the situation of change in topology, or there is very strong dependency for consecutive image spacing, all can not be well adapted for CTA figure coronarius.

Summary of the invention

It is an object of the invention to provide a kind of coronary artery CT angiographic image dividing method, with solve existing CTA image partition method can not well for the complex fine problem coronarius of segmenting structure.

The present invention is achieved in that a kind of coronary artery CT angiographic image dividing method, comprises the steps:

A, original corona arterial CT angiography: evaluation image is carried out pretreatment；

Specifically: the coronary artery CT angiographic image of all arranged in sequences is carried out two dimension median filter, adopt maximum variance between clusters that filtered image is carried out Threshold segmentation afterwards, obtain bianry image；

Adopt maximum variance between clusters that filtered image is carried out Threshold segmentation, specifically according to formula below, filtered image carried out Threshold segmentation:

$f_t(x,y)=\left\{\begin{array}{c}1,f(x,y)\&GreaterEqual;\mathrm{\&lambda;}t\\ 0,f(x,y)<\mathrm{\&lambda;}t\end{array}\right\}$

Wherein, t is segmentation threshold, and λ is for adjusting parameter；λ preferably 2.3；

B, all connected regions in pretreated bianry image are carried out labelling, and obtain the barycenter of each connected region, area and three characteristic informations of eccentricity；

C, from the first two field picture, calibrate target area, in the target area of the first two field picture, include aorta and two targets of coronary artery；

D, from the second two field picture, utilize whether the coronary artery that characteristic matching mechanism judges in present image according to the coronary artery in previous frame image divides；

E, choose the connected region that window inner area corresponding with coronary artery in previous frame image in present image is minimum, and judge that whether the area of the minimum connected region of this area is more than preset area threshold value C, if so, then execution step f；If it is not, then perform step g；

F, utilizing characteristic matching mechanism according to the target in previous frame image, the target in present image to be tracked, the target in described present image is coronary artery or coronary artery and aorta；

G, utilizing Kalman filter according to the target in previous frame image, the target in present image to be tracked, the target in described present image is coronary artery or coronary artery and aorta；

H, repetition step d～g, until completing the tracking to all image internal objects, namely achieve segmentation coronarius.

Also comprise the steps: after step h in each image using the barycenter of target that tracks as seed points, adopt Laplce's level-set segmentation methods that original corona arterial CT angiography: evaluation image is split, to extract coronary artery profile.

Step d utilizing, whether the coronary artery that characteristic matching mechanism judges in present image according to the coronary artery in previous frame image divides, specifically include following steps:

D1, choosing the window corresponding with coronary artery previous frame image in present image as the first trace regions, the area of this first trace regions is area coronarius slightly larger than in previous frame image；

D2, the area of all connected regions in area coronarius in previous frame image and present image the first trace regions is compared respectively, if there is at least one connected region, its area and the difference of area coronarius in previous frame image are less than the first area setting value, its barycenter and centroid distance coronarius in previous frame image are less than the first distance setting value, then the coronary artery judged in present image does not divide；Otherwise, then the coronary artery judged in present image there occurs division.

Step f utilizes characteristic matching mechanism according to the target in previous frame image, the target in present image is tracked, specifically includes following steps:

Whether f1, the target first determined whether in previous frame image comprise aorta, if not comprising aorta, then directly perform step f4；Otherwise perform step f2；

F2, choosing the window corresponding with aorta previous frame image in present image as the second trace regions, the area of this second trace regions is aortal area slightly larger than in previous frame image；If there is at least one connected region in present image, its area and the difference of aortal area in previous frame image are less than second area setting value, its barycenter and aortal centroid distance in previous frame image are less than second distance setting value, then judge that present image exists aorta, now perform step f3；Otherwise, it is determined that present image is absent from aorta, and perform step f4；

F3, choose in the second trace regions all connected regions and calculate each connected region cost function V (i, j):

V (i, j)=al (D (i, j))+bl (A (i, j))+cl (M (i, j)) (1)

Wherein, $l\left(D\right(i,j\left)\right)=\frac{D(i,j)}{\underset{m=1}{\overset{n}{\mathrm{\&Sigma;}}}D(i,m)}---\left(2\right)$

$D(i,j)=\sqrt{{(x_k^i-x_{k+1}^j)}^2+{(y_k^i-y_{k+1}^j)}^2}---\left(3\right)$

For the center-of-mass coordinate of kth frame image internal object,For the center-of-mass coordinate of jth connected region to be matched in kth+1 two field picture,For center-of-mass coordinate and the center-of-mass coordinate distance sum of all connected regions in kth+1 two field picture of kth frame image internal object, n is the number of all connected regions in kth+1 two field picture；(D (i, j)) represents jth connected region to be matched and kth frame image internal object similarity in position in kth+1 two field picture to l；

$l\left(A\right(i,j\left)\right)=\frac{\left|{\mathrm{Area}}_i-{\mathrm{Area}}_j\right|}{\underset{m=1}{\overset{n}{\mathrm{\&Sigma;}}}\left|{\mathrm{Area}}_i-{\mathrm{Area}}_m\right|}---\left(4\right)$

Area_iFor the area of kth frame image internal object, Area_jFor the area of jth connected region to be matched in kth+1 two field picture,For in area and kth+1 two field picture of kth frame image internal object the difference in areas of all connected regions and；(A (i, j)) represents jth connected region to be matched and kth frame image internal object similarity on area in kth+1 two field picture to l；

$l\left(M\right(i,j\left)\right)=\frac{M(i,j)}{\underset{m=1}{\overset{n}{\mathrm{\&Sigma;}}}M(i,m)}---\left(5\right)$

M (i, j)=L (i)/S (i)-L (j)/S (j) (6)

The major axis of L (i), S (i) respectively kth frame image internal object and minor axis length, the major axis of jth connected region to be matched and minor axis length in L (j), S (j) respectively kth+1 two field picture,Difference is asked to sue for peace respectively again for the ratio of semi-minor axis length of all connected regions in the ratio of semi-minor axis length of kth frame image internal object and kth+1 two field picture, (M (i, j)) represents jth connected region to be matched and kth frame image internal object similarity in eccentricity in kth+1 two field picture to l；

In formula (1), a, b and c are weight coefficient；Preferably, a is 0.5, b be 0.3, c is 0.2；

In kth+1 two field picture, (i, j) minimum connected region is the target tracked to cost function V, and the target tracked in this step is aorta；

F4, choosing the window corresponding with coronary artery previous frame image in present image as the first trace regions, the area of this first trace regions is area coronarius slightly larger than in previous frame image；

If coronary artery does not divide in present image, then coronary artery is tracked in the first trace regions of present image according to step f3；

If coronary artery there occurs division in present image, then find out the connected region that some barycenter are nearest from coronary artery barycenter in previous frame image, and the area of these connected regions and the difference of area coronarius is less than the 3rd area setting value with previous frame image, these connected regions found are the coronary artery tracked.

In step g, its system mode of Kalman filter is four dimensional vectors X (k)=(x_k,y_k,xv_k,yv_k)^T, externally input is U (k)=[xa_k,ya_k]^T；

Wherein: x_k,y_kThe respectively transverse and longitudinal coordinate of target centroid, xv_k,yv_kRepresent target speed in X-axis and Y direction, xa respectively_k,ya_kRepresent target acceleration in X-axis and Y direction respectively；The interval of adjacent two two field pictures is T and target does uniformly accelerated motion；

State-transition matrix A (k), systematic parameter B (k), observing matrix H (k) are respectively as follows:

$A\left(k\right)=\left[\begin{array}{cccc}1& T& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& T\\ 0& 0& 0& 1\end{array}\right]$ $B\left(k\right)=\left[\begin{array}{cc}0.5T^2& 0\\ T& 0\\ 0& 0.5T^2\\ 0& T\end{array}\right]$ $H\left(k\right)=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\end{array}\right]$

Covariance matrix is then respectively as follows:

$Q=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$ $R=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]$

Due to coronary morphology change complexity, there is subbranch structure unfixed, that curvature is complicated, therefore it is difficult with a kind of independent tracing algorithm and realize the whole of all targets are followed the trail of；It addition, the continuous division of target is also the big difficult point followed the trail of, and this is the independent insurmountable problem of Kalman (Kalman) algorithm.Therefore, by the mechanism of characteristic matching, the present invention judges whether the coronary artery in image (abbreviation arteria coronaria) there occurs division, judge whether arteria coronaria there occurs division only exactly, the accuracy that guarantee is followed the trail of.And, the present invention, when long-pending relatively big (more than preset area threshold value C) after initial (in about tens two field pictures) arteria coronaria does not divide and divides, adopts characteristic matching mechanism that target is tracked；And when amassing less after arteria coronaria divides, be converted to and be tracked by Kalman filter.This be due to: consider the change of target area form in tens two field pictures started, simultaneously in order to obtain initial parameter more accurately, tracking at about tens initial two field picture internal objects temporarily obtains target location without Kalman filter, but tries to achieve positional information by the mechanism of characteristic matching.When the area of image internal object is reduced to certain value, change in location amplitude increases, and now, is difficult to provide the parameter of a suitable control site again, such too for area, so can be substantially reduced the success rate of tracking.Based on this, the present invention gives preset area threshold value C, owing to area information role in this ring is obvious not as position feature, so using the C critical parameters as judgement target area (referred to herein as arteria coronaria) size.Once the area of target area is less than C, just tracking being transformed under Kalman filter and carry out, utilize previous frame Image Acquisition target position information to solve and obtain speed and acceleration information, the information finally these solved is as the initial value of Kalman filter.Hereafter, the mechanism of characteristic matching only plays the effect that target area splitting status judges, coordinates Kalman filter assistance tracking, compensates its defect, thus realizing coronarius following the trail of accurately.

The difficult problem that arteria coronaria form constantly divides can be solved by the present invention；And, based on the feature of arteria coronaria its own shape and motion, make characteristic matching mechanism and Kalman filter combine and target is tracked, it is possible to solve the tracking of target area under part specific form, improve the tracking accuracy rate of target area.

Accompanying drawing explanation

Fig. 1 is the method flow diagram of the present invention.

Detailed description of the invention

The embodiment of the present invention provides coronary artery sequences segmentation by 128 row's spiral CT scans, and has carried out three-dimensionalreconstruction after singulation, judge that the identification etc. of narrow location and degree, speckle provides basis for being used for observing later.The present invention is described in detail below.

As it is shown in figure 1, the present invention comprises the steps:

A, acquisition 128 row's CT spiral image.

128 rows' CT spiral image (image indicated by right side arrow is one of them image) obtained are original coronary artery CT angiographic image, and these CT angiogram seem arrange according to acquired time sequencing.First having to these images are carried out pretreatment, preprocessing process includes this step and following medium filtering and Threshold segmentation.In this step, it is necessary to change window width, window place value, this be due to: same CTA scans aspect, and the window width of setting, window place value difference can observe the gray scale image of different tissues, organ structure；And the window position of CTA machine default setting, window width value stick together when can cause following the trail of segmentation situation.Therefore the two parameter value is adjusted to 275,173 by original 50,350 by the present invention, thus generating new process data (in corresponding diagram 1 image) indicated by " 128 row's CT spiral image " right side arrow, and then coronary artery region is stripped out from atrium it is tracked segmentation again.

B, medium filtering.

For the impulsive noise suppressing CTA image scanning to produce after entering computer, make image more smooth, therefore CTA sequence image is carried out medium filtering by this step.According to image size, the present embodiment adopts 3*3 template that all images are carried out two dimension median filter.

C, Threshold segmentation.

Utilize maximum between-cluster variance (OTSU) method that the image after medium filtering is carried out Threshold segmentation, obtain bianry image.OTSU method is to seek segmentation result according to the following formula:

$f_t(x,y)=\left\{\begin{array}{c}1,f(x,y)\&GreaterEqual;\mathrm{\&lambda;}t\\ 0,f(x,y)<\mathrm{\&lambda;}t\end{array}\right\}$

Wherein, f_t(x, y) for the bianry image (image see indicated by " medium filtering and Threshold segmentation " right side arrow in figure) after segmentation, t is segmentation threshold.Owing to there is multiple connected region in image, and gray scale degree difference to some extent, therefore the present invention is multiplied by one before the segmentation threshold t tried to achieve and adjusts the result after parameter lambda as optimal threshold, and during λ=2.3, effect is best.

D, mark all connected regions and obtain the characteristic information of connected region.

First all connected regions in the bianry image obtained are carried out label, after the good all of connected region of labelling, the characteristic information of each connected region can be obtained, here include the barycenter that can characterize positional information, characterize the area of size and characterize eccentricity (such as can obtain eccentricity, focal length, major axis, the short axle etc.) feature etc. of shape.Utilizing these characteristic informations can judge whether target area (referring to arteria coronaria) divides exactly on the one hand, effectively utilizing these characteristic informations can also carry out information matches on the other hand, assists tracking coronarius.Meanwhile, the present invention, by the barycenter of the connected region tracked, mating area growth algorithm, is partitioned into target area, improves segmentation precision.So, the characteristic information obtaining all connected regions is by premise preparation and the Data safeguard of tracking and segmentation.

E, feature based matching mechanisms judge whether arteria coronaria divides.

By the mechanism of characteristic matching, the present invention judges whether the arteria coronaria in image there occurs division, judge whether arteria coronaria there occurs division only exactly, the accuracy that guarantee is followed the trail of.

Characteristic matching mechanism is utilized to judge whether arteria coronaria divides, specifically: choosing the window corresponding with coronary artery previous frame image in present image as the first trace regions, the area of this first trace regions is area coronarius slightly larger than in previous frame image；Area coronarius in previous frame image is compared respectively with the area of all connected regions in present image the first trace regions, if there is at least one connected region, its area and the difference of area coronarius in previous frame image are less than the first area setting value, its barycenter and centroid distance coronarius in previous frame image are less than the first distance setting value, then the coronary artery judged in present image does not divide；Otherwise, then the coronary artery judged in present image there occurs division.For not there is the coronary artery of division, area equation coronarius or be more or less the same in its area and previous frame image；For the coronary artery after dividing, its area is the area coronarius much smaller than in previous frame image necessarily, and centroid position also can change a lot.

The tracking of f, feature based coupling.

Consider the change of arteria coronaria form in tens two field pictures (such as the 12nd frame～the 15th frame) started, simultaneously in order to obtain initial parameter more accurately, tracking at tens initial two field picture internal objects temporarily obtains target location without Kalman filter, but tries to achieve positional information by the mechanism of characteristic matching.When the area of image internal object is reduced to certain value, change in location amplitude increases, and now, is difficult to provide the parameter of a suitable control site again, such too for area, so can be substantially reduced the success rate of tracking.Based on this, The present invention gives preset area threshold value C, owing to area information role in this ring is obvious not as position feature, so using the C critical parameters as judgement target area (referred to herein as arteria coronaria) size.Once the area S of target area is less than C, just tracking being transformed under Kalman filter and carry out, utilize previous frame Image Acquisition target position information to solve and obtain speed and acceleration information, the information finally these solved is as the initial value of Kalman filter.Hereafter, the mechanism of characteristic matching only plays the effect that target area splitting status judges, coordinates Kalman filter assistance tracking, compensates its defect, thus realizing coronarius following the trail of accurately.

Tracking to present image is to mate or Kalman by feature based, need first to choose the connected region that window inner area corresponding with coronary artery in previous frame image in present image is minimum, and judge that whether the area of the minimum connected region of this area is more than preset area threshold value C, if so, then select feature based matching mechanisms that target is tracked；If it is not, then select based on Kalman filtering, target to be tracked.

In this step, target is tracked by feature based matching mechanisms, specifically includes following steps:

Whether f1, the target first determined whether in previous frame image comprise aorta, if not comprising aorta, then directly perform step f4；Otherwise perform step f2.

Target in present image is tracked, is based on the target in previous frame image and carries out.For the first two field picture, target area thereon is artificial demarcation out.In the first two field picture, in target area, contain aorta and two targets of arteria coronaria.As time goes on, in image below, aorta can disappear, and arteria coronaria can divide.Therefore, in the process that aorta and arteria coronaria are tracked, should considering that aorta is either with or without disappearance in time, arteria coronaria is either with or without division.

This step first determining whether, whether the target in previous frame image comprises aorta, if just without aorta in previous frame image, then also necessarily without aorta in present image and later image, as long as follow-up, arteria coronaria is tracked, namely performs step f4；If aorta yet suffers from previous frame image, then perform step f2.

F2, choosing the window corresponding with aorta previous frame image in present image as the second trace regions, the area of this second trace regions is aortal area slightly larger than in previous frame image.Judge in present image, whether there is at least one such connected region, its area and the difference of aortal area in previous frame image are less than second area setting value, its barycenter and aortal centroid distance in previous frame image are less than second distance setting value, if there is at least one such connected region, then illustrate present image yet suffers from aorta, next perform step f3；If being absent from such connected region, then illustrate present image is absent from aorta, perform step f4.

F3, choose in the second trace regions all connected regions and calculate each connected region cost function V (i, j).Cost function V (i, expression j) is:

V (i, j)=al (D (i, j))+bl (A (i, j))+cl (M (i, j)) (1)

Wherein, $l\left(D\right(i,j\left)\right)=\frac{D(i,j)}{\underset{m=1}{\overset{n}{\mathrm{\&Sigma;}}}D(i,m)}---\left(2\right)$

$D(i,j)=\sqrt{{(x_k^i-x_{k+1}^j)}^2+{(y_k^i-y_{k+1}^j)}^2}---\left(3\right)$

For the center-of-mass coordinate of kth frame image (i.e. previous frame image) internal object,Center-of-mass coordinate for kth+1 two field picture (i.e. present image) interior jth connected region to be matched, (i j) is the Euclidean distance (namely in present image distance) between jth connected region barycenter to be matched and previous frame image internal object barycenter of two targets to be matched to D；For center-of-mass coordinate and the center-of-mass coordinate distance sum of all connected regions in kth+1 two field picture of kth frame image internal object, n is the number of all connected regions in kth+1 two field picture；L (D (i, j)) representing jth connected region to be matched and kth frame image internal object similarity in position in kth+1 two field picture, (D (i, j)) is more little for l, illustrate position closer to, the probability of coupling is more big.

$l\left(A\right(i,j\left)\right)=\frac{\left|{\mathrm{Area}}_i-{\mathrm{Area}}_j\right|}{\underset{m=1}{\overset{n}{\mathrm{\&Sigma;}}}\left|{\mathrm{Area}}_i-{\mathrm{Area}}_m\right|}---\left(4\right)$

Area_iFor the area of kth frame image internal object, Area_jFor the area of jth connected region to be matched in kth+1 two field picture,For in area and kth+1 two field picture of kth frame image internal object the difference in areas of all connected regions and；(A (i, j)) represents jth connected region to be matched and kth frame image internal object similarity on area in kth+1 two field picture to l；L (A (i, j)) is more little, illustrate area closer to, the probability of coupling is more big.

$l\left(M\right(i,j\left)\right)=\frac{M(i,j)}{\underset{m=1}{\overset{n}{\mathrm{\&Sigma;}}}M(i,m)}---\left(5\right)$

M (i, j)=L (i)/S (i)-L (j)/S (j) (6)

The major axis of L (i), S (i) respectively kth frame image internal object and minor axis length, the major axis of jth connected region to be matched and minor axis length in L (j), S (j) respectively kth+1 two field picture,Difference is asked to sue for peace respectively again for the ratio of semi-minor axis length of all connected regions in the ratio of semi-minor axis length of kth frame image internal object and kth+1 two field picture, l (M (i, j) jth connected region to be matched and kth frame image internal object similarity in eccentricity in kth+1 two field picture) are represented, l (M (i, j)) more little, illustrate eccentricity closer to, the probability of coupling is more big.

In the present invention, centroid position, area and eccentricity these three characteristic information characterize the ability of target area and inconsistent, it is therefore desirable to correspondingly give a weight coefficient in cost function, i.e. a, b and c in formula (1).The size of these weight coefficients is to be determined by the tracking effect tested, and by the tracking of several groups of sequence images be may determine that, when working as a=0.5, b=0.3, c=0.2, the situation of omission is minimum, and the accuracy of tracking is corresponding the highest.From these three weight coefficient it can also be seen that tracking coronarius, centroid position and area characteristic information is bigger compared to the effect that shape facility information plays in the matching process for the present invention.So, select suitable feature to mate according to target information, it is possible to farthest to ensure the accuracy rate followed the trail of.Meanwhile, the mechanism that division judges also is determined by these parameters.

At the cost function V obtaining all connected regions, (i, after j), (i, j) minimum connected region is the target tracked in present image, and the target tracked in this step is aorta to choose cost function V.

F4, choosing the window corresponding with coronary artery previous frame image in present image as the first trace regions, the area of this first trace regions is area coronarius slightly larger than in previous frame image.

If coronary artery does not divide in present image, then coronary artery is tracked in the first trace regions of present image according to step f3.

If coronary artery there occurs division in present image, then find out the connected region that some barycenter are nearest from coronary artery barycenter in previous frame image, and the area of these connected regions and the difference of area coronarius is less than the 3rd area setting value with previous frame image, these connected regions found are the coronary artery tracked.

Tracking under g, Kalman filtering.

Kalman filter is a linear filtering recurrence device.As the pith of Kalman filtering, by observation data, state estimation infers that stochastic variable finally realizes estimation and the prediction of real-time running state on the basis of the method such as least-squares estimation, linear minimum-variance estimation.

If the position that state parameter is target and speed.The system mode of definition Kalman filter is four dimensional vectors X (k)=(x_k,y_k,xv_k,yv_k)^T, externally input is U (k)=[xa_k,ya_k]^T.Wherein: x_k,y_kThe respectively transverse and longitudinal coordinate of target centroid, xv_k,yv_kRepresent target speed in X-axis and Y direction, xa respectively_k,ya_kRepresent target acceleration in X-axis and Y direction respectively；The interval of adjacent two two field pictures is T and target does uniformly accelerated motion.

Definition status shift-matrix A (k), systematic parameter B (k), observing matrix H (k) are respectively as follows:

$A\left(k\right)=\left[\begin{array}{cccc}1& T& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& T\\ 0& 0& 0& 1\end{array}\right]$ $B\left(k\right)=\left[\begin{array}{cc}0.5T^2& 0\\ T& 0\\ 0& 0.5T^2\\ 0& T\end{array}\right]$ $H\left(k\right)=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\end{array}\right]$

Covariance matrix is then respectively as follows:

$Q=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$ $R=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]$

Repeat step e～g, the tracking to all image internal objects can be completed, thus realizing segmentation coronarius.

H, dividing method based on Laplce's level set.

In step d by each connected component labeling of each image out, the characteristic information of all connected regions is namely obtained.Assuming P representative image matrix, the result after labelling is as follows:

$P=\left[\begin{array}{ccccccccc}1& 1& 1& 0& 0& 0& 0& 0& 0\\ 1& 1& 1& 0& 0& 0& 0& 0& 0\\ 1& 1& 1& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 2& 2& 2& 2\\ 0& 0& 0& 0& 0& 2& 2& 2& 2\\ 0& 0& 3& 3& 0& 0& 0& 0& 0\\ 0& 3& 3& 3& 3& 0& 0& 0& 0\\ 0& 3& 3& 3& 3& 0& 0& 0& 0\\ 0& 0& 3& 3& 0& 0& 0& 0& 0\end{array}\right]$

For convenient explanation, here respectively with a square, rectangle, circle as the connected region of institute's labelling, it is evident that in image array P, the square indicia in the upper left corner is No. 1 region, right row rectangle marked placed in the middle is No. 2 regions, and the circle mark near the lower left corner is No. 3 regions.If when following the trail of result and determining that new target area is No. 3 regions, other connected regions being reset, the two dimension segmentation result coronarius wanted can be obtained.

The segmentation result that but now obtains is based on the premise of global threshold segmentation and extracts and obtain, and segmentation precision is a bit weaker.Therefore, in this step after tracking new target area, ask for the barycenter in this region as new seed points, then the dividing method initial threshold segmentation method of replacement based on Laplce's level set extracts objective contour, just can obtain the segmentation result that precision is higher.

In this dividing method, coefficient of curvature CurSca is set to 1.0；Specific conductance is set to 2.0, TimeStep and is set to 0.125；The number of iterator selects 10, and threshold range is set to (150,600).Under the setting based on these parameters, Laplce's level set algorithm high-precision can be partitioned into target area.

I, obtain threedimensional model coronarius.

After obtaining comparatively complete accurate two dimension segmentation result coronarius, it is necessary to they are drawn and becomes threedimensional model.In the present invention, rendering algorithm selects MarchingCubes, i.e. MC algorithm.The ultimate principle of MarchingCubes algorithm is that first each voxel processes, and find out the voxel wherein intersected with contour surface.Being obtained the intersection point of whole cube limit and contour surface by the method for interpolation, be connected contour surface and cubical intersection point then in conjunction with the summit of voxel and the position of contour surface by certain form, thus forming contour surface, approaching expression as its in this cube.Finally obtain the draw about mating surface is drawn again after the relevant parameter of contour surface in volume data field and go out contour surface, the threedimensional model of arteria coronaria can be obtained.In the algorithm, DataExtent is set to (0,512,0,512)；DataSpacing is set to (0.4297,0.4297,0.6)；RadiusFactors is set to (222,222,100)；TargetReduction is set to 0.5；Iterations is set to 500；SpecularPower is set to 5.

Claims (8)

1. a coronary artery CT angiographic image dividing method, is characterized in that, comprises the steps:

A, original corona arterial CT angiography: evaluation image is carried out pretreatment；

Specifically: the coronary artery CT angiographic image of all arranged in sequences is carried out two dimension median filter, adopt maximum variance between clusters that filtered image is carried out Threshold segmentation afterwards, obtain bianry image；

B, all connected regions in pretreated bianry image are carried out labelling, and obtain the barycenter of each connected region, area and three characteristic informations of eccentricity；

C, from the first two field picture, calibrate target area, in the target area of the first two field picture, include aorta and two targets of coronary artery；

D, from the second two field picture, utilize whether the coronary artery that characteristic matching mechanism judges in present image according to the coronary artery in previous frame image divides；

E, choose the connected region that window inner area corresponding with coronary artery in previous frame image in present image is minimum, and judge that whether the area of the minimum connected region of this area is more than preset area threshold value C, if so, then execution step f；If it is not, then perform step g；

F, utilizing characteristic matching mechanism according to the target in previous frame image, the target in present image to be tracked, the target in described present image is coronary artery or coronary artery and aorta；

G, utilizing Kalman filter according to the target in previous frame image, the target in present image to be tracked, the target in described present image is coronary artery or coronary artery and aorta；

H, repetition step d～g, until completing the tracking to all image internal objects, namely achieve segmentation coronarius.

2. coronary artery CT angiographic image dividing method according to claim 1, it is characterized in that, also comprise the steps: after step h in each image using the barycenter of target that tracks as seed points, adopt Laplce's level-set segmentation methods that original corona arterial CT angiography: evaluation image is split, to extract coronary artery profile.

3. coronary artery CT angiographic image dividing method according to claim 1, is characterized in that, utilizes whether the coronary artery that characteristic matching mechanism judges in present image according to the coronary artery in previous frame image divides in step d, specifically:

D1, choosing the window corresponding with coronary artery previous frame image in present image as the first trace regions, the area of this first trace regions is area coronarius slightly larger than in previous frame image；

D2, the area of all connected regions in area coronarius in previous frame image and present image the first trace regions is compared respectively, if there is at least one connected region, its area and the difference of area coronarius in previous frame image are less than the first area setting value, its barycenter and centroid distance coronarius in previous frame image are less than the first distance setting value, then the coronary artery judged in present image does not divide；Otherwise, then the coronary artery judged in present image there occurs division.

4. coronary artery CT angiographic image dividing method according to claim 1, is characterized in that, utilizes characteristic matching mechanism according to the target in previous frame image, the target in present image to be tracked in step f, specifically:

Whether f1, the target first determined whether in previous frame image comprise aorta, if not comprising aorta, then directly perform step f4；Otherwise perform step f2；

F2, choosing the window corresponding with aorta previous frame image in present image as the second trace regions, the area of this second trace regions is aortal area slightly larger than in previous frame image；If there is at least one connected region in present image, its area and the difference of aortal area in previous frame image are less than second area setting value, its barycenter and aortal centroid distance in previous frame image are less than second distance setting value, then judge that present image exists aorta, now perform step f3；Otherwise, it is determined that present image is absent from aorta, and perform step f4；

F3, choose in the second trace regions all connected regions and calculate each connected region cost function V (i, j):

V (i, j)=al (D (i, j))+bl (A (i, j))+cl (M (i, j)) (1)

Wherein, $l\left(D\right(i,j\left)\right)=\frac{D(i,j)}{\underset{m=1}{\overset{n}{\mathrm{\&Sigma;}}}D(i,m)}---\left(2\right)$

$D(i,j)=\sqrt{{(x_k^i-x_{k+1}^j)}^2+{(y_k^i-y_{k+1}^j)}^2}---\left(3\right)$

For the center-of-mass coordinate of kth frame image internal object,For the center-of-mass coordinate of jth connected region to be matched in kth+1 two field picture,For center-of-mass coordinate and the center-of-mass coordinate distance sum of all connected regions in kth+1 two field picture of kth frame image internal object, n is the number of all connected regions in kth+1 two field picture；(D (i, j)) represents jth connected region to be matched and kth frame image internal object similarity in position in kth+1 two field picture to l；

$l\left(A\right(i,j\left)\right)=\frac{\left|{\mathrm{Area}}_i-{\mathrm{Area}}_j\right|}{\underset{m=1}{\overset{n}{\mathrm{\&Sigma;}}}\left|{\mathrm{Area}}_i-{\mathrm{Area}}_m\right|}---\left(4\right)$

Area_iFor the area of kth frame image internal object, Area_jFor the area of jth connected region to be matched in kth+1 two field picture,For in area and kth+1 two field picture of kth frame image internal object the difference in areas of all connected regions and；(A (i, j)) represents jth connected region to be matched and kth frame image internal object similarity on area in kth+1 two field picture to l；

$l\left(M\right(i,j\left)\right)=\frac{M(i,j)}{\underset{m=1}{\overset{n}{\mathrm{\&Sigma;}}}M(i,m)}---\left(5\right)$

M (i, j)=L (i)/S (i)-L (j)/S (j) (6)

The major axis of L (i), S (i) respectively kth frame image internal object and minor axis length, the major axis of jth connected region to be matched and minor axis length in L (j), S (j) respectively kth+1 two field picture,Difference is asked to sue for peace respectively again for the ratio of semi-minor axis length of all connected regions in the ratio of semi-minor axis length of kth frame image internal object and kth+1 two field picture, (M (i, j)) represents jth connected region to be matched and kth frame image internal object similarity in eccentricity in kth+1 two field picture to l；

In formula (1), a, b and c are weight coefficient；

In kth+1 two field picture, (i, j) minimum connected region is the target tracked to cost function V, and the target tracked in this step is aorta；

F4, choosing the window corresponding with coronary artery previous frame image in present image as the first trace regions, the area of this first trace regions is area coronarius slightly larger than in previous frame image；

If coronary artery does not divide in present image, then coronary artery is tracked in the first trace regions of present image according to step f3；

If coronary artery there occurs division in present image, then find out the connected region that some barycenter are nearest from coronary artery barycenter in previous frame image, and the area of these connected regions and the difference of area coronarius is less than the 3rd area setting value with previous frame image, these connected regions found are the coronary artery tracked.

5. coronary artery CT angiographic image dividing method according to claim 4, is characterized in that, in formula (1), a is 0.5, b be 0.3, c is 0.2.

6. coronary artery CT angiographic image dividing method according to claim 1, is characterized in that, in step g, its system mode of Kalman filter is four dimensional vectors X (k)=(x_k,y_k,xv_k,yv_k)^T, externally input is U (k)=[xa_k,ya_k]^T；

Wherein: x_k,y_kThe respectively transverse and longitudinal coordinate of target centroid, xv_k,yv_kRepresent target speed in X-axis and Y direction, xa respectively_k,ya_kRepresent target acceleration in X-axis and Y direction respectively；The interval of adjacent two two field pictures is T and target does uniformly accelerated motion；

State-transition matrix A (k), systematic parameter B (k), observing matrix H (k) are respectively as follows:

$A\left(k\right)=\left[\begin{array}{cccc}1& T& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& T\\ 0& 0& 0& 1\end{array}\right]$ $B\left(k\right)=\left[\begin{array}{cc}0.5T^2& 0\\ T& 0\\ 0& 0.5T^2\\ 0& T\end{array}\right]$ $H\left(k\right)=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\end{array}\right]$

7. coronary artery CT angiographic image dividing method according to claim 6, is characterized in that, covariance matrix is respectively as follows:

$Q=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]$ $R=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]$

8. coronary artery CT angiographic image dividing method according to claim 1, is characterized in that, adopts maximum variance between clusters that filtered image is carried out Threshold segmentation, specifically according to formula below, filtered image is carried out Threshold segmentation in step a:

$f_t(x,y)=\left\{\begin{array}{c}1,f(x,y)\&GreaterEqual;\mathrm{\&lambda;}t\\ 0,f(x,y)<\mathrm{\&lambda;}t\end{array}\right\}---\left(7\right)$

In formula (7), t is segmentation threshold, and λ is for adjusting parameter.