FR2623642A1 - Method of segmenting vascular angiographic images by vectorial tracking of axes and servo-controlled detection of contours - Google Patents

Abstract

The method of segmenting angiographic images according to the invention consists, in a preliminary phase, in extracting from the image a skeleton of the network formed by a series of related points corresponding to the axes of the vessels, then, from this skeleton, carrying out a servo-controlled detection of the contours of the vessels by using the fact that a contour exists on each side of an axis and that these contours are also series of related points directly linked to the series of points constituting the axes: a dynamic programming is implemented in order to search, from matrices of gradients with dimensions limited to search areas close to the axis, the paths corresponding to the contours in the corresponding cost matrices resulting from the summing of the values of gradients. Application to the effective segmentation of angiographic images, especially images obtained by digital radiography, for constructing 2D databases which can be used for the reconstruction of 3D networks.

Description

Angiographic image segmentation method
vascular by axis vector tracking
and slave contour detection
The invention relates to the field of medical image processing, and in particular to digital angiographic images obtained by digital radiology systems.

 The segmentation of images from a radiology system, including digital, consists of extracting the image of useful information items. For angiographic images of vascular networks, this segmentation must allow to find in the image the arteries and veins of the network.

In this field, several methods have already been developed to perform the segmentation of angiographic images. A first method - described by NGUYEN TV -SKLANSKY J. "A fast skeieton-Finder for coronary
Arteries ", IEEE Comp, Society, No. 742, vol.1, pages 481-483,
October 1986 - uses, to look for characteristic pixels, two windows of limited size, one of sides parallel to the axes of the image, the other of sides oriented at 450 of those of the first. These two windows are moved over the entire image, and the maxima, that is to say the highest gray levels, are searched in these windows in four directions (two directions per window). If the number of maxima following a given direction is sufficient, these maxima are validated as part of the "crest" of a vessel.

The ridge is what can also be called the axis of the vessel in the image. The contours are obtained by applying masks, and the comparison between the contour elements and the axis elements makes it possible to eliminate a certain number of false detections.

Another method - described by HOFFMANN K. R et al. In an article entitled "Automated tracking of the vascular tree in DSA images using a sqiiare-box region of search, algorlthm", SPIE vol.626, Medicine XIV pages 326-33, 1986 consists in detecting the points of a vessel by scanning the perimeter of a square centered on an already known point of a vessel. For initialization, the first point is indicated by the user. The size of the window is variable according to the gray level of the current point. Bifurcations are searched on the perimeter of the square, outside the structure being studied. The detected vessels are marked by variably spaced points, and the width of the vessels is obtained by comparison of the actual profile of the vessel. image of the vessel with the profile of the projection of a cylindrical element modeling the vessel
In a general way, these methods lead to the extraction of characteristic points, considered as axis points or contour points, which do not necessarily have links with each other, except when we look for a sequence of maxima in a given direction in a predetermined window. The disadvantage of these systems is that they therefore produce discontinuous structures, and that it is then necessary to treat them, by means of complex calculations, to make them continuous. Moreover, the edge detecting elements are generally unreliable. and it is necessary, considering the number of points, to seek coherences to arrive at an approximately correct definition of the coordinates of the points of the vessels, to reconstitute the vascular network.

 The subject of the invention is a method for segmenting angiographic images, which does not have the aforementioned disadvantages and which, by design, necessarily leads to a continuous and closed structure. This structure has no holes, even at bifurcations.

According to the invention, the main features of the method of segmentation of vascular angiographic images are as follows
- Vessel axes are determined by vector tracking producing a sequence of related points
- Contour tracking is directly related to the axes determined in the previous phase, starting from the assumption that, on either side of the axis, there is an outline; contour tracking from this center line also leads to related points suites
- The close link between the axes and the contours, each contour being assigned to a perfectly identified axis, allows in a final phase the correction of the position of the axes.

 This method makes it possible to obtain a very good accuracy in the extracted data, and thus makes it possible to constitute a reliable database, in particular for the three-dimensional reconstruction of the vessels from a small number of projection images.

 In addition, the method according to the invention automatically and accurately extracts the vessels from a digital radiography, without pretreatment of the image.

Of course, the segmentation method according to the invention uses; for some of its phases; In particular, the method according to the invention also uses a window for detecting points of the axis, but only for a set of starting points. Moreover, an algorithm for finding a contour by dynamic programming has already been described - for example by POPE. D. L et al, in an article entitled "Left ventricular
Birder Recognition Using a Dynamic Search Algorithm, Radiology No. 155, pp. 513-518, 1985 -. The method according to the invention also uses a dynamic programming for contour search but which relies on the central lines.

The method according to the invention is optimized, in particular
- by refinement of the criterion of the choice of the direction of the exploration of a vessel, the calculation, and the tests
- by the integration of an automatic procedure for the detection of seed points on any structure related to that studied
- by the constraint of maintaining contour tracking with respect to the central line
- by the connectivity obtained both for the axes of the vessels and for the contours
- by assigning any contour to a perfectly identified branch
- And the possibility of rectification of the position of the axis that follows.

According to the invention, a method of segmentation of digitized angiographic images for extraction of the axes and contours of vascular networks, is characterized in that it comprises
a first phase during which a skeleton of the network is extracted, by a vector tracking method of the central lines of the vessels, or axes, searching, from isolated points recognized as belonging to the network by verification of an Ex criterion and each having a direction, the following one or more related points, in a limited number of predefined directions depending on the direction assigned to the current point, and which satisfy a predefined criterion Cd function of the gray levels of a set of predefined points,
and a phase of dependent detection of contours based on the skeleton of the grating obtained in the preceding phase to search, on either side of the axes, an outline also formed of sequences of related points, corresponding to important transition points grayscale search in predefined areas of the image.

 The invention will be better understood and other characteristics will become apparent with the aid of the description which follows with reference to the appended figures.

FIGS. 1 and 2 schematize the various steps of the image segmentation method according to the invention;
FIG. 3 is an explanatory diagram showing the directions of investigation
FIG. 4 is another explanatory diagram illustrating the data useful for monitoring the vessel
- Figure 5 illustrates the criterion of choice of direction
FIG. 6 illustrates the vectors selected for establishing a gradient matrix
FIGS. 7a and 7b illustrate the shapes of the contour search zones, around an axis segment, according to the inclination of the segment
FIG. 8 illustrates, in the simplest case, the chaining of the contour elements
FIG. 9 also illustrates the chaining of the contour elements in a case where the search zones are not joined.
FIG. 10 illustrates the treatment of the first segment of a branch.

 The following description of the angiographic image segmentation method according to the invention given by way of non-limiting example refers to FIGS. 1 and 2 respectively for the treatment steps relating to the extraction of the "backbone" of the vascular network, and for the processing steps relating to contour tracking.

 An image obtained by digital angiography, resulting from the scanning of an area analyzed by an X-ray beam, is a series of points of variable gray level, for example between 0 and 256 if the gray levels are given on 8 bits.

These rays are indeed absorbed differently depending on the tissues crossed. This image is organized in rows, of variable index i, and in columns, of variable index j. In general, the segmentation of vascular angiographic images by vector tracking of the axes according to the invention consists in extracting from this digital angiographic image a "skeleton" which, by hypothesis, is continuous, and which corresponds to the axes of the vessels of the vascular network. . Then, starting from the assumption that, on both sides of each axis. there is a contour of the corresponding vessel, we search the associated contours by a slave detection.

 For the creation of the skeleton, an initialization phase makes it possible to look for points in the image having a good probability of belonging to the skeleton. Because of the analysis method, the points belonging to the skeleton are points of "maximum" gray level, that is to say whose neighboring points, during a scan of a line or a column, have gray levels that are lower.

 In the following description, the gray level associated with a point is either the absolute gray level of the corresponding pixel, or the average of the gray levels of the pixels of a neighborhood, for example 3 × 3, centered on this pixel; this second measurement is preferable to overcome noise defects.

 For initialization, the first phase consists of looking for maxima on a line (or on a column of the image). A first scan of a line i (or of a column j) makes it possible to detect the points of maximum gray level, which therefore have a certain probability of belonging to a vessel.

 The next phase is then the possible validation of these maxima detected after a test sequence: from the average of the gray levels in a less vascularized area of the image, the "local level" of gray is calculated for the different maxima detected in the previous phase by difference with this average, and possibly taking into account the levels of the neighboring points. The averaged average MOY can be for example 118 for a gray scale with 256 levels.

The local level of the detected point can be defined as follows

 M (i, j + k) being the pixel located on the line i and on the column j + k and N (i, j + k) being its absolute level of gray; we thus take into account the neighboring points of the same line in the previous and next columns.

The maximum gray level points are retained only if
- NIV1 is increasing on a predetermined number of consecutive points, on a line,
- NIV1 becomes decreasing afterwards,
there is no previously detected maximum in a predetermined VOIS neighborhood of this point.

A maximum retained belongs to the skeleton of the network only if it satisfies then a criterion of belonging resulting from the study of a neighborhood of this maximum, or rather of the study of the neighborhoods of this maximum in the various possible directions from this point: a maximum is retained as belonging to the axis line of the network if the average of the gray levels of h points in an investigation direction
o is greater than a threshold, S; "h" characterizes the horizon
o of investigation, and S is a function of the mean, MOY, of a s-predetermined correction coefficient in the image, and of the standard deviation g of the gray level in a zone of little vascularity of the image. of a maximum retained in the previous phase, M (i, j) of level N (i, j), the 8 possible directions of investigation are represented in FIG. 3 and denoted Do to D7 with their directional coefficients in the direction i, m # (lines) and in the direction j, nx (columns).

The skeleton criterion, calculated in the direction D (x = 0 to 7) of the mx and nx coefficients
x is then

M (i, j) is a maximum belonging to the skeleton if at least one of the Ex, for x = 0 to 7 is greater than the threshold
x
S = MOY + s <r.

 From these validated pixels, the initialization phase continues for determining the orientation of the vessel to which this validated point belongs.

Let DxO (mO, nO) be the direction that maximizes E; he is
x logic to hypothesize that the vessel follows in the immediate vicinity of the pixel M (i, j) the orientation of the line dx having as one of its directions Dxo and as another direction D 8) The tracking of the vessel is done then in two steps, following the two opposite directions of the right dx
These two directions are thus kept in memory, with the maximum validated as belonging to the network.

 From at least one validated maximum, the next step of the initialization phase is the creation of a so-called "seed" point file. Indeed to obtain the continuous skeleton of the vascular network, it is necessary to have at least one point "germ" per branch. This file of points "sprouts" is obtained in the following way: starting from a validated maximum, and by a systematic exploration of the perimeter of a window centered on this maximum, the treatment consists in looking for all the maxima, as in the first step of the initialization phase; the maxima detected must then be validated, using the same test with respect to the threshold.

 Any structure close to the maximum of origin is thus marked by a point "germ". Each of these is then taken as the center of a new window, etc. . This method, using the arborescent nature of the network guarantees the marking of all the branches of the network by at least one "seed" point from a single maximum detected in the first step of the initialization phase, during the analysis. one line, or one column.

 This file of points "seeds" being created, the next step of the initialization phase possibly consists in a hierarchical ranking of these points "germs". The purpose of this classification is to avoid stopping the exploration of an important vessel, following the meeting of a secondary branch, while this secondary branch has already been detected. Moreover, it is interesting to allow the user to show the axes of the main vessels before the axes of the secondary vessels. The classification criterion used for this. Hierarchy of points "germs" is the gray level of a neighborhood centered on each of the points "germs". This neighborhood has a predetermined dimension, that is to say it has a number of points close to the validated maximum, and the points "seeds" are then classified in a file according to the gray level of this neighborhood. Thus, at the end of this operation the final file includes an ordered list of points classified according to the decreasing values of the gray levels of the neighborhoods.

These points will be the successive origins used during the tracking phase of the vessels, as will be explained below.

 The initialization phase being completed, the tracking phase of the vessels by their axes, intended to establish the skeletal file of the vascular axis can then begin.

For this phase, illustrated in FIG. 4, the following property is used: the image of a vessel element is approximately equivalent to the projection of a cylinder portion, and is therefore characterized by a peak line of gray between two flanks of substantially equal slopes. For vessel tracking, that is, to move from a current point to the next point, the analysis uses this property. At each of the current points, two decisions can be taken
- either the stop of the exploration if one arrives at the end of a vessel,
- or further tracking the vessel, which imposes a criterion for choosing the direction of travel to arrive at the next point.

Let M (i, 3) P 1 be the pixel retained at the previous iteration, and let M (1,3) P be the pixel retained at the iteration just completed: the direction (Dxo) P 1 made it possible to move from the point determined in step p1 at the point determined in step p.Exploration of the neighborhood for the search of the next point of the vessel is performed only in three directions ~ the previous direction Dxo and directions D and Dxd respectivelyà
xg xd respectlveme # t left and right with respect to the previous direction Dxo, and whose coefficients are perfectly determined knowing Dxo: Dxd (md ## dJ and D # g (mg, ng) This limitation in the search holds the fact that there can not be a sudden change of direction along a ship from one pixel to the next pixel.

In the step of tracking the vessel, the point M (i, j) can only be followed by one of the three neighboring points situated in these three directions, ie M (i + m0, j + n0), M (i + md, j + n #) and M (l + m, j + n), respectively labeled M1, M2 and
gg
M3 in Figure 4.

The following criteria are therefore used to determine vessel tracking
- criterion for stopping the exploration: for each of the three directions x defined previously, Dxo, Dxg and DXd, an average amplitude of the gray levels of a certain number of points situated in this direction, ie ho points (c ' that is to say up to the horizon ho), is calculated according to the same criterion Ex as previously, in the initialization step (see formula 1 above). The stop of the exploration is decided if the averages calculated in the three directions are all lower than the threshold S.The stop of the exploration is also decided if there exists a point already explored and different from the point ## MP 1 (i,) in a neighborhood of dimensions 3x3 centered on the point MP (i, j), that is to say if one encounters a bifurcation or a crossing with an axis already detected.

 choice of direction: the criterion chosen for the choice of direction, illustrated by FIG. 5, takes into account the property of the vessels to project along a peak line of the gray levels comprised between two flanks of substantially equal slopes; the decision criterion consists of a gray-level function component of the vessel axis and a function component of the difference between the gray levels on either side of the axis, on the edges of the vessels.

For this, the criterion is a function of the gray levels of the points situated in the direction of investigation and in two directions which are parallel to it at a predetermined distance from the axis.

FIG. 5 illustrates, for the first investigation direction Dx, the two orthogonal directions D (x + 2) mod 8 and D (x + 6) mod 8 and the parallels at a distance from the assumed axis.

où - h est, comme précédemment, l'horizon d'investigation
o dans la direction d'investigation, variable selon le niveau de gris du point courant ou d'un voisinage, c'est-à-dire fonction de l'épaisseur du vaisseau, et lo est un paramètre: 1 = i à ho;
o o
- m et n sont les coefficients directeurs de la direction pour laquelle le critère est calculé,
- i' = i+m'.np et j' = j+n'.np ; i" = i+m".np et jll = i+m".np, où m', n' et m", n" sont respectivement les coefficients directeurs selon les deux directions orthogonales à la direction testée,
- h1 est l'horizon fixe d'investigation sur les bords supposés du vaisseau de part et d'autre de l'axe à une distance np de l'axe, et 11 est un paramètre: li 11 à hi. The criterion calculated for the choice of direction is then

where - h is, as before, the horizon of investigation
o in the direction of investigation, variable according to the gray level of the current point or of a neighborhood, that is to say function of the thickness of the vessel, and lo is a parameter: 1 = i to ho;
oo
m and n are the directional coefficients of the direction for which the criterion is calculated,
i '= i + m'.np and j' = j + n'.np; i "= i + m" .np and jll = i + m ".np, where m ', n' and m", n "are respectively the directional coefficients according to the two directions orthogonal to the tested direction,
- h1 is the fixed horizon of investigation on the supposed edges of the vessel on either side of the axis at a distance np from the axis, and 11 is a parameter: li 11 to hi.

 k is a predetermined weighting coefficient.

This algorithm is called vector tracking algorithm because it takes into account pixels located on vectors in the same direction as the direction tested for the ship, along the assumed axis and along the edges. h
o the horizon of exploration in the direction of the sxe is a function of the gray level of the current point, and h1 represents the modulus of the vector on the edges, to verify the symmetry of the vessel with respect to the axis, is a fixed horizon comprising, for example four to six pixels.

 The chosen direction, among the 3 possible, is that which maximizes the criterion Cd.

 This algorithm is not very sensitive to noise and irregularities due for example to the poor distribution of the contrast medium in the vessels. It makes it possible to fill information gaps in the case of fine and low-contrast vessels, and it explores the structure on a horizon of several points, 4 to 10 points for example, before deciding on the displacement of a step on the picture ; the selected pixel does not necessarily have the largest gray level among those of the three possible pixels but, given the elaborstion of the criterion, it is the one most probably located on the axis of the vessel. The procedure is repeated for finding the next point on the vessel axis.

 When the stop of the exploration on a branch is decided, the next point "germ", in the ordered list of the points "seeds", is taken as base for the exploration of a new branch of the vascular network, if an axis already detected is not in its vicinity. The direction taken at the beginning of the tracking of a new branch from a seed point is one of those (opposite) stored in memory with the seed point, in the initialization phase.

 Thus, when the file of seed points is exhausted, all the pixels constituting the skeleton, that is to say the axis points of the vascular network in the image.

 This database being constituted, the next phase of the method consists in structuring this database so as to make the best use of the data resulting from the structure of the network, namely the connectivity and the weak variation of direction on each branch forming the network. An important property of the vector tracking of the axes is indeed, as indicated above, to produce a series of related points. The path of a branch requires only the exploration of a neighborhood restricted to a 3x3 window insofar as, by hypothesis, the axes are only one pixel thick and consist of the line joining a sequence of related pixels.

 The structuring of the database used for the search of contours is a decomposition into branches.

Each branch ends at the extremity of a vessel, that the current point has no successor, and begins at a junction, that is to say, where the current point has two successors. The decomposition is a little more complex than the decomposition into segments (a segment is here an element of axis between two bifurcations or crossing): the choice of the path to meet a bifurcation or a crossing is a function of a criterion taking into account the relative direction of the next segment and the gray level of the beginning thereof.

But since the branches are of longer length than the segments, contour tracking without breaking is easy. The decomposition process into branches can be detailed as follows
in an initialization phase, the origin of the common core is detected by searching in a predetermined area of the image of the highest gray level skeleton (the angiographies are taken at known angles, the tree on a view therefore always has a similar position);
- by traversing the branch including the common core, bifurcations are then detected: each is the beginning of a branch
- at the end of the course of the first branch the first fork is taken again for the course of the next branch and each new fork is memorized again
- by iterating this process, the whole tree is decomposed into branches from near to near.

 The next phase of the vascular angiographic image segmentation method according to the invention is then the controlled detection of contours, by dynamic programming, from the axes formed by the skeleton. Indeed, the axis of a vessel provides vital information about the contours: the axis is necessarily between two contours. The contour search area may therefore be restricted to the vicinity of the skeleton. The classical methods, ie the application of masks, of the first derived or second derivative type, are not sufficient to avoid interruptions in the contours. On the other hand, the use of dynamic programming guarantees the continuity of the contours, as the process used in the previous phase ensured the continuity of the axes.

 The dynamic programming, defined in general, is an algorithmic method allowing by successive approaches to find the optimal solution of a part of the global problem, the widening of the analyzed part leading to a new optimal path starting from the previous solution. For the controlled detection of contours, the dynamic programming algorithm is based on the search for an optimal path in a matrix of gradients whose dimension depends on the size of a segment, and which includes values representative of the variations of the levels. of gray in directions compatible with those of the axis of the analyzed branch. From this matrix of gradients, a cost matrix is created whose values result from the accumulation of information on paths defined in the gradient matrix, in moving from one line to the next by the route that maximizes the cumulative cost; then the optimal path in the cost matrix is detected, it corresponds to the contour.

The first step of this phase of contour detection is thus the creation of gradient matrices: MGRAD (x, y), from the stored information for a branch as defined above, a gradient matrix being defined for each segment of the branch
A branch is cut into elementary segments of length N5 pixels, NS being variable according to the curvature of the branch. The start and end coordinates of each segment are stored. The "direction" of each segment is defined by the Cartesian coordinates of its extremities. This division into segments makes it possible to obtain contour search zones in which the direction of the axis varies little and which have dimensions that are suitable for subsequent calculations, as explained below.

 Given the structure of the network, a branch has its origin in a bifurcation point detected on the central line of the network. Depending on the orientation of this vessel in three-dimensional space with respect to the image plane, the actual bifurcation may actually be behind, forward, or on the side of the actual vessel giving rise to the branch. Consequently, the detection of edges, or contours, at the birth of the branch is particularly delicate and will be performed as explained below in a particular manner.

 Regardless of this initialization problem, the gradient matrix has the following parameters: since the contour is always substantially parallel to the vessel's sue, the gradient is a function of the variation of the gray level on the one hand, and of the skeleton direction on the other hand.

To meet the first condition, the gradient calculation is performed on an elementary space having three side pixels.

The direction of this gradient can vary in significant proportions, particularly because of the presence of bifurcations or vessel crosses. The values of this parameter are then no longer representative of the variations in the gray levels induced by the studied branch but are in fact a combination of two vectors
a vector V1 linked to the branch under study
a vector V2 representative of the contour of the vessel outside this branch.

 The dimension of the domain used for the calculation of the gradient (3 pixels x 3 pixels) renders illusory the search for a great precision on the direction of the gradient obtained, and therefore on the exact contribution of the vector related to the branch studied. It is however interesting to reduce the influence of the vector V2 which tends to make the contour deviate with respect to the axis. To solve this problem, the values of the gradient matrix are the algebraic values of the projections of the gradients on the perpendicular to the segment representative of the axis. FIG. 6 shows the axis of the vessel and a segment AB on this axis, the gradient G at the point M (1, j), and the projection GP of this gradient on the perpendicular to the axis of the vessel, whose algebraic value is stored in the gradient matrix.

 This parameter being defined, another preliminary step consists in defining the search zone of the contour, starting from a given segment. This search area is different according to the inclination of the segment analyzed in the plane of the image. Indeed this zone is defined from the components of the segment and from the maximum width of a vessel of the vascular network being studied. Figures 7a and 7b show the two possible configurations of the search spaces of the contour around the analyzed elementary segment: in FIG. 7a, iB-3A is greater than iB-iA, whereas in FIG. 7b, iB # iA is greater than 3B-iA. If EC is the maximum deviation between the axis and the contour, that is the width of the search space on one side of the segment, and if LG is the length of the search space, the x and y coordinates of the gradient matrix MGRAD (x, y) are such that O <y <LG and O <x <EC.LG is equal to (3B-3A). possibly with a correction in the case of Figure 7a and LG is equal to (iB-iA) possibly with a correction in the case of Figure 7b. The correction will be defined below.

 The matrix of gradients MGRAD (x, y) contains as many values as there are points in each contour search zone, and the coordinates (i, j) of the center pixels of the elementary spaces (3x3) of calculation of the gradients in the image plane are stored.

 For each side of the segment, we have a matrix of gradients of LG lines and EC columns containing function values of the probability of presence of a contour at the point (x, y) corresponding to the coordinates (i, 3) of the image. The problem to solve is then to define a method to determine the most likely contour. Dynamic programming breaks down the search for the best path through the gradient matrix into a series of basic decision-making processes that take the step to be taken. to move from one row of the matrix to the next.

The first step of the algorithm is therefore to construct a cost matrix, MCOUT (x, y), from the gradient matrix MGRAD (x, y). This construction is carried out as follows
elementary paths are established between the first and the second line of the gradient matrix with a maximum allowed lateral deviation of the pixel
each point in the second row of the cost matrix takes the cumulative value corresponding to the path of the highest cost between the three possible paths
the process is repeated from line 2 to line 3 and so on until the last line of the gradient matrix of this search space.

The next step is then to construct the optimal path of a contour element. This construction is in the opposite direction of the accumulation of values in the cost matrix
- from the last line of the cost matrix we look for the maximum corresponding to the optimal path
- The search for the path to the previous line is done by selecting from one to the next the maximum among the three nearest points
- the procedure is repeated up to the first line of the cost matrix. It is stopped if the outer column of the cost matrix is reached
- the optimum end-of-path coordinates are stored.

The operations described above make it possible to search in a predefined search area for the optimal path of an outline element. There then remains a subsequent step of connecting together the contour elements thus detected. This operation is the chaining of the contour elements. Several cases are to be distinguished according to the rank of the segment in the branch. In order to simplify the description, the case of an intermediate segment of the branch is first described below.
The division of a branch into ordered segments was intended to reduce the dimensions of the search graph of the optimal path and to facilitate the maintenance of the constraint on the relative positions of the axes and contours. Indeed, if the segments are too long, that is to say, if the decomposition is insufficient, the assumption that a contour element exists on either side of the axis is no longer necessarily verified, especially when the direction of the vessel varies significantly; if the segment is too long the direction of the segment, only defined by its extremities, would be such that the two outlines would be on the same search area. But, because of this decomposition, it is necessary then to "chain" the contours without creating discontinuities at the connection points.

 Because of the method of detection of the axes constituting the skeleton of the vascular network, any segment has a successor except if one is at the end of a vessel. Therefore, the end coordinates of a segment correspond to the start coordinates of the next segment. Therefore, it will be the same for the contours, a contour element necessarily having an end coordinate identical to the start coordinate of the contour element corresponding to the next segment.

 The simplest case for the chaining of the segments and the corresponding contours is shown in FIG. 8: the contour search zones are joined because their widths are parallel (in FIG. 8 in the direction of the axis of the contours). i). In this case for the successive segments Si and the edge start of the contour associated with a search zone on one side of Si + 1 is the end of the contour associated with the previous search zone on the same side of Si: DSi + 1 = FSI and it is the same for the search areas on the other side of the segments.

Cumulative values in the cost matrix,
MCOUT (x, y) is normally in the opposite direction of the path of the branch, that is to say starting from the end of the branch. Thus, in FIG. 8 the sum of the values will be FS. to DS. (from the end to the beginning of the segment), the route of the
I i branch being performed from DS. at FS; it is the same for the
I i segment Si + 1 (cumulative FSi + 1 to DSi + i and goes from DSi + 1 to FS1 + 1). The last line of the matrix MGRADi + i linked to the segment Si + 1 corresponds to the vertical passing through DSi + 1.

Given the procedure used, the maximum on the last line of the MCOUT matrix. is not necessarily confused with the end of the contour associated with the previous segment, FC .. The maximum on the line ymax-l of MCOUTi + 1 can also be more than 1 point to the left or right of FC.

But the continuity of the contour is preferred: the path in the cost matrix characterizing the contour can therefore for two successive segments be slightly different from the optimal path, the connectivity constraint of the points of the path being maintained throughout the path.

On the other hand, when the direction of the segment varies with respect to that of the preceding segment Si so that we pass from #i> bj for Si, to ai <#j for Si + 1, as represented in FIG. or conversely from # i <# j to #i>#j), the two search zones passing through the ends of the segments are no longer contiguous parallel limits å j for Si and parallel to i for Si + 1. In that case
the end point of contour associated with Si, FCi is not on the first line of the gradient matrix MGRADi, which is the line passing through FSi,
on one side of the segment, the position of the contour point associated with Si, FCi, beginning of the contour associated with Si + 1, DCi + 1, leads to the reduction of the number of lines of MGRADi + i of L
LG = (jmax-jmin) + AL with bL <O.

On the other hand, on the other hand, the position of the end of contour point associated with Si, F'C1, requires an increase in the number of lines of MGRAD 'of a value fl L', so that
F'c. is included in the search box associated with the following segment Si + 1:
LG = (jmax-jmln) + ss L with tL> O.

To solve this problem of continuity of contours, the following improvements are made
- During the accumulation phase of the values in MCOUT if, on the FS line. there is no sufficiently consistent point with the presence of an outline, this line is rejected. The same is true for the following lines until a projected gradient value exceeds a threshold: all the following lines are then taken into consideration.

 - When traversing the path, if it approaches the end of a line (x = xmax - 1) in MCOUT (x, y), the contour is considered as leaving the field because of a curvature. The last validated point is taken as the end point of the contour Ci, FCi and as the initial point of the contour Ci + 1, DCi + 1, and this point also defines the last line of the matrix MGRAD. (as shown in Figure 8).

 The above description explains contour tracking for any segments of the branches except the first ones. Indeed their origins are bifurcations whose positions are defined imprecisely by the tracking of the axes when obtaining the skeleton; the actual start of the contours can be either at the edge of the vessel giving rise to a branch, or superimposed on a part of this vessel, this according to the relative directions of the vessel in the 3D space and in the projection plane. As a result, contour detection is more reliable at a distance from bifurcations.

 For this reason, and only for the first segment of each branch, the above method is retained, except for the direction of cumulation of the values and direction of the path in the cost matrix that are inverted with respect to those described above, as shown in FIG. 10, which represents a branch Bn with its axis and contours and a branch Bm originating at the point of bifurcation Bif on the axis of Bn: for S1 the accumulation of values is done at from the end of the segment FS i and this operation starts only from the first line where MGRAD contains at least one projected gradient value greater than a threshold and therefore possibly characterizing a contour point. path in MCOUT is carried out from the beginning of contour associated with S 2, DC2 even if this one does not coincide exactly with FC1 (cumulation maximum), for the reasons explained above.

 Thus, due to the method of the continuous contours directly related to the axes defined in the prior phase.

 A possible final step may then be to correct the position of the axes according to the contours thus detected, especially if in a subsequent image processing the data relating to the axes are used, for example for a three-dimensional reconstruction of the vascular network from projection images, the useful data of these projection images being limited to the axes. For this final step, it is sufficient to shift the axis of a number of pixels such that it is in the middle position with respect to the two previously defined outlines.

 Such a method achieves a reliable and highly efficient segmentation of X-ray images for all vascular network analysis applications, which is assumed to be continuous.

Claims (12)

 1. A method of segmentation of digitized angiographic images for extraction of the axes and contours of vascular networks, characterized in that it comprises  a first phase during which a skeleton of the network is extracted by a vector tracking method of the central lines of the vessels, or axes, searching, from isolated points recognized as belonging to the network by verification of an Ex criterion and affected each of a direction, the following one or more related points, in a limited number of predefined directions depending on the direction assigned to the current point, and which satisfy a predefined criterion Cd function of the gray levels of a set of predefined points,  and a phase of controlled detection of contours based on the skeleton of the grating obtained in the preceding phase to search, on either side of the axes, an outline also formed of sequences of related points corresponding to transition points gray level, searched in predefined areas of the image.  2. Method according to claim 1, characterized in that the phase of extraction of the skeleton comprises an initialization phase during which is formed a set of germ points clearly marking all the branches of the network, these points being recognized as belonging to the backbone of the network if, in at least one direction, defined by the tested point and one of its neighbors, the criterion Ey is

 where h, horizon of investigation parameterized by 1, is equal to a predefined number of pixels in the direction of investigation D  x of directional coefficients mx and nx, i and j are the coordinates of the current point M (i, j) in the image, and N (i, j) its gray level; the direction which maximizes Ex being assigned to the connected point of the current point in that direction recognized as belonging to the skeleton of the network.  3. Method according to claim 2, characterized in that the test points for the constitution of the set of seed points are, in an initial phase of the pixels corresponding to maxima of gray level along a line in the image, then the pixels corresponding to maxima of gray level ~ on the perimeter of windows successively centered on the germ points previously retained, all the branches of the network thus being marked at least by a seed point.  A method according to claim 3, characterized in that a current point is retained as corresponding to a gray level maximum along an image line or column portion for a maximum of the local level of gray, this local level being defined as the average of the absolute levels, of gray of three consecutive pixels on the line or the column, possibly diminished of the average level of a zone little vascularized of the image, this local level being increased on a predetermined number consecutive pixels, then becoming decreasing after the pixel retained.  5. Method according to one of claims 2 to 4, characterized in that the set of seed points is classified by decreasing values of neighboring gray levels centered on these seed points, to form an ordered file of seed points.  6. Method according to claim 1, characterized in that the tracking of the axes for obtaining the skeleton is achieved by selecting, from a previously selected point, only three possible search directions: the direction Do afffectée at the previous point and the two nearest directions on both sides of the direction Do, Dg and Dd, the following related point, if any, necessarily being in one of these three directions, the sequence of related points being interrupted by recognition of a

 is less than a threshold S for the three possible search directions, where h is the investigation horizon equal to  a predefined number of pixels in the direction of investigation D of direction coefficients mx and n I i and j the coordinates  x from the current point and N (l + mx1, i + nx1) the gray level of successive pixels of coordinates (i + m # l, j + n # 1).  7. Method according to claim 6, characterized in that the continuation of the follow-up, when the criterion E is greater  x at the threshold S for at least one of the three directions, leads to the selection of the next connected point among the three possible points, using a criterion of choice function on the one hand the persistence of a significant gray level characterizing the axis of a ship's crest in the direction of investigation and, on the other hand, of a sufficient symmetry on either side of the direction of investigation, on the edges of the vessels, parallel to the direction of investigation.  8. Method according to claim 7, characterized in that the next connected point is held in the direction which maximizes the criterion Cd such that

 h being the investigation horizon parameterized by 1 in the  o o direction of investigation of directional coefficients m and n, h being a fixed horizon parametrized by 11; (m ', n') and (m ", n") being the directional coefficients of the directions orthogonal to the direction of investigation, and i ', j', i ", j" being such that  i '= i + m'.np; i "= i + m" .np;  jt = j + n.np; j "= j + n" .np where np is the number of pixels between the axis line and the edge lines, where i and j are the coordinates of the starting point of the test.  9. Method according to claim 1, characterized in that for each of the points of the skeleton of the network, a database has its coordinates, the gray level calculated by the average of the levels in a neighborhood centered on this point, and the direction assigned to it, and in that this database is structured in branches defined between one end and one origin for the main branch and between one end and a branch point for the other branches, the criterion for determining the continuation of a branch at a junction or crossing point taking into account the directions assigned to the following points in relation to the preceding direction, and the gray levels of these following points.  10. Method according to claim 9, characterized in that the slave detection contours based on the skeleton of the network implements a dynamic programming for the search of these contours in areas of varying sizes located on either side axes of the network.  11. The method of claim 10, characterized in that, for the implementation of a dynamic programming, each branch is divided into segments of variable length according to the local curvature of the branch, and known start and end coordinates ,  in that the search zones on either side of the axis have the length of the segments as length and the width as a variable width according to the average gray level of the segment and greater than the maximum distance between the segment; axis and contour,  in that a matrix of gradients representative of the gray variations of the points, or their neighborhoods, relative to the direction of the segment is established for each search area,  in that a cost matrix resulting from the accumulation of gradient values, along the elementary paths of one line of the gradient matrix to the next, according to the path of the highest cost and with a maximum authorized lateral deviation of a point is established,  and in that the optimal path is followed from the maximum of the last row of the cost matrix to the first line, by selecting at each step the maximum of the three points closest to the preceding line, the procedure being stopped before the first line when the point is at the edge of the search area.  12. Method according to claim 11, characterized in that for any segment of a branch, with the exception of a first branch segment originating from a bifurcation point, the direction of accumulation in the gradient matrix for constructing the cost matrix starts from the end of the segment to go to its origin, the direction of construction of the optimal path characterizing the contour starting from the line of the matrix associated with the beginning of the segment to arrive at the line of the matrix associated with its end, these senses of accumulation and search of the optimal path being reversed for the first segments of the branches to take into account the uncertainties related to the bifurcations.