CN101301207A - Vascular angiography three-dimensional rebuilding method under dynamic model direction - Google Patents

Abstract

The present invention provides a dynamic model directed angiography three-dimensional reconstruction method, pertaining to an intersectional field of digital image processing and medicine imaging, for satisfying special requests of auxiliary detection and surgery guidance for cardiovascular disease in clinical medicine. The invention includes steps of angiography image preprocessing, vascellum segmentation, vascellum skeleton and radii extraction, model directed vascellum base element recognition, vascellum matching and vascellum three-dimensional reconstruction. The invention also provides a cardiovascular dynamic model construction method including steps of cardiovascular slice data extraction, cordis static and dynamic model building, and static and dynamic model building. According to the invention, good angiography three-dimensional reconstruction result may be obtained, which effectively assists detection and surgery guidance for cardiovascular disease, thereby satisfying clinical requests.

Description

Angiography three-dimensional rebuilding method under dynamic model instructs

Technical field

The invention belongs to the crossing domain of Digital Image Processing and medical imaging, be specifically related to the angiography three-dimensional rebuilding method under a kind of dynamic model guidance.This method can solve by the various visual angles angiogram carries out reliably a three-dimensional reconstruction difficult problem automatically, satisfies the application requirements of clinical medicine cardiovascular disease auxiliary detection and surgical navigational.

Background technology

The vascular tree three-dimensional reconstruction is a process of recovering the blood vessel three-D space structure by corresponding image information in the X ray two-dimensional projection image of different visual angles.It and general visible 3-dimensional reconstruction have a great difference.X-ray imaging is the image that forms on fluorescent screen after the decay of X ray through different tissues in the human body of projection human body, the value of each pixel be by on the X ray path the stack of decay in a organized way, and noise is very strong, the background complexity, it is very big to obtain three-dimensional cardiovascular tree difficulty from x-ray imaging figure reconstruction.

At present hospital mostly usefulness be that x-ray imaging is done to patient by X ray single armed radiography system, obtain a radiography graphic sequence by rotation radiography arm corresponding to different radiography angles.The single armed radiography can easily carry out the radiography of different angles to patient, be exactly the radiography figure of our different visual angles that can't obtain synchronization but a shortcoming is arranged concerning rebuilding, and this brings very big difficulty to reconstruction.

Rebuild the true three-dimension space structure of vascular tree, need obtain the projection information of at least two different angles of vascular tree.Traditional method at first extracts the skeleton of blood vessel, then by different visual angles space constraint relation, the blood vessel pixel of different visual angles projected image is correctly mated and rebuilds, little when the blood vessel geometric deformation, when geometrical relationship is tangible, could recover the three-D space structure of blood vessel substantially.

Generally form by following 3 committed steps: the extraction of (1) vascular skeleton based on the three-dimensional reconstruction technology of the vascular tree of two width of cloth monoplane contrastographic pictures; (2) identification of characteristic point and coupling; (3) match of the estimation of blood vessel spatial point and reconstruction and blood vessel primitive.

Vascular tree can be regarded the crooked tubular system that extends in the space on the whole as, and its skeleton is the continuous space curve with tree, has reflected the configuration feature of vascular tree.Detect by the blood vessel profile, realize that it is the problem of a difficulty that vascular skeleton and radius extract accurately automatically.At present, the more feasible method in this aspect has ditetragon range searching (the double square box region of search) method of people's propositions such as Hoffman, be used for from the motion tracking blood vessel, obtain the parameters such as local size, position, axis of blood vessel simultaneously.People such as Coatrieux and Collorec is used for the extraction on vessels axis and border with the vector tracking algorithm, and it is few to have a manual intervention, can realize that blood vessel layering from coarse to fine is gone forward one by one to detect and save time etc. advantage.But it is high that these methods require original image quality, because x light contrast imaging quality is limited at present, no matter is definition, or resolution, utilizes existing image pre-processing method it to be strengthened poor effect.In addition, wanting blood vessel detected accurately with mating also needs a large amount of anatomical knowledge, will increase the complex nature of the problem and cost if introduce specialist system.Adopt the method for manual detection, man-machine interaction to carry out the extraction of vascular skeleton in clinical mostly.

For identification of the structure of cardiovascular double vision angle contrastographic picture and matching problem, have following difficult point: the blood vessel topological structure in the first, two contrast imaging face often differs bigger; The second, the existence that the projection of projected image medium vessels intersects is having a strong impact on the identification and the description of blood vessel topological structure.

Traditional method is to come the puncta vasculosa among two width of cloth radiography figure is mated with the method for outer polar curve constraint+linear programming, but we know, outer polar curve utilization be space geometry relation between the subpoint, its supposition be two subpoint correspondences be same three dimensions point.We can't obtain the radiography figure corresponding to synchronization in the single armed radiography, and because the motion puncta vasculosa of heart also can move, the match point that finds with outer polar curve constraint just has very mistake like this.Up to the present, also rest on for the most frequently used method of different visual angles blood vessel pixel coupling in the Reconstruction of vessel and manually choose the right method of matching characteristic point.This not only needs operator to have rich experiences, but also can take a long time.It is right to want effective recognition to go out the characteristic point of the coupling in the contrastographic picture of double vision angle, and set up fast and the effectively automatic algorithms of different visual angles blood vessel coupling, need take all factors into consideration the diameter information of seriality, blood vessel of vessel directions vector and the constraint of blood vessel topological structure or the like factor.

In view of the problem that exists in traditional blood vessel structure identification and the matching process, external also once the someone proposed to instruct the method for radiography figure coupling with the vascular tree model, wherein that most of usefulness is a general coronary artery MODEL C oronix by the special coronary artery modelling of 37 individualities according to the Dodge proposition.Because it is static being used for instructing the vascular tree model of radiography figure coupling, Many researchers both domestic and external mostly is in the reference diagram of the radiography figure of heart movement synchronization in the cycle as reconstruction by selection in sequence image, such as diastasis.But have two problems like this: how whether (1) be in the cardiac cycle the same moment to two width of cloth radiography figure is adjudicated.When usually judging the outside diastole of blood vessel to maximum (at this moment blood vessel wrapping heart size maximum) is diastasis, vasoconstriction is an end-systole to minimum (at this moment blood vessel wrapping heart size minimum) time, but with human eye judge still exist uncertain; (2) can slattern like this at radiography pictures, may lose a lot of Useful Informations except these selected constantly other moment.And do not consider like this effect of radiography figure yet, such as contrast agent due in, blood vessel block, noise etc., be difficult to guarantee that the angiogram of selected these different directions constantly is suitable for doing three-dimensional reconstruction.

After the correct corresponding relation that finds different views medium vessels pixel, just can recover the 3 d space coordinate of blood vessel corresponding pixel points by space geometry relations act or homogeneous coordinate transformation method.Recover the space coordinates of blood vessel pixel from limited X-ray projection view, this method can only be that a kind of similarity is rebuild, thereby often make up energy function with its Euclidean distance between corresponding image planes medium vessels structure by the projection of reconstruction blood vessel structure in all views, pass through optimization method, find the value of an optimum, be reconstructed results.

When carrying out blood vessel 3 D reconstructing, the general three-dimensional coordinate that does not need to reconstruct each point on the blood vessel, can obtain some sample points to sampling on each vessel segment among the radiography figure earlier, then these sample points are carried out three-dimensional reconstruction and obtain some three dimensions points, the method according to interpolation or curve fitting obtains whole section blood vessel from these three dimensions points structures again.The method of curve fitting is varied, can make up suitable cubic fitting model according to the specific form of blood vessel, describes the spatial shape of vascular skeleton more accurately.In the traditional method, curve can by comparing the fitting degree of curve in the actual view, be revised space vascular tree model with its back projection to image planes after spatial fit again, but this process need manual intervention, and accuracy is not high.Sophisticated now algorithm has B batten, Snake model etc., can express the vascular skeleton information between the different blood vessel primitive more accurately.For spline-fit, be actually on the basis of the local flatness of supposing curve and between sample point, carry out interpolation, match point is many more, and precision is high more; And, need structure suitable internal force and forcing function for the Snake model, and the effect of external force is to make matched curve approach projection as far as possible, the effect of internal force then is to make curve level and smooth as far as possible.Comparatively speaking, the Snake model can obtain higher precision, but internal force and forcing function structure are relatively difficult; The method of spline-fit is fairly simple, needs bigger amount of calculation but improve precision.

Summary of the invention

The angiography three-dimensional rebuilding method that provides a kind of dynamic model to instruct of the present invention, purpose is by a heart and cardiovascular four-dimensional dynamic model, instructs the blood vessel primitive projection coupling among the radiography figure of double vision angle, carries out three-dimensional reconstruction automatization, robust.

A kind of angiography three-dimensional rebuilding method provided by the invention based on dynamic model, its step comprises:

(1) chooses the X-ray coronarogram of two width of cloth different visual angles, they are carried out pretreatment, at first adopt the enhanced method of frequency domain to strengthen blood vessel to each radiography figure image sequence, adopt morphology Bottom-Hat conversion further to emphasize blood vessel then, eliminate interference noise, keep and outstanding required image information;

(2) pretreated contrastographic picture is carried out the segmented extraction angiosomes, the refinement angiosomes obtains vascular skeleton, select for use eight connection chain codes to carry out the cardiovascular skeleton and follow the tracks of, extract vessel radius, adopt the topological structure of binary tree structure storage vascular tree according to following process;

(2.1) the blood vessel axis in angiogram extracts in the image, begin the blood vessel axis is carried out binary tree search from the vascular tree root node, and the blood vessel pixel of searching for made marks, show that when searching the blood vessel pixel that labelling crosses circulation appears in search, judging has blood vessel to intersect, and stop to search for forward, forward step (2.2) to; The appearance if search procedure does not circulate forwards step (2.4) to;

(2.2) circulation time occurs when binary tree search,, extract the blood vessel that identifies intersection in the image at the blood vessel axis of angiogram according to the topology information and the geological information of vascular tree;

(2.3) get back to step (2.1);

(2.4) search finishes;

(3) respectively the blood vessel that the blood vessel axis of the angiogram of two width of cloth different visual angles extracts in the image is carried out labelling by following process:

(3.1) from different directions the vascular tree model of choosing is carried out the model projection figure that projection obtains different visual angles by the x-ray imaging mode;

(3.2) ask the blood vessel axis of each projection and angiogram to extract the similarity of image T by following process, find the projection P that wherein similarity is the highest;

(3.2.1) the blood vessel topological structure among the blood vessel axis extraction image T of vascular tree model projection figure and angiogram is kept at respectively among binary tree T1 and the T2;

(3.2.2) from binary tree T2, extract the subtree T3 that mates with the T1 topological structure; T1 and T3 are mated, utilize following formula to calculate the cost function of following T1 and T3:

f＝a ₁·f ₁(angle)+a ₂·f ₂(curve)+a ₃·f ₃(length)

F wherein ₁(angle) represent the similarity of blood vessel topological structure between the two, f ₂(curve) represent the similarity of blood vessel structure curve shape between the two, f ₃(length) similarity of two sections length of vessel of expression is expressed as the function of length of vessel length; a ₁, a ₂, a ₃Be the weight coefficient of these three kinds of similaritys, a ₁+ a ₂+ a ₃=1, the blood vessel of getting cost function minimum wherein is to being the coupling blood vessel;

(3.2.3) get wherein that the subtree of cost function minimum is the matched children M of T1, with the inverse of this minimum cost function as the similarity between model projection figure and the angiogram;

(3.3) obtain extracting the highest matched children M of projection P in T of similarity of image T by the above-mentioned step of asking similarity with the blood vessel axis of angiogram; Name with the projection of P medium vessels primitive is carried out labelling to the blood vessel primitive projection among the subtree M, as two radiography figure that choose not during the synchronization in cardiac cycle, the vascular tree model with these two moment carries out labelling to the blood vessel primitive projection of corresponding radiography figure respectively;

(4) carry out the blood vessel coupling by following process:

(4.1) find the blood vessel primitive projection of same tag among two width of cloth radiography figure as the projection of coupling blood vessel primitive;

(4.2) contrast two width of cloth radiography figure find a pair of blood vessel primitive projection of the identical blood vessel primitive projection of topological structure as coupling according to being limited in of topological structure in the blood vessel primitive projection that does not have labelling among two width of cloth radiography figure; Its process is:

(4.2.1) establishing this two width of cloth radiography figure medium vessels topological structure exists respectively among left_T and the right_T;

(4.2.2) find and do not have the blood vessel of labelling UnName1 among the left_T, whether the search of same position place also has unlabelled blood vessel structure UnName2 in right_T, may there be corresponding blood vessel structure in explanation if just have, and does not judge then that UnName1 is pseudo-blood vessel;

(4.2.3) if UnName2 exists, whether correspondingly one by one judge between UnName1 and the UnName2; If corresponding one by one, then directly UnName1 and UnName2 are ordered with same tag; If not corresponding one by one, then need to mate between any two, the blood vessel that finds the cost function minimum is to being the coupling blood vessel, and makes same tag;

(4.2.4) repeating step (4.2.1)-(4.2.3) is up to the whole labellings of the blood vessel in the vascular tree or till being judged as pseudo-blood vessel;

(5) instruct the blood vessel axial point coordinate transform of will extract among two width of cloth radiography figure to synchronization by dynamic model; Blood vessel primitive projection to topological structure coupling among two width of cloth radiography figure is further mated each last pixel, obtains the matched pixel point, carries out three-dimensional reconstruction.

The inventive method is considered the Topology Similarity of human body coronary artery system structure, respectively the blood vessel primitive projection among the x-ray imaging figure of two different visual angles is discerned with dynamic cardiovascular tree-model, thereby instruct the coupling of two width of cloth radiography figure medium vessels primitive projections, improve the reliability and the accuracy of blood vessel structure characteristic matching, increase the precision of blood vessel 3 D reconstructing.Particularly, the present invention has the technique effect of following three aspects:

(1) uses dynamic model to instruct, strengthen the robustness of blood vessel structure coupling under the different projection angles, compensated the influence of heart movement simultaneously.In three-dimensional reconstruction algorithm in the past, mostly the method for the outer polar curve of utilization is carried out the coupling of blood vessel, but the mistake matching rate of this method is very high, and the radiography figure that obtains when different angles is not when synchronization obtains especially.In order to overcome the deficiency of outer polar curve matching process, someone has proposed to instruct with the vascular tree model method of radiography figure coupling, though can partly improve the precision of blood vessel coupling, but still there is a problem in it: promptly because the vascular tree model that is used for instructing radiography figure to mate is static, coming that with it difference angiogram is constantly mated the result who obtains often has very big error, particularly vasomotion and causes blood vessel to occur overlapping and when intersecting.Here we use dynamic vascular tree model that the blood vessel coupling is instructed, respectively the blood vessel among the radiography figure is carried out labelling by the model projection figure in the difference moment and the topological structure coupling of choosing between the angiogram that is used for rebuilding, overcome the error that vasomotion brings.

On the other hand, because two width of cloth radiography figure that rebuild can not rebuild them with the general three dimension reconstruction formula not at synchronization.Therefore we have utilized the blood vessel dynamic model that radiography figure medium vessels is carried out motion compensation, make vessel position in two width of cloth contrastographic pictures corresponding to the same moment as far as possible, and then utilize the blood vessel 3 D reconstructing formula that it is rebuild, can reduce the error that vasomotion brings to reconstruction like this.

(2) set up multiple dimensioned vascular tree model, instruct the blood vessel coupling better.Also useful abroad model instructs blood vessel labelling and Matching Algorithm, but the model that they use all is a fixed model, does not have the variation of yardstick, can not satisfy the requirement of Different Individual for details.Here we replenish model in by the process that instructs radiography figure coupling in double vision angle with model simultaneously, the details blood vessel that does not have before can adding.

We propose the notion with the multiple dimensioned vascular tree model of vascular tree progression equivalence, and vascular tree is multiple dimensioned to be the progression of relative vascular tree.According to the anatomical knowledge of vascular tree, vascular tree can be divided into the blood vessel of a lot of different progression, and first order blood vessel is the blood vessel of root, be exactly to connect aortal coronary artery in the arteria coronaria system, its branch is exactly the secondary blood vessel, in like manner, three grades of blood vessels are branches of secondary blood vessel, successively down.When rebuilding, we can select to use the vascular tree model of different scale or different progression that radiography figure is mated guidance, thereby can more effectively instruct the registration and the reconstruction of blood vessel according to the practical situation of rebuilding requirement and radiography figure.Only require third level blood vessel such as rebuild, we just only need provide three grades of vascular tree models to instruct coupling.Cut apart perhaps that getable maximum detail can only reach three grades among the radiography figure of back, we are same only need to instruct coupling with three grades of vascular tree models.By setting up the robustness that multiple dimensioned vascular tree model has improved blood vessel coupling and three-dimensional reconstruction, improved the reconstruction precision.

(3) the blood vessel matching algorithm is not too dependent on the geological information of model, can well solve the problem of individual difference.In a lot of algorithms abroad, instruct the blood vessel coupling all to utilize a lot of information with model about model geometric, coordinate such as the projection medium vessels, blood vessel primitive length etc., and these can be subjected to the influence of individual differences such as blood vessel size, distortion and skew, and it is used in the robustness that can reduce algorithm in the algorithm greatly; By contrast, topological structure is just much stable than geological information, and it can not be subjected to the influence of picture size, has rotational invariance simultaneously, and not influenced by translation.Given this, our method has only been utilized the topological structure and the part geological information not too responsive with rotation to picture size of blood vessel, as blood vessels adjacent primitive length ratio, blood vessel primitive curvature (be defined as on the blood vessel primitive had a few average curvature) etc., instruct the coupling of blood vessel, thereby make this method can reduce sensitivity, make process of reconstruction become robust and more high automation degree more individual difference.

Description of drawings

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

Fig. 2 is the FB(flow block) of cardiovascular modeling part;

Fig. 3 is the FB(flow block) that x-ray imaging is rebuild;

Fig. 4 is a coordinate system of setting up heart and vascular tree model;

Fig. 5 (a) is the heart sectioning image;

Fig. 5 (b) is the image that extracts cardiac component in the heart sectioning image, and wherein regional A2 is a left ventricle, and regional A3 is a right ventricle, and regional A1 is a left atrium, and regional A4 is a right atrium;

Fig. 6 is a cardiac silhouette sampling point diagram;

Fig. 7 is the heart three-dimensional model diagram, and wherein outmost profile is a pericardium, and area B 1 is a left ventricle, and area B 2 is rights ventricle, and area B 4 is left atriums, and area B 3 is right atriums;

Fig. 8 (a)～8 (c) is the cardiod diagram that has the individual specificity by the position simulation that changes the control point;

Fig. 8 (a) is normal heart figure;

Fig. 8 (b) is heart differentially expanding figure;

Fig. 8 (c) is a contract drawing in the heart part;

Fig. 9 is that cardiac cycle is schemed each period

Figure 10 is left ventricular volume linear change figure;

Figure 11 (a)～11 (g) heart is at the cardiac cycle illustraton of model in each period;

Figure 11 (a) is that heart is in the isovolumic contraction period illustraton of model;

Figure 11 (b) is that heart is in the phase of maximum ejection illustraton of model;

Figure 11 (c) is that heart is in the slow ejection period illustraton of model;

Figure 11 (d) is appearance filling period illustratons of model such as heart is in;

Figure 11 (e) is that heart is in the phase of rapid filling illustraton of model;

Figure 11 (f) is that heart is in the slow filling period illustraton of model;

Figure 11 (g) is that heart is in the Atrial systole illustraton of model;

Figure 12 (a) is the 48th sectioning image in the heart section sequence;

Figure 12 (b) is the vascular cross-section image that (a) extracted;

Figure 12 (c) is the 50th sectioning image in the heart section sequence;

Figure 12 (d) is the vascular cross-section image that (c) extracted;

Figure 13 is the model sketch map through blood vessel primitive labelling;

Figure 14 (a)～14 (g) is that coronary artery is in the cardiac cycle illustraton of model example in each period;

Figure 14 (a) is that the left coronary artery blood vessel is in heart isovolumic contraction period illustraton of model example;

Figure 14 (b) is that the left coronary artery blood vessel is in heart phase of maximum ejection illustraton of model example;

Figure 14 (c) is that the left coronary artery blood vessel is in heart slow ejection period illustraton of model example;

Figure 14 (d) is that the left coronary artery blood vessel is in appearance filling period illustraton of model examples such as heart;

Figure 14 (e) is that the left coronary artery blood vessel is in heart phase of rapid filling illustraton of model example;

Figure 14 (f) is that the left coronary artery blood vessel is in heart slow filling period illustraton of model example;

Figure 14 (g) is that the left coronary artery blood vessel is in heart Atrial systole illustraton of model example;

Figure 15 is left coronary artery blood vessel and the variation tendency of background noise point gray scale in a radiography graphic sequence of choosing, wherein curve 1 denotation coordination is (252, the gray value of blood vessel pixel 134), curve 2 denotation coordinations are the gray value of the background dot of (362,95).As can be seen from the figure the grey scale change of puncta vasculosa is bigger than background noise point, blood vessel and background area can be separated according to these characteristics;

Figure 16 (a) is the left coronary artery angiographic image;

Figure 16 (b) is that the left coronary artery angiography strengthens image;

Figure 16 (c) is a left coronary artery angiography Bottom-Hat modified-image;

Figure 17 is that the circulation sketch map appears in the blood vessel search;

Figure 18 is a vessel projection intersection sketch map;

Figure 19 is relatively sketch maps of two width of cloth radiography figure medium vessels topological structure, wherein if vessel segment E1, E2 respectively and F1, the F2 correspondence judges that then the topological structure of vessel segment UE and UF is identical, otherwise different;

Figure 20 (a) instructs the ramose labelling sketch map of rebuilding of vascular tree model medium vessels;

Figure 20 (b) instructs the result schematic diagram that contrastographic picture medium vessels branch is marked with model;

G1 among Figure 20 (a), G2 ..., G7 respectively with Figure 20 (b) in H1, H2 ... the H7 labelling is identical;

Figure 21 is two unmarked blood vessels same topology location sketch maps in vascular tree, vessel segment I1 wherein, and I2, I3 are same tag, therefore judge unmarked vessel segment U1, the topological structure of U2 is identical;

Figure 22 is the end product that respectively blood vessel in the angiogram of two width of cloth different visual angles is carried out labelling and coupling, the vessel segment J1 among radiography Figure 22 (a) wherein, J2 ..., J13 respectively with 22 (b) in vessel segment K1, K2 ..., K13 mates corresponding;

Figure 23 (a)～Figure 23 (i) be by model instruct the coronary artery three-dimensional reconstruction process and with the comparison of other reconstructing blood vessel algorithm;

Figure 23 (a) is LCA LOOK LEFT (26.8 ,-a 27.2) contrastographic picture;

Figure 23 (b) is LCA LOOK RIGHT (50.8, a 30.2) contrastographic picture;

Figure 23 (c) is LCA LOOK LEFT blood vessel segmentation and the figure as a result that extracts axis;

Figure 23 (d) is LCA LOOK RIGHT blood vessel segmentation and the figure as a result that extracts axis;

Figure 23 (e) is the model projection figure that mates with LCA LOOK LEFT topological structure;

Figure 23 (f) is the model projection figure that mates with LCA LOOK RIGHT topological structure;

Figure 23 (g) be with our method to (a), (b) in radiography figure carry out the result of blood vessel 3 D reconstructing;

Figure 23 (h) is to use the method for conventional polar curve constraint to (a), and (b) middle radiography figure carries out blood vessel 3 D reconstructing result, reconstruction failure;

Figure 23 (i) be method by manually choosing match point to (a), (b) in radiography figure carry out the angiography three-dimensional reconstruction result;

Figure 24 (a) is the image that 23 (g) medium vessels three-dimensional reconstruction result is pressed radiography angle projection among Figure 23 (a);

Figure 24 (b) is the image that 23 (g) medium vessels three-dimensional reconstruction result is pressed radiography angle projection among Figure 23 (b);

Figure 24 (c) and (d) and Figure 23 (c) and (d) identical;

Figure 25 (a)～(c) is the multiple dimensioned model sketch map of blood vessel;

Figure 25 (a) is a left coronary artery one-level vascular tree model sketch map;

Figure 25 (b) is a left coronary artery secondary vascular tree model sketch map;

Figure 25 (c) is three grades of vascular tree models of left coronary artery sketch map;

Figure 26 (a) uses the inventive method to rebuild to the radiography figure in the accompanying drawing 23 to obtain the blood vessel three dimensional structure;

Figure 26 (b) is the result who Figure 26 (a) medium vessels axis is added radius;

Figure 26 (c) is the vascular tree model of setting up according to visual human's slice of data;

Figure 26 (d) carries out subsidiary details result afterwards with reconstructed results in (a) to model in (c), and wherein the vessel segment that identifies with circle is the blood vessel primitive that adds to thin yardstick in the model;

Figure 27 is the imaging sketch map of radiography system two different angles.

The specific embodiment

The present invention is further described below in conjunction with accompanying drawing and example:

As shown in Figure 1, the inventive method is to utilize heart and cardiovascular dynamic model to carry out blood vessel 3 D reconstructing from the X-ray contrast image of double vision angle, heart that is utilized and cardiovascular dynamic model can adopt original known heart and cardiovascular tree-model, also can adopt each one self-built heart and cardiovascular dynamic model, hereinafter provide a kind of method of setting up heart and cardiovascular dynamic model, with for referencial use.

(1) sets up the cardiovascular tree-model

Be illustrated in figure 2 as and set up cardiovascular tree-model flow chart.

When the three-dimensional static of setting up heart and blood vessel and dynamic model, to set up coordinate system: x, y axle and initial point as follows and overlap with x, y axle and the initial point of the slice map of the apex of the heart, the z axle is perpendicular to slice plane, the direction at the bottom of from the apex of the heart to the heart.(as shown in Figure 4)

(1.1) set up the heart dynamic model

The purpose of setting up the heart dynamic model is to set up the four-dimensional dynamic model of coronary artery in order to help, and its step comprises:

(1.1.1) extract heart section profile;

According to the picture of anatomical knowledge and the cross-section dissection of human body, in every width of cloth original image, tell pericardium, left ventricle, right ventricle, five parts of left atrium and right atrium.

The initial data picture comes from the U.S.'s visual people's plan (VHP:Visible Human Project) data set.The picture that we choose is the cardiac component of the anatomic image of Visible Man data centralization.Fig. 5 (a) is the original heart slice map that we obtain.

The structure more complicated of cardiac component, we cut apart by photoshop.Fig. 5 (b) is our segmentation result to the heart slice map among Fig. 5 (a).Outmost profile is a pericardium, and regional A2 is a left ventricle, and regional A3 is a right ventricle, and regional A1 is a left atrium, and regional A4 is a right atrium.

(1.1.2) above-mentioned heart section profile is sampled, obtain sampled point.

Because data can increase amount of calculation too much, and do not need too many data to carry out curve fitting under most of situation, therefore at first initial data is sampled, obtain some representational sampling points.Consider the characteristics of heart chamber and surface smoothing, all average sampling at heart section contour direction (laterally) and on perpendicular to slice direction (vertically) both direction.

Uniform sampling transversely is the center (centre of form of profile) of at first finding out profile, then from center emission n bar ray, wherein adjacent interradial angle is the 360/n degree, and the intersection point of choosing these rays and profile is a sampled point, so the number of sampled point also is n.Fig. 6 is 12 uniform sampling points that obtain when getting n=12.

Owing to vertically on (z axle) 138 width of cloth slice maps are arranged, can obtain some sampling pictures and carry out modeling by get a width of cloth slice map every several sections.

What should be noted that a bit is, sample point is few more, and the model that obtains is level and smooth more; Sample point is many more, and the details of reservation is many more, but amount of calculation is also big more.

(1.1.3) utilize above-mentioned sample point, carry out the B-spline surface match, set up heart three-dimensional static model according to following formula:

The B-spline surface fitting formula is $P(u,v)=\underset{i=1}{\overset{N_u}{\mathrm{\Σ}}}\underset{j=1}{\overset{N_v}{\mathrm{\Σ}}}{C}_{\mathrm{ij}}{B}_i\left(u\right)B_j\left(v\right),$ P is the point on the model, and u, v are corresponding this P laterally (U to) of order and vertical (V to) parameter, N _u, N _vBe respectively U in the B batten model, the control point number on the V direction.B _i(u), B _j(v) then be respectively at U, the B spline base function on the V direction, C _IjBe the control point of B-spline surface, constitute N _u* N _vThe control point grid, determined the shape of B-spline surface.

B-spline surface can be regarded as the curved surface after the B-spline curves by both direction carry out tensor product.Therefore the surface fitting process can be divided into following two steps:

1) u to the B-spline curves match;

The B-spline curves fitting formula is $P\left(u\right)=\underset{i=1}{\overset{N_u}{\mathrm{\Σ}}}{C}_i{B}_i\left(u\right),$ The fitting data here is the sample point on each heart section profile.These points are carried out B-spline curves match obtain the cutting into slices control point C ' of contour curve _Ij(j=1,2 ..., N _u).

2) v to the B-spline curves match.

The fitting data in this step is the control point C ' of the section contour curve that obtains of previous step _Ij(j=1,2 ..., N _v).They are carried out the control point C that the B-spline curves match obtains whole B-spline surface _Ij(i=1,2 ..., N _u, j=1,2 ..., N _v).

Behind the control point that obtains B-spline surface, substitution B-spline surface equation just obtains the parameter model of pericardium and each chamber.Specify a parameter value (u ₀, v ₀), just can calculate its some P (u corresponding in the space by the surface fitting formula ₀, v ₀) coordinate (x ₀, y ₀, z ₀).Can be by changing control point position change heart shape, simulation has the heart of individual character.

Fig. 7 lumps together the three-dimensional static model that constitutes heart with heart chamber and surface.Outmost profile is a pericardium, and area B 1 is a left ventricle, and area B 2 is rights ventricle, and area B 4 is left atriums, and area B 3 is right atriums.

Fig. 8 is the heart by the mimic individual specificity of having in position who changes the control point.

(1.1.4) set up the heart four-dimensional dynamic model: set up dynamic heart model by the motion model of setting up all chambers of the heart chamber and pericardium.

The motion of heart is very complicated, experimental result according to medical observation and biomedical engineering proves, ventricular wall motion by cardiac muscle interior to contractile motion, heart move horizontally and heart form in addition along three kinds of main motion modes of rotation of axle reverse in addition, local motion such as local stretching, extension.In different forms of motion, as the source of the power of heart transportation blood, the expansion motion has nearly accounted for 90%, so we mainly consider this motion mode of heart.We do following hypothesis earlier before modeling:

1. the synchronized movement of left atrium and ventricle;

2. pericardium is along with chamber moves;

3. ignore the influence of valvular motion to heart chamber;

4. the direction of cardiac muscle on the heart wall is taken as the direction of the middle cardiac muscle of heart wall, promptly along the heart wall circumferencial direction;

5. the lip-deep point of heart chamber and pericardium all moves along normal direction;

6. we are divided into 7 time periods with cardiac cycle, and heart chamber and pericardium all are approximately linear movement in each time period of reasonable assumption.

According to the physiological property of heart movement, we are divided into seven time periods to cardiac cycle, analyze chamber pressure and the volume-variation of each time period in each cardiac cycle.Fig. 9 is each pressure variation in period and ventricular volume variation diagram of cardiac cycle, and wherein C1 represents atrial systole, and C2 represents isovolumic contraction period, C3 represents fast rapid fire blood, and C4 represents slow ejection period, and C5 represents protodiastole, C6 represents isovolumic relaxation phase, C7 represents quick filling phase, and C8 represents slow filling phase, and C9 represents aortic pressure, C10 represents ventricular volume, C11 represents atrial pressure, and C12 represents ventricular pressure, and C13 represents electrocardiogram.Figure 10 is that we are with the linearizing result of the volume-variation of left ventricle.

Based on the hypothesis of front and approximate, the general steps of dynamic heart modeling can be divided into for three steps: at first set up the motion model of heart chamber, then set up the motion model of pericardium, at last they are combined the four-dimensional dynamic model that constitutes whole heart.Owing in the process of setting up the chamber dynamic model, considered the effect between the adjacent chamber, situation about clashing therefore can not occur combining.The four-dimensional dynamic model of whole heart can be shown with VTK at last:

(1.1.4.1) make up the chamber dynamic model: the main thought of chamber dynamic modeling step is exactly to release the deformation of chamber from the variation of chamber volume, promptly $V\left(t\right)\⇒r\left(t\right),$ Here V (t) represent chamber volume over time, on behalf of the deformation of chamber, r (t) change the displacement of the point on the curved surface (be over time) in time.At first determine that according to the form of the static models of heart it is at a heart movement residing position t in the cycle ₀, set up the dynamic model of four chambers of heart then respectively according to following step:

1) according to the kinetic characteristic of heart, the heart one-period is divided into 7 different periods, chooses the time point of the starting and ending in these seven periods and put t as 7 reference times _i(i=1 ..., 7), t wherein ₁Be the next time of t0.

2) at time t _iThe time, obtain with respect to time t according to curve V (t) _I-1Change in volume Δ V.

3) the surface area S on calculating left ventricle surface (not comprising the chamber junction).So-called chamber junction is exactly atrioventricular valves (A V valves) and these places of interventricular septum.

4) calculate from time t _I-1To t _iThe displacement along this normal direction that middle each point produces is d=Δ V/S.

5) each sample point P to (not comprising the chamber junction) on the curved surface calculates the surface normal direction at a P place

Order $\stackrel{\→}{{\mathrm{OP}}^{\′}}=\stackrel{\→}{\mathrm{OP}}-d\·\stackrel{\→}{n},$ Wherein O is a zero, and some P ' is the reposition of the sample point after the deformation, in left ventricle,

$d=\left\{\begin{array}{cc}210t/S\left(t\right),& t\∈\left[\mathrm{0,0.1}s\right)\\ 0,& t\∈\left[\mathrm{0.1,0.15}\right)\\ 470(0.15-t)/S\left(t\right),& t\∈\left[\mathrm{0.15,0.25}\right)\\ \frac{460}{3}(0.25-t)/S\left(t\right),& t\∈\left[\mathrm{0.25,0.40}\right)\\ 0,& t\∈\left[\mathrm{0.40,0.47}\right)\\ \frac{4600}{11}(t-0.47)/S\left(t\right),& t\∈\left[\mathrm{0.47,0.58}\right)\\ \frac{150}{11}(t-0.58)/S\left(t\right)& t\∈\left[\mathrm{0.58,0.80}\right)\end{array}\right.,,$

S (t) is at t heart chamber surface area constantly.

6) modeling procedure with the front comes new sample point reconstruct heart model.Store new heart model and time corresponding t _i

7) move on to next time point t _I+1, the 2nd to the 6th step above repeating finishes up to a cardiac cycle.

8) will be corresponding to putting t access time _iThe heart chamber model couple together the dynamic model that constitutes chamber, make the point on the heart chamber between adjacent time point be approximately linear movement;

(1.1.4.2) the pericardium dynamic model makes up: the modeling of pericardium is based on the dynamic model of chamber and the contracting model of cardiac muscle, and step is listed below.Based on another hypothesis, promptly myocardial volume is constant during the course.Make that heart is V shrinking forward and backward volume ₁And V ₂, d ₁And d ₂Be the thickness of cardiac muscle between pericardium and chamber, L ₁And L ₂Be the length of cardiac muscle fiber.According to hypothesis before, we can obtain following formula:

$\begin{array}{c}{L}_1{\mathrm{<figure-callout\; id="1"\; label="d"\; filenames="A20081004785300444.png"\; state="\{\{state\}\}">d</figure-callout>}}_1=L_2{d}_2\\ \frac{V_1}{V_2}=\frac{{\mathrm{<figure-callout\; id="1"\; label="L"\; filenames="A20081004785300444.png"\; state="\{\{state\}\}">L</figure-callout>}}_1^3}{L_2^3}\end{array}\&DoubleRightArrow;\frac{d_1}{d_2}={\left(\frac{V_1}{V_2}\right)}^{1/6}$

Concrete steps are as follows:

1) from time t ₀Beginning.

2) to each sample point P of pericardium, calculate normal in this point Find normal

Cross point Q with chamber.

3) at time t=t _iThe time, obtain the displacement vector of a Q according to chamber dynamic modeling algorithm

4) according to the thickness d of the myocardium THICKNESS CALCULATION heart wall between pericardium and the chamber=| PQ| obtains new heart wall thickness then $d^{\&prime;}=d\&CenterDot;{\left(\frac{V_i}{V_{i+1}}\right)}^{1/6}.$

5) order $\stackrel{\&RightArrow;}{{\mathrm{OP}}^{\&prime;}}=\stackrel{\&RightArrow;}{\mathrm{OP}}+(\stackrel{\&RightArrow;}{s}\&CenterDot;\stackrel{\&RightArrow;}{v}+d-d^{\&prime;})\&CenterDot;\stackrel{\&RightArrow;}{n}$ Be the sample point reposition vector after the deformation.

6) rebuild the pericardium model according to new sample point, preserve model and time corresponding t thereof _i

7) move to next time point t _I+1, the step 2 of repetition modeling) finish up to a cardiac cycle to step 6).

8) will be corresponding to putting t access time _iThe pericardium model couple together the dynamic model that constitutes pericardium, make the point on the pericardium between adjacent time point be approximately linear movement;

(1.1.4.3) foundation of heart dynamic model: the dynamic model of heart is to be made of one group of different threedimensional model constantly in the cardiac cycle.We have been divided into 7 sections with cardiac cycle in the superincumbent hypothesis, and the heart in supposing every section is approximately linear movement.Because our initial data is the anatomical data of human body, the volume size of calculating in physiological property after stopping to beat about heart according to related data and the heart static models infers that the original static model that the front is set up is in Atrial systole.Therefore we can set up the transient state threedimensional model of 8 time points earlier, describe transient process between them with linear change then, have so just constituted the motion in whole cycle of heart.These eight time points corresponding respectively 7 sectional first and last time points of the initial time of dynamic model (being original static model residing moment in cardiac cycle) and cardiac cycle.Here should be noted that any be exactly interval between the adjacent time point be not to equate, but according to actual cardiac cycle the time come mutually to determine.An intervalometer is set, responds different heart model data constantly, can in VTK, show dynamic heart.

To be us comprise the heart dynamic model VTK result displayed in 7 stages with the cycle of setting up to Figure 11, and every width of cloth subgraph respectively is an example in each stage.

(1.2) the four-dimensional dynamic model of coronary artery is set up

Setting up the four-dimensional dynamic model of coronary artery is in order to instruct the coronary artery three-dimensional reconstruction from several different projection angle x-ray imaging figure, to the steps include:

(1.2.1) Vascular Slice profile extraction step:, in every width of cloth original image, tell blood vessel according to the picture of anatomical knowledge and the cross-section dissection of human body.Figure 12 is original heart sectioning image and vessel extraction thereof a comparison diagram as a result.

(1.2.2) extract blood vessel axial point and radius step: set up axial point and two information of radius that the coronary artery model will extract blood vessel.Axial point is got the centre of form that is decided to be the Vascular Slice profile, supposes that simultaneously blood vessel is circular, and the radius of blood vessel is the inscribed circle radius of profile.

(1.2.3) rebuild blood vessel three-dimensional framework step: each ramose axial point of coronary artery that previous step is extracted is carried out the three-dimensional framework model that match obtains blood vessel with B-spline curves.

(1.2.4) construction step of the blood vessel three-dimensional model of band radius: on described blood vessel three-dimensional framework model, use the GC model to add radius information, set up static coronary artery model.

Generalized cylinder represents and makes general cone represent that it is that a kind of pushing away swept model.Form by axis, cross section and scanning rule.Axis is a three-dimensional space curve, and the cross section is a plane domain, and it scans at a certain angle along axis, and the shape in cross section can change in scanning process, thereby can represent different bodies.Scanning rule has defined the geometrical relationship between cross section and the axis.Object surfaces is exactly the union of these cross section closed curves, and pushing away the volume of sweeping formation is exactly generalized cylinder.Many objects with axial symmetry can be well described in this expression, and similar object has similar axle and cross section.

(A, E α) define blood vessel to generalized cylinder (GC) with a ternary formula.(Z) (s) is axis to A for X, Y, and it is a space curve by the parameter s definition; α is the inclination angle, and it is by plane, cross section place and locates the angle of axis tangent line at A (s); E=(t s) is a plane curve, it with in the planar parameter-definition of A (s) cross section, can be divided into two functions: E (s, t)=and r (s) C (t), wherein profile function C (t) has described the shape of cross section, and function of radius r (s) has described its size.Here we suppose that the cross section of blood vessel is circular C (t)=(cos2 π t, sin2 π t), and perpendicular to the tangential direction of axis (being that inclined angle alpha is 90 degree).

Utilize this model that blood vessel is carried out said three-dimensional body and describe, think that vascular tree is the shape of a series of cross sections and all constantly set of tubulose (cylinder) structure of variation of direction of size and axis.

(1.2.5) at last to each blood vessel primitive in the model according to Dodge et al. (Dodge JT, Brown BG, Bolson EL, Dodge HT: " Intrathoracic spatial location of specified coronary segments on thenormal human heart. " Circulation 78:1167-1180,1988.) She Ji nomenclature is named.Figure 13 is the result that we name model medium vessels primitive.

(1.2.6) set up the four-dimensional dynamic model step of arteria coronaria: motion coronarius is kinetic by heart, and it is the same with heart carries out cycle movement.Suppose that blood vessel keeps tubulose, radius does not change in motor process.The motion modeling of blood vessel was divided into for two steps.

1) sets up the motion model of the axis of blood vessel

Some points above the sampling blood vessel axis make them according to the motion of pericardium laws of motion, gain knowledge according to heart physiological, and cardiac cycle is divided into 7 stages, supposes that reasonably these sample points are approximately linear movement respectively in these 7 stages.To the blood vessel sample point by setting up the motion model in this 7 stages respectively, obtain they whole cardiac cycle motion model.Then these sample points are carried out the B-spline curves match in the position of synchronization and obtain corresponding to shaft model in each blood vessel constantly, in the one-period in all blood vessels constantly shaft model just formed the periodic movement model of blood vessel axis.Figure 14 is the blood vessel axis motion model in one group of seven moment of foundation.

2) on the motion model of axis, add radius information

Here on dynamic vascular skeleton model, add radius information by the GC model equally, obtain dynamic coronary artery model.

(2) at least two width of cloth X-ray contrast figure from different perspectives carry out blood vessel 3 D reconstructing (being illustrated in figure 3 as the blood vessel 3 D reconstructing flow chart):

(2.1) angiogram pretreatment: in order to eliminate interference noise, keep and outstanding required image information, finally from contrastographic picture, extract blood vessel structure exactly, at first adopt the enhanced method of frequency domain to strengthen blood vessel, adopt morphology Bottom-Hat conversion further to emphasize blood vessel then image sequence.

(2.1.1) frequency domain strengthens

Because heart has the period of motion, what obtained by the radiography system imaging is that a coronarogram is as sequence.Study this image sequence, we find that common rib, spinal cord, lung and some interior tissues can be regarded as the signal of constant or slow variation, yet coronary artery presents the kinestate that similar frequencies is arranged with heart along with the contraction and the expansion of heart.Figure 15 is the variation tendency of gray scale in the radiography graphic sequence of left coronary artery blood vessel and background noise point, and what wherein choose is a patient's of Beijing ZhaoYang Hospital radiography graphic sequence.Wherein curve 1 denotation coordination is the gray value of the blood vessel pixel of (252,134), and curve 2 denotation coordinations are the gray value of the background dot of (362,95).As can be seen from the figure the grey scale change of puncta vasculosa is bigger than background noise point, can remove background noise by high-pass filtering according to these characteristics.

If a frame number is the image sequence f of T ₀(t), the single image size is M * N for x, y, 0≤x＜M, and 0≤y＜N, 0≤t＜T, x, y, t is integer.Adopting following method to carry out frequency domain to image sequence strengthens:

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 ₀(x, y, k), and 0≤k＜T, the discrete Fourier transform (DFT) formula is as follows: wherein j is an imaginary unit

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

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

2) to frequency-region signal F ₀(x, y k) adopt following formula to carry out M * N high-pass filtering along the k axle,

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

The high-pass filtering function be H (x, y, k)=(1-e ^{-β t}), wherein β is a constant given in advance, in order to suppress the background noise of some frequency.In theory, the span of β is [0, ∞], when β → 0, and (1-e ^{-β t}) → 0, (k) → 0, all signals comprise that coronary artery all is removed to F for x, y, make that the inhibition noise is excessive; As β → 0, (1-e ^{-β t}) → 1, and F (x, y, k) → F ₀(k) all signals comprise that coronary artery all is retained for x, y, do not play the effect of removing noise.Suitable β can remove noise to the full extent and keep the arteria coronaria signal.Through repeatedly experiment, the preferable range of β value is 1/8～1/5.

3) utilize the following formula inverse discrete fourier transform to act on filtered frequency-region signal, obtain new image sequence:

$f(x,y,t)=\underset{k=0}{\overset{T-1}{\mathrm{\&Sigma;}}}\stackrel{\&OverBar;}{F}(x,y,k)e^{j2\mathrm{\&pi;k}/T}$

(2.1.2) multiple dimensioned morphology Bottom-Hat conversion

We adopt morphology Bottom-Hat conversion further to strengthen blood vessel.The Bottom-Hat conversion is the morphological operator that the dark structure of image is extracted in a kind of common being used to, and is defined as g=f-(fb), and wherein, f is an input picture, and b is the structural element function, $f\&CenterDot;b=(f\&CirclePlus;b)\mathrm{\&Theta;b}$ Represent closed operation, by expansive working

Operation is formed with corrosion Θ.Closed operation can be removed structures of interest darker in the image of being represented by structural element b, and keeps than bright pixel.In contrastographic picture, the gray value of target object (blood vessel) is low than background, and therefore, the result of fb can be regarded as background, can obtain vascular site by former figure subtracting background image.

In realization, we select two different-diameter yardstick D ₁, D ₂(D ₁＞D ₂) the disc structure operator, than major diameter D ₁Need be slightly larger than the wideest coronary artery, than minor diameter D ₂Be slightly larger than the narrowest coronary artery and choose the not high small artery of contrast to concentrate.(x, y) carrying out yardstick respectively is D to strengthening image f for we ₁And D ₂The Bottom-Hat conversion, obtain g ₁(x, y) and g ₂(x, y).

Through the image that obtains after the above-mentioned conversion is the grayscale image of outstanding blood vessel, but still contains background noise, needs respectively to g ₁(x, y) and g ₂(x y) carries out following steps to obtain bianry image.

1) relative threshold relatively.The Bottom-Hat conversion has been strengthened the dark structures of interest in the image and has been weakened other parts, make the vasculature part pixel value of radiography figure after the Bottom-Hat conversion be higher than the vasculature part pixel value among the radiography figure before the conversion, the background parts pixel value is lower than the background parts pixel value before the conversion.By choosing appropriate threshold T ₁, with pixel value and the T of figure after the pixel value of figure before the Bottom-Hat conversion and the Bottom-Hat conversion ₁Compare, can draw which point is target, and which point is a background.Because widely different between the image, distributing automatically according to the gradation of image value, calculated threshold has the suitability widely.We call this method " relative threshold is relatively ", and it can access the binary map that blood vessel and background segment are opened:

At first, adopt following formula relatively g (x, y) and f (x, y) obtain differential chart d (x, y),

$d(x,y)=\left\{\begin{array}{cc}g(x,y)-f(x,y)& \mathrm{if}(g(x,y)-f(x,y)>0)\\ 0& \mathrm{otherwise}\end{array}\right.,$

Then, (x, y) self adaptation is calculated threshold value T according to d ₁, T ₁For d (x, the y) meansigma methods of the non-zero pixels value among the figure,

${\mathrm{<figure-callout\; id="1"\; label="T"\; filenames="A20081004785300444.png"\; state="\{\{state\}\}">T</figure-callout>}}_1=\frac{\underset{d(x,y)>0}{\mathrm{\&Sigma;}}d(x,y)}{N_d},$

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

At last, if d (x, y) 〉=T ₁, then this is the blood vessel picture element, if d (x, y)＜T ₁, then this is a background dot.With g ₁(x, y) and g ₂(x y) carries out above-mentioned processing respectively and obtains binary map B ₁(x, y) and B ₂(x, y).

2) integrate binary map.Stack binary map B ₁(x, y) and B ₂(x, y), purpose is to strengthen terminal blood vessel and little blood vessel, removes some areas simultaneously and is lower than threshold value T ₂The block distortion zone.At this moment, we have just obtained the whole extraction result of coronary artery.

Figure 16 strengthens primitive vessel radiography figure and Bottom-Hat conversion result afterwards through frequency domain.Figure 16 (a) is the 22nd frame of original radiography sequence chart, and Figure 16 (b) is the result after the process frequency domain strengthens, and Figure 16 (c) carries out Bottom-Hat variation result afterwards with strengthening image.

(2.2) extract vascular skeleton and radius: pretreated image is carried out the segmented extraction angiosomes, adopt morphology methods refinement angiosomes to obtain vascular skeleton.Select for use eight connection chain codes to carry out the tracking of cardiovascular skeleton, extract vessel radius, with the topological structure of binary tree storage vascular tree.

(2.2.1) be communicated with chain code and carry out the tracking of cardiovascular skeleton by eight:

Step1. search for the upper left corner of starting point from edge graph, from top to bottom, direction search from left to right.Get dir=7, and preserve a new profile origin coordinates (r, c).(wherein dir is the scanning direction variable, arrives the moving direction of working as fore boundary point along previous boundary point in the record previous step, and following table has been listed the corresponding different moving direction of different dir values).

dir	0	1	2	3	4	5	6	7
									X is to variation	1	1	0	-1	-1	-1	0	1
Y is to variation	0	1	1	1	0	-1	-1	-1

Step2. by 3 * 3 neighborhoods of counterclockwise searching for current pixel, its initial direction of search is set as follows:

If dir is an odd number, get (dir+7) mod 8;

If dir is an even number, get (dir+6) mod 8.

The pixel that first that searches in 3 * 3 neighborhoods is identical with current pixel value just is new boundary point A _n, more new variables dir is new direction value simultaneously.After current marginal point being carried out chain code coding, the pixel of current marginal point is changed to background pixel value, make and detect the edge next time and can not search again.

There has not been marginal point to exist if Step3. detect 3 * 3 neighborhoods of a marginal point, then this some distal point that is profile.Be stored in subsequently in the corresponding profile chain code.

After extracting vascular skeleton, measure the cardiovascular diameter according to the following step: use the disk template along skeleton direction finding diameter, extract the skeleton point place of diameter at needs, the diameter of disk template increases in suitable scope, when as long as a points of tangency occurring, this moment, the diameter of disk template was this blood vessel diameter when disk and marginal point.Utilize the seriality of vessel diameter change,, increase within the specific limits or dwindle blood vessel diameter, can reduce the hunting zone, improve search efficiency along the direction of skeleton by the skeleton chain code.

(2.2.2) x-ray imaging figure medium vessels topological structure is extracted

Represent the topological structure of blood vessel with tree, at first only need manually to specify the position of the root node of blood vessel, the tree that obtains blood vessel by binary tree search is represented then.In this process two problems can appear: first, because angiogram obtains by X ray transmission human body, may on projection, form intersection by the different blood vessel in the projection three-dimensional space, must identify intersecting blood vessels wherein, otherwise the result that binary tree search obtains not is correct blood vessel primitive projection; The second, because the decay that influence of the various factors in the actual x-ray imaging process and tissue produce X ray can be introduced a lot of noises in radiography figure, thereby may cause the appearance of pseudo-blood vessel, this also must be eliminated.

(2.2.2.1) identification of intersecting blood vessels

Ideally intersecting blood vessels if the blood vessel hop count that point connects on the blood vessel axis is 4, just judges that blood vessel occurs intersecting in this point in a bit.If only this situation occurs, also can remove intersection at an easy rate, promptly earlier vessel segment is disconnected, more corresponding vessel segment is coupled together.Yet under the practical situation, in the blood vessel axis extraction image of angiogram, the cross section of blood vessel may be a bit of blood vessel rather than a point, so just can not judge with top method simply.

Extract the image and can see from the blood vessel axis of angiogram, when intersecting blood vessels occurring, circulation can appear in binary tree search.The circulation sketch map appearred in binary tree search when Figure 17 was intersecting blood vessels.

Utilize these characteristics to judge, circulation time occurs in the binary tree search process, then explanation has blood vessel that intersection has taken place.

After in finding search procedure, circulation occurring, need to find the vessel segment that intersection takes place below:

The ring that the circulation vessel segment is formed takes place in record, and according to the characteristics of binary tree search, the blood vessel (i.e. the cross section of two vessel segments) that intersects should or be connected on the ring on ring.Blood vessel on the search ring or be connected blood vessel structure on the ring (promptly be connected on the vessel segment on the ring but not on ring) one by one, when they meet the following conditions, with it as intersecting the candidate of blood vessel:

Number(children)≥2，

Promptly the sub-blood vessel number when it just might be the blood vessel primitive projection that intersects more than or equal to 2, and cost function is set

f＝a ₁·Length(seg)+a ₂·min{fabs(angle(seg1，seg3)-π)+fabs(angle(seg2，seg4)-π)，

fabs(angle(seg1，seg4)-π)+fabs(angle(seg2，seg3)-π)}

A wherein ₁, a ₂Be weight coefficient, choose according to the normalization coefficient of blood vessel primitive length and vessel branch angle.Seg is current candidate vessels structure, and seg1 and seg2 be for being connected the blood vessel structure of seg on the ring, the blood vessel structure that seg3 and seg4 are connected for the seg another side, as shown in figure 18.

The candidate vessels structure of choosing cost function minimum wherein is as intersecting the projection of blood vessel primitive.Carry out flatness according to following function then and judge,

g ₁＝fabs(angle(seg1，seg3)-π)+fabs(angle(seg2，seg4)-π)

g ₂＝fabs(angle(seg1，seg4)-π)+fabs(angle(seg2，seg3)-π)

If g ₁＜g ₂, then seg1 and seg3 belong to same blood vessel, and seg2 and seg4 belong to same blood vessel; If g ₁〉=g ₂, then seg1 and seg4 belong to same blood vessel, and seg2 and seg3 belong to same blood vessel.

Two sections blood vessels that will belong to same blood vessel couple together and obtain a complete blood vessel, thereby remove the intersecting blood vessels that projection forms.

(2.2.2.2) identification of pseudo-blood vessel and removal

The pseudo-blood vessel discussed among the present invention is primarily aimed at and occurs in radiography figure, because the own organizational structure of human body (as rib), contrast agent distribution is inhomogeneous or inject the conduit of contrast agent process and reasons such as noise, and some similar blood vessels in image, have been stayed but not the figure of blood vessel.In the cutting procedure these figures are thought by mistake blood vessel segmentation comes out, formed the pseudo-blood vessel of mentioning in these chapters and sections.

Because the model among the present invention is just expressed the basic structure of blood vessel, there is not the yet unlikely individual difference of considering, therefore can't differentiate pseudo-blood vessel with model.Differentiate pseudo-blood vessel so utilize the topological structure of blood vessel to concern: generally speaking, topological structure was inequality when pseudo-blood vessel appeared in the radiography view of two angles; According to perspective geometry knowledge, projection can't change the topological structure (being meant qualitative but not quantitative annexation here) of blood vessel in addition.Based on above-mentioned condition, by relatively whether existing the different blood vessel structure of topological structure to judge whether this blood vessel structure is pseudo-blood vessel among two width of cloth radiography figure, as shown in figure 19.

(2.2.2.3) performing step

Consider top problem, the present invention has designed following steps:

1) begins search according to the method for binary tree search from the vascular radicle node, and the blood vessel pixel of searching for is made marks, when searching the blood vessel pixel that labelling crosses, then stop to search for forward, so just can avoid taking place endless loop;

2) search for whole tree, whether detection circulation occurs, as if circulation occurring, then finds out intersection blood vessel structure and correction according to the method for mentioning among the 3.1.1;

3) repeating step 2, till circulation does not occur.

(2.3) angiogram labelling

Blood vessel by among coupling angiogram and the vascular tree model projection figure carries out labelling to the blood vessel in the angiogram.The vascular tree model projection figure here is the mode of blood vessel three-dimensional model being copied the angioradiographic system projection, adopt to form parameter in the corresponding angiogram and carry out that projection obtains, therefore can judge whether they mate by the Topology Similarity between two sections blood vessels relatively.

Because individual difference, blood vessel quantity in angiogram and the vascular tree model projection and vascularity are also inequality.Because the vascular tree model of setting up is only to have got the significant blood vessel of arteria coronaria systematic comparison, can represent most of people's general character, suppose reasonably that therefore the blood vessel in the vascular tree model all has corresponding blood vessel in radiography figure.

At first need the selected yardstick that is used for instructing the vascular tree model of radiography figure medium vessels labelling,, select for use secondary vascular tree model that it is carried out labelling here owing to cut apart the blood vessel that obtains and few among selected patient's radiography figure.Blood vessel structure by in coupling angiogram and the vascular tree model carries out labelling to the blood vessel in the angiogram.

The blood vessel markers step:

(I) from different directions the vascular tree model is carried out the projection that projection obtains different visual angles by the x-ray imaging mode.

(II) ask the blood vessel axis of each projection and angiogram to extract the similarity of image T, find the projection P that wherein similarity is the highest.

Similarity is asked for by following steps:

A, the blood vessel topological structure that the blood vessel axis of vascular tree model projection figure and angiogram is extracted among the image T are kept at respectively among binary tree T1 and the T2;

B, from T2, extract the subtree T3 with T1 topological structure coupling.T1 and T3 are mated, calculate following cost function:

f＝a ₁·f ₁(angle)+a ₂·f ₂(curve)+a ₃·f ₃(length)

F wherein ₁(angle) represent the similarity of blood vessel topological structure between the two, be expressed as the function of angle angle between the vessel branch;

f ₂(curve) represent the similarity of blood vessel structure curve shape between the two, be expressed as the function of curvature of curve curve;

f ₃(length) similarity of two sections length of vessel of expression is expressed as the function of length of vessel length;

a ₁, a ₂, a ₃Be the weight coefficient of these three kinds of similaritys, choose according to the normalization coefficient of angle, blood vessel curvature and the length of vessel of vessel branch.The blood vessel of getting cost function minimum wherein is to being the coupling blood vessel.

C, get wherein that the subtree of cost function minimum is the matched children of T1, with the inverse of this minimum cost function as the similarity between model projection figure and the angiogram.

(III) obtain extracting the highest matched children M of projection P in T of image T similarity by the above-mentioned step of similarity of asking with the blood vessel axis of angiogram.Name with the projection of P medium vessels primitive is carried out labelling to the blood vessel primitive projection among the subtree M.Figure 20 is the PRELIMINARY RESULTS of carrying out labelling, (a) G1 among the figure, G2 ..., G7 respectively with (b) figure in H1, H2 ... the H7 labelling is identical.

(IV), be respectively labelling is carried out in the blood vessel primitive projection of corresponding radiography figure with the vascular tree model in these two moment as two radiography figure that choose not during the synchronization in cardiac cycle.

(2.4) blood vessel coupling

The blood vessel matching process comprises following two steps:

1) finds a pair of blood vessel primitive projection that same tag is arranged blood vessel primitive projection in two width of cloth radiography figure as coupling;

2) contrast two width of cloth radiography figure find a pair of blood vessel primitive projection of the identical blood vessel primitive projection of topological structure as coupling according to being limited in of topological structure in the blood vessel primitive projection that does not have labelling among two width of cloth radiography figure;

For the blood vessel structure that does not also have labelling among two width of cloth radiography figure, can look for the blood vessel of coupling right by the following method:

A, establish this two width of cloth radiography figure medium vessels topological structure and exist respectively among left_T and the right_T.

B, find and do not have the blood vessel of labelling UnName1 among the left_T, whether the search of same position place also has unlabelled blood vessel structure UnName2 in right_T, and may there be corresponding blood vessel structure in explanation if just have, and does not judge then that UnName1 is pseudo-blood vessel;

C, if UnName2 exists, whether one by one need judge between UnName1 and the UnName2 correspondence.Consider that in some cases may there be a more than unmarked blood vessel structure in the same position place among the radiography figure, be divided into two kinds of situations again and discuss:, then directly UnName1 and UnName2 are ordered with same tag if corresponding one by one; If not corresponding one by one, as shown in Figure 21, then need to mate between any two, the blood vessel that finds the cost function minimum is to being the coupling blood vessel, and makes same tag.

D, continuation step a, b, c is up to the whole labellings of the blood vessel in the vascular tree or till being judged as pseudo-blood vessel.

E, Figure 22 are to the last matching result of two width of cloth radiography figure, the vessel segment J1 among radiography Figure 22 (a) wherein, J2 ..., J13 respectively with 22 (b) in vessel segment K1, K2 ..., K13 mates corresponding.

Experimentize by the radiography figure to a lot of patients, three grades of blood vessels can correct labeling with interior trunk, and the following little blood vessel labelling accuracy rate of level Four is not very high, and the details of this and picture quality and model all has relation.Along with further replenishing of model detail, and the further optimization of model, the accuracy rate of little blood vessel labelling will improve gradually.

(5) three-dimensional reconstruction: further seek the correspondence between each pixel in the projection of coupling blood vessel primitive, blood vessel is carried out three-dimensional reconstruction according to reconstruction formula.

(5.1) before rebuilding, have any to consider earlier: when carrying out blood vessel 3 D reconstructing, two width of cloth radiography figure that choose must corresponding synchronization.And the image that we obtain generally is a radiography sequence image that obtains by single armed radiography system, all corresponding different time of the contrastographic picture in single armed radiography system.In order to satisfy two width of cloth radiography figure choose in unified requirement constantly, researcheres mostly adopt to be chosen ED angiogram and does three-dimensional reconstruction, but can slattern the radiography picture in other moment like this, may lose a lot of information.And do not consider like this effect of radiography figure yet, such as blood vessel block, noise etc., be difficult to guarantee that the angiogram of selected these different directions constantly is suitable for doing three-dimensional reconstruction.If can make other angiogram constantly all available, will increase the available information amount greatly, improve the precision of rebuilding.

As two width of cloth radiography figure that choose during not at synchronization, because their not corresponding immobilized hearts can not be used the general three dimension reconstruction formula, otherwise have big error.Therefore need elder generation that the blood vessel among this two width of cloth figure is unified as far as possible to the same moment, promptly carry out heart motion compensated:

Suppose to have chosen now two width of cloth radiography figure I ₁And I ₂, correspond respectively to the t in the cardiac cycle ₁, t ₂Constantly, I ₂Projected angle be θ.If I ₂In the blood vessel pixel be P, t ₁The time be engraved in the θ direction radiography figure be I ₃, P is at I ₃In corresponding point be P ', the process that is transformed to P ' by P is exactly motion compensation.

If heart model is at t ₁, t ₂The time radiography figure that is engraved on the angle θ be respectively J ₁, J ₂Point P is at J ₁In corresponding point be Q, the some P ' at J ₂In corresponding point be Q ', O is an initial point.Then can be by following formula compute motion compensated:

$\stackrel{\&RightArrow;}{{\mathrm{OP}}^{\&prime;}}-\stackrel{\&RightArrow;}{\mathrm{OP}}=\mathrm{scale}\&CenterDot;(\stackrel{\&RightArrow;}{{\mathrm{OQ}}^{\&prime;}}-\stackrel{\&RightArrow;}{\mathrm{OQ}}),$

Wherein scale is the normalization coefficient when radiography figure medium vessels is normalized to the vascular tree model.

With I ₂After motion compensation, obtain the new position I of blood vessel ₃Utilize I ₁And I ₃Carry out blood vessel 3 D reconstructing, wherein I ₁And I ₂The blood vessel matching relationship to I ₁And I ₃Stand good.

(5.2) based on the three-dimensional reconstruction of corresponding point

This section is utilized the discrete point reconstruction of three-dimensional coronary artery 3 d-dem point on the single face coronary angiography image of two different angles, briefly, just according to the corresponding point on the projected image of two width of cloth different angles, calculates the three-dimensional coordinate of spatial point.

(5.2.1) coordinate transform on the coronary angiography system three dimensions

Accompanying drawing 27 is radiography system imaging sketch maps two different angles.Utilize the rudimentary knowledge of coordinate transform, among Figure 27 from coordinate system x ₁y ₁z ₁s ₁To x ₂y ₂z ₂s ₂The geometric transformation campaign be three-dimensional rigid motion, can be with rotatablely moving and translational motion is described.And x in the motor process ₁y ₁s ₁Plane and s ₁z ₁The relative position of axle does not change, so only need make s by coordinate transform ₁z ₁Axle and s ₂z ₂Axle overlaps, and can make coordinate system x ₁y ₁z ₁s ₁And x2y ₂z ₂s ₂Overlap.

Describe for convenient, do following agreement:

(1) (x _1i, y _1i, z _1i) and (x _2i, y _2i, z _2i) an expression point P respectively _iAt three-dimensional system of coordinate x ₁y ₁z ₁s ₁And x ₂y ₂z ₂s ₂In coordinate;

(2) (u _1i, v _1i) and (u _2i, v _2i) an expression point P respectively _iIn projection plane coordinates is U ₁V ₁O ₁And U ₂V ₂O ₂In coordinate;

(3) x-ray source S ₁And S ₂And the distance between the center of rotation (isocenter) is respectively L ₁And L ₂

Be transformed into x by the XYZO coordinate system ₁y ₁z ₁s ₁Coordinate system makes ZO axle and s ₁z ₁Overlap, can be by twice rotation transformation, a translation transformation realization: (1) be around the OY axle α that turns clockwise ₁(2) be rotated counterclockwise β around the OX axle ₁(3) move to s from the O point ₁The point.

Promptly

[x ₁，y ₁，z ₁] ^T＝R _X(β ₁)·R _Y(-α ₁)·[X，Y，Z] ^T+T ₁ (1)

So

${[X,Y,Z]}^T=R_Y^{-1}(-{\mathrm{\&alpha;}}_1)\&CenterDot;R_X^{-1}\left({\mathrm{\&beta;}}_1\right)\&CenterDot;({[x_1,{\mathrm{<figure-callout\; id="1"\; label="y"\; filenames="A20081004785300444.png"\; state="\{\{state\}\}">y</figure-callout>}}_1,z_1]}^T-T_1)---\left(2\right)$

Can realize that the XYZO coordinate is tied to x by twice rotation transformation and a translation transformation equally ₂y ₂z ₂s ₂The conversion of coordinate system: (1) is around the OY axle α that turns clockwise ₂(2) around the OX axle β that turns clockwise ₂(3) move to s from the O point ₂The point.Promptly

[x ₂，y ₂，z ₂] ^T＝R _X(β ₂)·R _Y(-α ₂)·[X，Y，Z] ^T+T ₂ (3)

So:

${[X,Y,Z]}^T=R_Y^{-1}(-{\mathrm{\&alpha;}}_2)\&CenterDot;R_X^{-1}\left({\mathrm{\&beta;}}_2\right)\&CenterDot;({[x_2,y_2,z_2]}^T-T_2)---\left(4\right)$

Composite type (2) and (3):

[x ₂，y ₂，z ₂] ^T＝R([x ₁，y ₁，z ₁] ^T-t) (5)

Wherein $R=R_X\left({\mathrm{\&beta;}}_2\right)R_Y\left({\mathrm{\&alpha;}}_2\right)R_Y^{-1}\left({\mathrm{\&alpha;}}_1\right)R_X^{-1}(-{\mathrm{\&beta;}}_1)t=-R^{-1}{T}_2+{\mathrm{<figure-callout\; id="1"\; label="T"\; filenames="A20081004785300444.png"\; state="\{\{state\}\}">T</figure-callout>}}_1$

T ₁＝(0，0，L ₁) ^T T ₂＝(0，0，L ₂) ^T

3D method for reconstructing (5.2.2)

In Figure 27, false coordinate is x ₁y ₁z ₁s ₁The three-dimensional coordinate of middle coronary arterial tree skeleton point P is made as (x _1i, y _1i, z _1i), P is at image A coordinate system U ₁V ₁O ₁On subpoint p ₁Coordinate is (u _1i, v _1i); Coordinate system x ₂y ₂z ₂s ₂The three-dimensional coordinate of middle P is made as (x _2i, y _2i, z _2i), P is at image B coordinate system U ₂V ₂O ₂On subpoint p ₂Coordinate is (u _2i, v _2i).D ₁And D ₂Represent x-ray source S respectively ₁And S ₂Arrive the vertical dimension of projection plane separately.

If known (u _1i, v _1i) and (u _2i, v _2i), by following derivation, can be according to (R t) finds the solution the three-dimensional coordinate (x of P _1i, y _1i, z _1i):

According to the geometrical relationship of perspective projection, as can be known

$\left[\begin{array}{c}{x}_{1i}\\ y_{1i}\\ z_{1i}\end{array}\right]=z_{1i}\&CenterDot;\left[\begin{array}{c}{\mathrm{\&xi;}}_{1i}\\ {\mathrm{\&eta;}}_{1\mathrm{<figure-callout\; id="1"\; label="i"\; filenames="A20081004785300444.png"\; state="\{\{state\}\}">i</figure-callout>}}\\ 1\end{array}\right],$ $\left[\begin{array}{c}{x}_{2i}\\ y_{2i}\\ z_{2i}\end{array}\right]=z_{2i}\&CenterDot;\left[\begin{array}{c}{\mathrm{\&xi;}}_{2i}\\ {\mathrm{\&eta;}}_{2i}\\ 1\end{array}\right]---\left(6\right)$

Wherein

${\mathrm{\&xi;}}_{1i}=\frac{u_{1\mathrm{<figure-callout\; id="1"\; label="i\; D"\; filenames="A20081004785300444.png"\; state="\{\{state\}\}">i</figure-callout>}}}{D_{}},$ ${\mathrm{\&eta;}}_{1i}=\frac{v_{1i}}{D_1}---\left(7\right)$

${\mathrm{\&xi;}}_{2i}=\frac{u_{2i}}{D_2},$ ${\mathrm{\&eta;}}_{2i}=\frac{v_{2i}}{D_2}---\left(8\right)$

Further can obtain

$\frac{x_{1i}}{z_{1i}}=\frac{u_{1\mathrm{<figure-callout\; id="1"\; label="i\; D"\; filenames="A20081004785300444.png"\; state="\{\{state\}\}">i</figure-callout>}}}{D_{}}={\mathrm{\&xi;}}_{1i},$ $\frac{y_{1i}}{z_{1i}}=\frac{v_{1\mathrm{<figure-callout\; id="1"\; label="i\; D"\; filenames="A20081004785300444.png"\; state="\{\{state\}\}">i</figure-callout>}}}{D_{}}={\mathrm{\&eta;}}_{1i},$ $\frac{x_{2i}}{z_{2i}}=\frac{u_{2i}}{D_2}={\mathrm{\&xi;}}_{2i},$ $\frac{y_{2i}}{z_{2i}}=\frac{v_{2\mathrm{<figure-callout\; id="1"\; label="i\; D"\; filenames="A20081004785300444.png"\; state="\{\{state\}\}">i</figure-callout>}}}{D_{}}={\mathrm{\&eta;}}_{2i}---\left(9\right)$

Can obtain by formula (2)

$\left[\begin{array}{c}{x}_{2i}\\ y_{2i}\\ z_{2i}\end{array}\right]=R(\left[\begin{array}{c}{x}_{1i}\\ y_{1i}\\ z_{1i}\end{array}\right]-\stackrel{\&RightArrow;}{t})=\left[\begin{array}{c}({\mathrm{<figure-callout\; id="11"\; label="r"\; filenames="A20081004785300444.png"\; state="\{\{state\}\}">r</figure-callout>}}_{11}\&CenterDot;x_{1i}+r_{12}\&CenterDot;y_{1i}+r_{13}\&CenterDot;z_{1i})-({\mathrm{<figure-callout\; id="11"\; label="r"\; filenames="A20081004785300444.png"\; state="\{\{state\}\}">r</figure-callout>}}_{11}\&CenterDot;{\mathrm{<figure-callout\; id="1"\; label="t"\; filenames="A20081004785300444.png"\; state="\{\{state\}\}">t</figure-callout>}}_1+{\mathrm{<figure-callout\; id="12"\; label="r"\; filenames="A20081004785300521.png"\; state="\{\{state\}\}">r</figure-callout>}}_{12}\&CenterDot;t_2+{\mathrm{<figure-callout\; id="13"\; label="r"\; filenames="A20081004785300444.png,A20081004785300465.png"\; state="\{\{state\}\}">r</figure-callout>}}_{13}\&CenterDot;t_3)\\ ({\mathrm{<figure-callout\; id="21"\; label="r"\; filenames="A20081004785300444.png"\; state="\{\{state\}\}">r</figure-callout>}}_{21}\&CenterDot;x_{1i}+r_{22}\&CenterDot;y_{1i}+r_{23}\&CenterDot;z_{1i})-({\mathrm{<figure-callout\; id="21"\; label="r"\; filenames="A20081004785300444.png"\; state="\{\{state\}\}">r</figure-callout>}}_{21}\&CenterDot;{\mathrm{<figure-callout\; id="1"\; label="t"\; filenames="A20081004785300444.png"\; state="\{\{state\}\}">t</figure-callout>}}_1+r_{22}\&CenterDot;t_2+{\mathrm{<figure-callout\; id="23"\; label="r"\; filenames="A20081004785300444.png"\; state="\{\{state\}\}">r</figure-callout>}}_{23}\&CenterDot;t_3)\\ (r_{31}\&CenterDot;x_{1i}+r_{32}\&CenterDot;y_{1i}+r_{33}\&CenterDot;z_{1i})-(r_{31}\&CenterDot;{\mathrm{<figure-callout\; id="1"\; label="t"\; filenames="A20081004785300444.png"\; state="\{\{state\}\}">t</figure-callout>}}_1+r_{32}\&CenterDot;t_2+{\mathrm{<figure-callout\; id="33"\; label="r"\; filenames="A20081004785300394.png"\; state="\{\{state\}\}">r</figure-callout>}}_{33}\&CenterDot;t_3)\end{array}\right]---\left(10\right)$

$R=\left[\begin{array}{ccc}{r}_{11}& r_{12}& r_{13}\\ r_{21}& r_{22}& r_{23}\\ r_{31}& r_{32}& r_{33}\end{array}\right]$ $\stackrel{\&RightArrow;}{t}=\left[\begin{array}{c}{t}_1\\ t_2\\ {\mathrm{<figure-callout\; id="3"\; label="t"\; filenames="A20081004785300444.png"\; state="\{\{state\}\}">t</figure-callout>}}_3\end{array}\right]$

Simultaneous formula (9) (10) can solve:

$\left[\begin{array}{ccc}1& 0& -{\mathrm{\&xi;}}_{1\mathrm{<figure-callout\; id="0"\; label="i"\; filenames="A20081004785300394.png,A20081004785300401.png"\; state="\{\{state\}\}">i</figure-callout>}}\\ 0& 1& -{\mathrm{\&eta;}}_{1i}\\ r_{11}-r_{31}\&CenterDot;{\mathrm{\&xi;}}_{2i}& r_{12}-r_{32}\&CenterDot;{\mathrm{\&xi;}}_{2i}& r_{13}-{\mathrm{<figure-callout\; id="33"\; label="r"\; filenames="A20081004785300394.png"\; state="\{\{state\}\}">r</figure-callout>}}_{33}\&CenterDot;{\mathrm{\&xi;}}_{2i}\\ r_{21}-r_{31}\&CenterDot;{\mathrm{\&eta;}}_{2i}& r_{22}-r_{32}\&CenterDot;{\mathrm{\&eta;}}_{2i}& r_{23}-{\mathrm{<figure-callout\; id="33"\; label="r"\; filenames="A20081004785300394.png"\; state="\{\{state\}\}">r</figure-callout>}}_{33}\&CenterDot;{\mathrm{\&eta;}}_{2i}\end{array}\right]\&CenterDot;\left[\begin{array}{c}{x}_{1i}\\ y_{1i}\\ z_{1i}\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ \stackrel{\&RightArrow;}{a}\&CenterDot;\stackrel{\&RightArrow;}{t}\\ \stackrel{\&RightArrow;}{b}\&CenterDot;\stackrel{\&RightArrow;}{t}\end{array}\right]---\left(11\right)$

Wherein:

$\stackrel{\&RightArrow;}{a}=\left[\begin{array}{ccc}{r}_{11}-r_{31}\&CenterDot;{\mathrm{\&xi;}}_{2i}& r_{12}-r_{32}\&CenterDot;{\mathrm{\&xi;}}_{2i}& r_{13}-{\mathrm{<figure-callout\; id="33"\; label="r"\; filenames="A20081004785300394.png"\; state="\{\{state\}\}">r</figure-callout>}}_{33}\&CenterDot;{\mathrm{\&xi;}}_{2i}\end{array}\right]---\left(12\right)$

$\stackrel{\&RightArrow;}{b}=\left[\begin{array}{ccc}{r}_{21}-r_{31}\&CenterDot;{\mathrm{\&eta;}}_{2i}& r_{22}-r_{32}\&CenterDot;{\mathrm{\&eta;}}_{2i}& r_{23}-{\mathrm{<figure-callout\; id="33"\; label="r"\; filenames="A20081004785300394.png"\; state="\{\{state\}\}">r</figure-callout>}}_{33}\&CenterDot;{\mathrm{\&eta;}}_{2i}\end{array}\right]---\left(13\right)$

Formula (10) can be abbreviated as

A·C＝B

$A=\left[\begin{array}{ccc}1& 0& -{\mathrm{\&xi;}}_{1\mathrm{<figure-callout\; id="0"\; label="i"\; filenames="A20081004785300394.png,A20081004785300401.png"\; state="\{\{state\}\}">i</figure-callout>}}\\ 0& 1& -{\mathrm{\&eta;}}_{1i}\\ r_{11}-r_{31}\&CenterDot;{\mathrm{\&xi;}}_{2i}& r_{12}-r_{32}\&CenterDot;{\mathrm{\&xi;}}_{2i}& r_{13}-{\mathrm{<figure-callout\; id="33"\; label="r"\; filenames="A20081004785300394.png"\; state="\{\{state\}\}">r</figure-callout>}}_{33}\&CenterDot;{\mathrm{\&xi;}}_{2i}\\ r_{21}-r_{31}\&CenterDot;{\mathrm{\&eta;}}_{2i}& r_{22}-r_{32}\&CenterDot;{\mathrm{\&eta;}}_{2i}& r_{23}-{\mathrm{<figure-callout\; id="33"\; label="r"\; filenames="A20081004785300394.png"\; state="\{\{state\}\}">r</figure-callout>}}_{33}\&CenterDot;{\mathrm{\&eta;}}_{2i}\end{array}\right]---\left(14\right)$

$C=\left[\begin{array}{c}{x}_{1i}\\ y_{1i}\\ z_{1i}\end{array}\right]$ $B=\left[\begin{array}{c}0\\ 0\\ \stackrel{\&RightArrow;}{a}\&CenterDot;\stackrel{\&RightArrow;}{b}\\ \stackrel{\&RightArrow;}{b}\&CenterDot;\stackrel{\&RightArrow;}{t}\end{array}\right]$

By formula (14) as can be seen, this equation group is made up of 4 linear equation, finds the solution 3 unknown quantity x _1i, y _1i, z _1i, therefore being one transfinites and decides equation group, and its least square solution is:

C＝(A ^T·A) ^-1·A ^T·B (15)

Be coordinate system x ₁y ₁z ₁s ₁Three-dimensional coordinate (the x that middle P is ordered _1i, y _1i, z _1i).According to formula (15), just three-dimensional coordinate point can have been reconstructed.

Consider the defective of outer polar curve constraint, utilize the method for extract minutiae, the each point in a pair of blood vessel primitive projection of mating among two width of cloth radiography figure is mated, it is right to obtain each matched pixel point, and carries out three-dimensional reconstruction and obtain each space three-dimensional point.The space three-dimensional point that obtains with these reconstructions of B spline-fit obtains three-dimensional vascular tree.

In the accompanying drawing 23 (g) and (h) shown the angiogram three-dimensional reconstruction result that instructs algorithm for reconstructing and the reconstruction of outer polar curve matching algorithm to obtain with model respectively, be in two width of cloth radiography figure, manually to get the reconstructed results that corresponding point obtain (i).

In the accompanying drawing 24 (a) and (b) be that the blood vessel three dimensional structure that will rebuild with this paper method be Figure 23 (g) carries out back projection respectively by Figure 23 (a) and radiography angle (b) result, (c) and (d) be two axis images that original radiography figure extracts, respectively with Figure 23 (c) and (d) identical to being used for rebuilding.By comparison diagram (a) and (c), (b) and can see that (d) result of reconstruction still is reasonable.

(2.6) new model more: after reconstructing the blood vessel three dimensional structure, carry out perfect to model with the correct three-dimensional reconstruction result that obtains.Comprise employing all means available of following two aspects:

(2.6.1) three-dimensional reconstruction result and model are merged, make more generalization of model.Here merge and comprise following step:

1) finding the blood vessel primitive corresponding with the model medium vessels in the three-dimensional blood vessel structure that reconstruction obtains is the reference vessel primitive.Because corresponding blood vessel primitive is also so can determine on the three-dimensional blood vessel structure that reconstruction obtains for the corresponding relation of the projection that must be the blood vessel primitive when instructing the projection mark of radiography figure medium vessels primitive with model and the projection of model medium vessels, the blood vessel primitive on the model;

2) with the model be reference, the reference vessel primitive carried out corresponding blood vessel primitive mates in affine transformation and the model;

3) coordinate of the corresponding point on the corresponding blood vessel primitive in the reference vessel primitive after the affine transformation and the model is averaged obtains new vascular tree model.

(2.6.2) model detail is replenished.In model occurring during non-existent blood vessel structure, if it has all occurred in the result that the radiography figure to different patients rebuilds, then it is added in the model, and blood vessels all in the model is stored by progression, such as the one-level blood vessel, secondary blood vessel or the like.Set up the multiple dimensioned model of a vascular tree thus.According to different situations, select different yardstick models.

Accompanying drawing 25 is multiple dimensioned model sketch maps of our left coronary artery set up.Wherein: Figure 25 (a) is an one-level vascular tree model sketch map; Figure 25 (b) is a secondary vascular tree model sketch map; Figure 25 (c) is three grades of vascular tree model sketch maps;

Accompanying drawing 26 is we carry out subsidiary details to model with reconstructed results sketch maps.Wherein, be that we obtain left coronary artery blood vessel three dimensional structure by accompanying drawing 23 (a) two width of cloth radiography figure reconstruction (b) (a), different with Figure 23 (g) they are that Figure 26 (a) has used three grades of vascular tree models that blood vessel is instructed reconstruction.Figure 26 (b) is the result who Figure 26 (a) medium vessels axis is added radius, (c) be the vascular tree model of setting up according to visual human's slice of data, (d) be model in (c) to be carried out result after the subsidiary details with reconstructed results in (a), the blood vessel primitive of vessel segment that mark of circle wherein for replenishing.

Claims (5)

1, the angiography three-dimensional rebuilding method under a kind of dynamic model instructs, its step comprises:

(1.1) choose the X-ray coronary artery radiography figure of system of two width of cloth different angles, they are carried out pretreatment, at first adopt the enhanced method of frequency domain to strengthen blood vessel to each radiography figure image sequence, adopt morphology Bottom-Hat conversion further to emphasize blood vessel then, eliminate interference noise, keep and outstanding required image information;

(1.2) pretreated contrastographic picture is carried out the segmented extraction angiosomes, the refinement angiosomes obtains vascular skeleton, select for use eight connection chain codes to carry out the cardiovascular skeleton and follow the tracks of, extract vessel radius, adopt the topological structure of binary tree structure storage vascular tree according to following process;

(1.2.1) the blood vessel axis in angiogram extracts in the image, begin the blood vessel axis is carried out binary tree search from the vascular radicle node, and the blood vessel pixel of searching for made marks, show that when searching the blood vessel pixel that labelling crosses circulation appears in search, judging has blood vessel to intersect, and stop to search for forward, forward step (1.2.2) to; The appearance if search procedure does not circulate forwards step (1.2.4) to;

(1.2.2) circulation time occurs when binary tree search,, extract the blood vessel that identifies intersection in the image at the blood vessel axis of angiogram according to the topology information and the geological information of blood vessel;

(1.2.3) get back to step (1.2.1);

(1.2.4) search finishes;

(1.3) respectively the blood vessel axis extraction image of two width of cloth angiogram is carried out labelling by following process:

(1.3.1) from different directions the vascular tree model of choosing is carried out the projection that projection obtains different visual angles by the x-ray imaging mode;

(1.3.2) ask the blood vessel axis of each projection and angiogram to extract the similarity of image T, find the projection P that wherein similarity is the highest by following process;

(A1) the blood vessel topological structure among the blood vessel axis extraction image T of vascular tree model projection figure and angiogram is kept at respectively among binary tree T1 and the T2;

(A2) from binary tree T2, extract the subtree T3 that mates with the T1 topological structure; T1 and T3 are mated, utilize following formula to calculate the cost function of following T1 and T3:

f＝a ₁·f ₁(angle)+a ₂·f ₂(curve)+a ₃·f ₃(length)

F wherein ₁(angle) represent the similarity of blood vessel topological structure between the two, f ₂(curve) represent the similarity of blood vessel structure curve shape between the two, f ₃(length) similarity of two sections length of vessel of expression is expressed as the function of length of vessel length; a ₁, a ₂, a ₃Be the weight coefficient of these three kinds of similaritys, a ₁+ a ₂+ a ₃=1, the blood vessel of getting cost function minimum wherein is to being the coupling blood vessel;

(A3) get wherein that the subtree of cost function minimum is the matched children M of T1, with the inverse of this minimum cost function as the similarity between model projection figure and the angiogram;

(1.3.3) obtain extracting the highest matched children M of projection P in T of similarity of image T by the above-mentioned step of asking similarity with the blood vessel axis of angiogram; Name with the projection of P medium vessels primitive is carried out labelling to the blood vessel primitive projection among the subtree M, as two radiography figure that choose not during the synchronization in cardiac cycle, the vascular tree model with these two moment carries out labelling to the blood vessel primitive projection of corresponding radiography figure respectively;

(1.4) carry out the blood vessel coupling by following process:

(1.4.1) find the blood vessel primitive projection of same tag among two width of cloth radiography figure as the projection of coupling blood vessel primitive;

(1.4.2) contrast two width of cloth radiography figure find a pair of projection of the identical blood vessel primitive projection of topological structure as coupling blood vessel primitive according to being limited in of topological structure in the blood vessel primitive projection that does not have labelling among two width of cloth radiography figure; Its process is:

(B1) establishing this two width of cloth radiography figure medium vessels topological structure exists respectively among left_T and the right_T;

(B2) find and do not have the blood vessel of labelling UnName1 among the left_T, whether the search of same position place also has unlabelled blood vessel structure UnName2 in right_T, and may there be corresponding blood vessel structure in explanation if just have, and does not judge then that UnName1 is pseudo-blood vessel;

(B3), whether correspondingly one by one judge between UnName1 and the UnName2 if UnName2 exists; If corresponding one by one, then directly UnName1 and UnName2 are made same tag; If not corresponding one by one, then need to mate between any two, the blood vessel that finds the cost function minimum is to being the coupling blood vessel, and makes same tag;

(B4) repeating step (B1)-(B3) is up to the whole labellings of the blood vessel in the vascular tree or till being judged as pseudo-blood vessel;

(1.5) instruct the blood vessel axial point coordinate transform of will extract among two width of cloth radiography figure to synchronization by dynamic model; Seek the correspondence between the pixel in the projection of coupling blood vessel primitive, blood vessel is carried out three-dimensional reconstruction according to reconstruction formula.

2, angiography three-dimensional rebuilding method according to claim 1 is characterized in that: in completing steps (1.5), reconstruct the blood vessel three dimensional structure after, carry out perfect with rebuilding the correct blood vessel three dimensional structure obtain to model by following process:

(2.1) three-dimensional reconstruction result and model are merged, make the model generalization;

(2.2) model detail is replenished, comparison model and reconstructed results, when occurring in the model non-existent blood vessel structure in the reconstructed results, if its correctness of checking earlier correctly then it is added in the model, is improved the multiple dimensioned model of vascular tree.

3, angiography three-dimensional rebuilding method according to claim 2 is characterized in that: the fusion in the step (2.1) comprises following step:

(2.1.1) finding the blood vessel primitive corresponding with the model medium vessels in the three-dimensional blood vessel structure that reconstruction obtains is the reference vessel primitive;

(2.1.2) with the model be reference, the reference vessel primitive carried out corresponding blood vessel primitive mates in affine transformation and the model;

(2.1.3) coordinate of the corresponding point on the corresponding blood vessel primitive in the reference vessel primitive after the affine transformation and the model is averaged obtains new vascular tree model.

4, angiography three-dimensional rebuilding method according to claim 1 is characterized in that: the vascular tree model that step (1.3.1) is chosen can make up according to following process:

(c) set up the heart four-dimensional dynamic model according to following process:

(c1) from the heart sectioning image, be partitioned into the profile of pericardium and four chambers of heart;

(c2) set up the three-dimensional surface static models of pericardium and four chambers of heart respectively by following process:

(c21) pericardium or the heart chamber section profile that extract in the step (c1) are taken a sample, obtain sample point;

(c22) sample point to pericardium or heart chamber section profile carries out the B-spline surface match, sets up pericardium and heart chamber three-dimensional static model;

(c3) set up the motion model of all chambers of the heart chamber at first respectively, the motion of the motor control pericardium by all chambers of the heart chamber obtains the four-dimensional dynamic model of whole heart again;

(d) set up the four-dimensional dynamic model of coronary artery according to following process:

(d1) each ramose Vascular Slice profile of extraction coronary artery system from the heart sectioning image;

(d2) extract each ramose blood vessel axial point and radius of coronary artery system from the Vascular Slice profile;

(d3) each ramose axial point of coronary artery system of extracting is carried out the three-dimensional framework model that match obtains blood vessel with B-spline curves, and each blood vessel primitive on the model is named;

(d4) control motion coronarius by the periodic movement of heart, at first set up the motion model of coronary artery axis, on the motion model of axis, add radius information then and set up complete motion model coronarius.

5, angiography three-dimensional rebuilding method according to claim 4 is characterized in that: step (c3) is set up the heart four-dimensional dynamic model according to following process:

(c31) at first determine its residing position t0 in a cardiac cycle according to the form of the heart of the static models of setting up, then according to the kinetic characteristic of heart, the heart one-period is divided into 7 different periods, chooses the time point of the starting and ending in these seven periods and put t as 7 reference times _i(i=1 ..., 7), t wherein ₁Be the next time of t0,, make up the chamber dynamic model according to following step from t0:

(e1) make i=1;

(e2) at time t _iThe time, obtain with respect to time t according to curve V (t) _I-1Change in volume Δ V; The volume of the time dependent chamber of V (t) expression;

(e3) calculate the surface area S that each chamber surfaces does not comprise the chamber junction, described chamber junction is meant the position of atrioventricular valves (A V valves) and interventricular septum;

(e4) from time t _I-1To t _i, calculate each point because the displacement d=Δ V/S of the normal orientation that deformation produces;

(e5) to not comprising each sample point P of chamber junction on the curved surface, calculate surface normal direction at a P place

Order $\stackrel{\&RightArrow;}{{\mathrm{OP}}^{\&prime;}}=\stackrel{\&RightArrow;}{\mathrm{OP}}-d\&CenterDot;\stackrel{\&RightArrow;}{n},$ Wherein O is a zero, and some P ' is the reposition of the sample point after the deformation;

(e6) to new sample point reconstruct heart model, store new heart model and time corresponding t _i

(e7) move on to next time point t _I+1, repeating step (e2) is finished up to a cardiac cycle to (e6);

(e8) will be corresponding to putting t access time _iThe heart chamber model couple together the dynamic model that constitutes chamber, making the point on the heart chamber between adjacent time point is linear movement;

(c32) make up the pericardium dynamic model according to following step:

(f1) from time t ₀Beginning to each sample point P of pericardium, is calculated the normal in this point

Find normal

Cross point Q with chamber;

(f2) at time t=t _iThe time, obtain the displacement vector of a Q according to chamber dynamic modeling algorithm

(f3) according to the thickness d of the myocardium THICKNESS CALCULATION heart wall between pericardium and the chamber=| PQ| obtains new heart wall thickness then $d^{\&prime;}=d\&CenterDot;{\left(\frac{V_i}{V_{i+1}}\right)}^{1/6};$

(f4) order $\stackrel{\&RightArrow;}{{\mathrm{OP}}^{\&prime;}}=\stackrel{\&RightArrow;}{\mathrm{OP}}+(\stackrel{\&RightArrow;}{s\&CenterDot;}\stackrel{\&RightArrow;}{v}+d-d^{\&prime;})\&CenterDot;\stackrel{\&RightArrow;}{n}$ Be the sample point reposition vector after the deformation;

(f5) rebuild the pericardium model according to new sample point, preserve model and time corresponding t thereof _i

(f6) move to next time point t _I+1, the step (f1) that repeats modeling finishes up to a cardiac cycle to step (f5);

(f7) will be corresponding to putting t access time _iThe pericardium model couple together the dynamic model that constitutes pericardium, making the point on the pericardium between adjacent time point is linear movement;

(c33) dynamic model with pericardium and four heart chambers is combined into the heart dynamic model.

Images