Tracking of Moving Coronary Artery Segment in Sequence of X-Ray Images

The variety of motion artefacts like situations when the patient is not able to hold his or her breath, natural organ movements and others cause great problems in current medical imaging techniques, such as image subtraction technique (eg. Digital Subtraction Angiography), registration technique, morphological and functional analysis, e.g. flow or perfusion estimation, and support of therapeutic intervention, eg. surgical tools and organ movement tracking. There are numerous reports in the literature mentioning possible types of motion artefacts. Many techniques are described to avoid and to correct such artefacts by means of digital image processing. The most important tasks utilized in all these techniques are the construction of a 2D geometrical transformation, which refers to original 3D deformation, finding correspondence between two images using only grey-level information, and efficient computer application of this correspondence [Meijering, 2000; Wang et al., 1999]. The procedure of finding the correspondence or estimation of the level of similarity between two images is the most important part of the above mentioned techniques, which has been applied either in static or dynamic approaches. This work presents an application of a correspondence finding technique to the densitometric analysis performed for coronary flow estimation on the image of a moving artery segment in X-ray image sequences of coronary arteries (Fig. 1).

techniques, such as image subtraction technique (eg.Digital Subtraction Angiography), registration technique, morphological and functional analysis, e.g.flow or perfusion estimation, and support of therapeutic intervention, eg.surgical tools and organ movement tracking.There are numerous reports in the literature mentioning possible types of motion artefacts.Many techniques are described to avoid and to correct such artefacts by means of digital image processing.The most important tasks utilized in all these techniques are the construction of a 2D geometrical transformation, which refers to original 3D deformation, finding correspondence between two images using only grey-level information, and efficient computer application of this correspondence [Meijering, 2000;Wang et al., 1999].The procedure of finding the correspondence or estimation of the level of similarity between two images is the most important part of the above mentioned techniques, which has been applied either in static or dynamic approaches.This work presents an application of a correspondence finding technique to the densitometric analysis performed for coronary flow estimation on the image of a moving artery segment in X-ray image sequences of coronary arteries (Fig. 1).Coronary image sequences are collected during coronarographic examination with X-ray contrast medium injection into coronary arteries for visualization purposes.Based on coronarographic sequences, both the morphometry of the arteries and the velocity of washing out of the indicator can be assessed [Nissen & Gurley, 1991;Goszczyńska & Rewicki, 2008].For densitometric analysis the measurements of the image gray levels of the pixels forming a measuring window (called the cross-sectional line) are performed for each image of a coronarographic sequence.The measurements are performed along the same cross-sectional line with its adjacent background.Fig. 2 shows the cross-sectional line positioned on the left coronary artery (LCA).where x is the coordinate within the cross-sectional line.Fig. 3b presents an example of the densitometric curve A(t) for a sequence of coronarographic images.Densitometric curve A(t) is a basis for coronary flow estimation in indicator dilution method (Appendix A).In coronary flow estimation method based on indicator dilution presented in [Goszczynska & Rewicki, 2008] the number of frames needed for the calculation of the densitometric curve A(t) depends on contrast medium injection parameters and the coronary flow rate, being usually larger than 100.Sections 1.1 covers problems concerning data collection from coronarographic image sequences.In Section 1.2 the techniques of object movement estimation in image sequences are briefly described.Section 2 presents a proposed method of automatic tracking of an artery segment and Section 3 brings the results obtained with this method.The error of the method is estimated in Section 4, and other solutions for the problem of artery movement tracing are mentioned in Section 5.

Issues related to automatic drawing of the cross-sectional line
The estimation of coronary blood flow, besides the difficulties concerning densitometric measurements, is affected by an additional problem of the cyclic movement of an artery segment or part of the myocardium.In the following sections we describe some issues related to the automatic and manual drawing of the artery cross-sectional line.In the manual case, the operator marks a cross-sectional line at the same point of the artery on each frame of the sequence.Unfortunately the phrase "the same point" is not precise due to at least two reasons: • the myocardium movement results in the changes of the location and shape of the particular part of the artery which makes ensuring that the cross-sectional line on the frame is in the same place as in the previous frame difficult, • the subjective nature of manual positioning of the same point, in the absence of some characteristic points, is very error-prone, even when the location and the shape of the object remain unchanged.The above issues suggest that the process of drawing the cross-sectional line should be automated.The manual solution, besides being inaccurate, is also very exhaustive and boring for the operator -the cross-sectional line has to drawn on at least hundred of consecutive frames for a single measurement.For a clinical study of the flow measurement method presented in [Goszczynska & Rewicki, 2008;Goszczynska, 2008] the analysis of more than ten thousand images was required, and that is why certain type of automatic aid for cross-sectional line positioning has been strongly suggested.The algorithms performing the task mentioned above are described as movement correction (the 3D transformation of each frame in sequence resulting in apparent lack of heart movement) or tracing the object (relocation of the ROI, so that the particular detail of the heart image -the artery fragment -is traced).The movement correction algorithms may be applied to local areas of the whole frame.For heart images processing movement correction for the whole frame is used to reduce the distortions made by patient movement.This makes it impossible to precisely determine a mask frame within the DSA technique [Meijering, 1999;Meijering, 2000].
The crucial difficulty preventing an effective incorporation of those algorithms to coronarographic frames analysis is a low densitometric contrast of the heart areas of interest.As an example an artery fragment image can be considered, across the whole sequence of the frames, also covering frames before contrast medium injection, and during contrast medium washing.In the instances when the area of interest covers only a part of the myocardium with tiny arteries, the frame analysis results in the analysis of the areas with nearly homogenous brightness, similar to the surrounding area.In such situation the movement detection can be performed indirectly, by tracing the movement of the nearest artery segment and applying a 2D interpolation.In the analysis of movement in frames containing the image of clearly distinguishable arteries another difficulty appears.Local examination of continuous movement for contour in the plane perpendicular to the image detection plane may be hardly possible.The local movement analysis of a continuous contour had first been described by Stumpf in 1911 [Hildreth, 1987].Obviously the observed movement of an artery using 2D systems lacks component of the real movement perpendicular to the image plane (the observed movement is barely visible).Computing any corrections of the event on the basis of only single plane images is a difficult task.

Movement estimation techniques
Artery segmentation and/or contour and central line detection methods are frequently implemented as a preliminary, but not mandatory step, in the movement detection algorithms described in literature [Coatrieux, 1995;Rus, 1995;Benayoun & Ayache, 1998;Konrad, 2000].These methods can be divided into two categories: • those limited to the area being a subject of segmentation only, tracing the artery assuming the continuity of the arteries tree structure, and • those analyzing the complete image and sharpening the structures investigated, so that the background can be removed.The first step of the analysis utilizing the methods of the first category is manual positioning of the point of the center line or on the contour of the element investigated.Then the operator proceeds in certain direction according to some criteria, and positions the next points of the line or contour.An example of such a criterion may be the maximum gradient of the grayscale.The methods of the second category are generally two types of filtration -linear and nonlinear.The linear filters utilized for contour detection are based on testing two-dimensional derivatives of the first degree (Sobel or Roberts operators) or the second degree (Marr-Hildreth and Laplace operators) for the brightness change function within the selected neighborhood of an indicated point or, in a more general way, on the operation analysis between the defined operator mask and the neighborhood of a certain point.The adaptive linear filters, meaning those in which a mask is adapted to the shape of the object tested, maximize the output value for the central line of the object.The morphological filters belonging to the nonlinear filter class are based on the comparison of the object tested with a structuring element.For segmentation purposes, opening operation i.e. erosion and dilation is usually used.Some other methods may be mentioned here like image decomposition by eigenvalues, a curvature analysis, vector detection of the central line and the method of geometric moments.The non-subtractive images, in opposite to subtractive ones are distinguished by high nonhomogeneity of the background, which causes the artery segmentation to be a very difficult task.Due to this difficulty the main criterion for estimation of the results provided by filters is their ability to equalize the background level.The artery segmentation phase, setting their contours and central lines increases the effectiveness of movement detection methods implemented in the solutions, which are more sensitive to the structure and less sensitive to the brightness.Unfortunately it also increases the time of the movement analysis process.The local movement detection methods or movement measurement methods in images may be generally divided into two categories: an optical flow procedure, based on the brightness gradient analysis, and a template matching procedure, based on image similarity measures.The first procedure, suggested by Horn and Schnuck [Horn, 1981], is based on the continuity principle, also known in physics, in the movement analysis of the frame sequence is described by the formula: where I is the brightness of the pixel tested, v x is the x component of the tested pixel velocity, and v y is the y component of the pixel tested velocity.The pixel velocity is determined according to the spatial and time resolution of the image tested sequence.By such processing we obtain the estimation of the apparent movement at every point of the frame.Application of differential techniques based on this method results in a relatively fast algorithm.There are two important conditions (considered disadvantages) that have to be met: the brightness of the moving object should be constant in time and the movement of the object between two consecutive frames should be relatively small.In the case of arteries, both of the mentioned conditions are not fulfilled for the entire frame sequence.Based on some presumptions, this technique may provide acceptable results for the estimation of the general myocardium shape change during its working cycle [Meunier et al, 1989].
In the next chapter we will describe a method of automatic movement trajectory estimation for an artery segment based on template matching.

Method
In the frame sequence studied for marked artery segment three stages can be indicated: • before contrast medium appearance -T a , • with contrast medium filling -T c and • the contrast medium wash-out period -T p .The movement analysis of a fragment of the artery under investigation was divided into three steps: • movement estimation for the marked fragment of the artery during contrast medium filling period T c , • correction of faulty detections, • movement estimation in stages T a and T p .

Estimation of movement of an artery segment in the indicator filling period T c
Below a method of automatic data collection based on a popular block-matching technique is described [Gonzalez & Wintz, 1987;Rong, 1989;Maragos & Pessoa, 2000;Buzug, 1998].This method was used in the similarity analysis of a few selected fragments of original frames sequence (without performing any pre-processing) with a chosen template, i.e. the fragment of an artery as it was shown on the frame where the arteries are clearly visible (which means the end of the contrast medium injection phase).The sum of difference squares C SSD was used (R, S denotes the compared fragments of the images in the area of m, n dimensions): In order to find the area in frame I(t+1) which is similar to the given area R in frame I(t), the area S is defined in frame I(t+1) with the same coordinates and size as the area R in frame I(t).The next step is the area enlargement by values k max and l max which are greater or equal to the maximum movement of the selected fragment between the two frames (Fig. 4).
The maximum value C max of this coefficient tells the location of the area most similar to the reference area.The indexes k and l corresponding to the maximum value of the coefficient determine the coordinates of the area investigated (Fig. 4).
Other measures of similarity have been also proposed: • variance normalized correlation: This measure refers to a nonlinear signal correlation.As the following equations is true |a-b|=a+b-2min(a,b), then under some conditions (i.e. when the correlation value is normalized by dividing it by the average brightness value for windows R and S), the minimum value of C SDA is equal to the maximum value of the morphological mutual correlation CM defined by equation ( 7): 11 00 m i n ( (, ) , (, ) ) • mutual information method , incorporating entropy E: 0 () l o g() where P(g i ) is the probability of occurrence for the i th grade on the grayscale in the range of 0-I.The mutual information coefficient denoted by C v is defined by equation ( 9): 00 (, ) (, ) l o g () () Pg g Cv P g g P g P g where P(g i ,g j ) is the probability of occurrence of the pixel with the i-level gray in one frame and the pixel with the same coordinates and j-level gray in the range of 0-J in the second frame.

Correction of faulty detections and movement estimation in periods T a and T p
In the measurements of coronary flow with automatic data collection method [Goszczyńska, 2004] manual correction was applied for faulty artery movement detections during the contrast period T c (medium presence) using time-related interpolation.To draw the artery movement trajectory in the period T a (before contrast medium injection) and during the period T p (contrast medium washing-out) a manual method of time extrapolation was applied.

Algorithm and technical solution
Fig. 5 describes the semiautomatic procedure of establishing a cross-sectional line for an artery in the N-frame sequence.It starts with the operator establishing frame n r where he positions window R containing the template and determines the size S of the area to be investigated in the process (see Fig. 4).Next the procedure is initiated to find the values for indexes k, l (for each frame in sequence), for which the similarity measure C of the window R and subarea of window S reaches the maximum value.After the operator positioned the cross-section line l with k pixels in the frame n r , the procedure determining the coordinates of cross-sectional lines on N frames of the sequence is initiated according to the index table computed in the previous stage.The next stage is the automatic brightness data collection along the preset lines (Fig. 3a -line D a (x)).
The MS Windows application YFlow (Fig. 6) implementing the automatic movement tracking of a chosen artery fragment using the method of template matching has been developed (the C++ language, authors: H. Goszczyńska and L. Kowalczyk) [Goszczynska et al., 2006].

Results
Fig. 7 shows an example of a proper matching performed with mutual correlation method for all frames in the interesting range of frames, that is from the moment of contrast appearance until it is washed-out from the chosen artery segment.

1-25 26-50
51-75 76-100 Fig. 7. Parts of 100 images from a negative coronarographic film with an automatically marked cross-sectional line There were seventy coronarographic sequences analyzed in [Goszczynska, 2004] (for 16 patients, 4-5 sequences per patient).A single case of faulty matching in the analyzed material was observed.In [Goszczynska, 2008], examples of mismatch of template R with window S in a coronarographic sequence (Fig. 8) and the comparison of the matching results calculated for the matching concordance measurements listed in Section 2.1 are presented (Fig. 9).

Discussion
The occurrence of faulty matching and computational complexity are the main problems of the method of automatic tracing of the artery segment movement.The most frequent reason for mismatching is too small size of searching area S. The proper p o s i t i o n i n g o f t h i s a r e a m u s t b e d o n e a c c o rding to the results of the analysis of the maximum movements of the R mask with a chosen artery segment.Such faults can easily be avoided during the analysis by choosing a bigger searching area (taking care not to include more objects similar to the one which is necessary for the detection).However this operation requires more time for the analysis, so for calculation of C max S could be kept small but displaced slightly for each consecutive frame, based on the previously estimated movement.This solution, unfortunately falls into the same trap as that described below.The changing template algorithm has also been tested.As a primary template a fragment of an artery from the frame at the end of injection was selected.In such case, it works as a template only for comparison with the previous and the next frame.For the remaining frames the template is a result of the previous searching.Such an algorithm, intuitively better, still has a trap: due to an accidental mismatching the template may change entirely.Various solutions of this problem are proposed in the literature, mainly based on the reliability of the match [Guest, 1999].
Besides the problem of finding the best template or templates when choosing an optimal window size (m, n) containing the template and an optimal size of the searching area k max and l max , the detection faults may also be reduced by adopting a stronger criterion of similarity.The tests we have conducted showed that the best correlation coefficient is a normalized one.A morphological mutual correlation (a sharper matching peak) works better than a linear mutual correlation in detection of the area R within area S.Moreover, the morphological correlation algorithm (the sum of minima) is faster than the linear correlation algorithm (sum of products).Yet it is not determined by our tests what are the advantages of one algorithm over the other in regards to the precision of the matching.
The results obtained with our method should be compared with those based on other criteria, when for example two templates are used or other similarity measures, i.e. those based on histogram analysis applied for the difference (subtraction) of two adjacent frames [Buzug et al., 1998].
Our algorithm inability to analyze fragments of the object with no characteristic morphological details such as bifurcations or curvatures is a known issue.The solution here would be the incorporation of 3D interpolation methods [Schrijver, 2002].
Regarding the problem of faulty matching correction a 2D interpolation for the template positioning based on the neighbor proper matching seems to be a better solution than the one suggested in the present study.Yet specifying the point for which the automatic template matching was performed, the process of replacing densitometric data was easier and led to a smaller error (this is due to the continuity of the process of changing the indicator concentration in the artery cross-section investigated).
Controlling of the faulty matching can be computer-assisted, e.g. by comparing the values of indexes (l n ,k n ) obtained with the range limit, and finding the maximum movement of the artery fragment during the heart cycle, or by the analysis of movement direction, i.e. the direction of the vector defined by indexes (l n-1 ,k n-1 ) and (l n ,k n ).
A proper matching for the T p period is essential in the case when a part of the myocardium perfusion is estimated, which movement, on the basis of the closest artery displacement, is calculated.The degree of perfusion can be estimated from the image brightness changes of the part of the myocardium.The analysis of an artery movement trajectory may also allow to estimate changes in duration of the heart cycle phases and may serve as a basis for modeling epicardial strains [Young & Hunter, 1992].In the next section described is the estimation of the error in the densitometric curve calculation caused by the use of the automatic procedure for data collection.

Method of error estimation
The presented here application of the correspondence finding technique was used in coronary flow estimation by computer analysis of angiographic image sequences.For the estimation of the error caused by the automatic artery tracking method one may define as the reference either the values of the coordinates obtained during manual cross-sectional line drawing or other values calculated during the estimation of the flow rate.The following options for the reference value were analyzed: • the coordinates of the P 1 and P 2 ends of the cross-sectional line (Fig. 11a) The following section discusses the analysis of the error caused by the automatic drawing of the cross-sectional line for the three above listed reference values.

The coordinates of cross-sectional line ends
The coordinates of the cross-sectional line end values were used for the estimation of the slope angle of the positioned cross-sectional line and for the construction of the movement trajectory.For image sequence of RCA (Fig. 10b), the changes in the slope angle of a manually positioned cross-sectional line with the axis of the artery segment were measured.It allowed to determine the error level of the automatic analysis caused by the fact that for the implemented algorithm the cross-sectional line has the same slope.For the example tested the difference between the two methods was not greater than 14° (Fig. 12a).
There was also the influence of those changes on the length of the cross-sectional line for one pixel (Fig. 12b) calculated (with a unit being the pixel size).As presented in Fig. 12c for the three positions of the line its length expressed in pixels is constant.This is due to the line generation algorithm in the discrete space (for the simplest position of the artery axis on the digitized plane).Thus in our case the fact that a constant slope of the cross-sectional line is assumed has no influence on the calculation results, that is on the range of the integration in equation ( 1).Despite the number of pixels in the segment remaining constant, the pixels that create the segment are not identical, which is also a source of error.
The only exception to the above results would be such a positioning of the cross-sectional line that the range of change of 14° covers the line drawn at 45° with the artery axis, which in practice means a wrong positioning of the cross-sectional line.

The value of field A
The changes in the background of the artery fragment investigated make the value of the A field sensitive to the slope angle of the cross-sectional line and the artery axis and to the interpolation method of the background lines [Goszczynska, 2004].Fig. 13 shows two densitometric curves for the same cross-section of a catheter filled with a contrast medium, computed with two background line interpolation algorithms.Although method I (Fig. 13a) gives larger errors compared with method II (Fig. 13b), the final result, that is for the field under the A(t) curve, is computed with significantly smaller error.The conclusion based on the above results was to use the calculation of the field under the A(t) curve as the reference values for the estimation of accuracy of automatic cross-sectional line positioning.

The area under the densitometric curve A(t)
To estimate the error resulting from automatic data collection, densitometric analyses of the same artery segment were performed for automatic and manual positioning of the crosssectional lines.The lines were positioned in three ways: automatically correctly (i.e.perpendicularly to the artery axis) on the pattern image, automatically purposely incorrectly (i.e.not perpendicularly to the artery axis) and manually correctly (Fig. 14).where max c , min c are the maximum and minimum values of the area for "correct" lines, min m is the minimum area value for "manual" lines.Error of type 1 was estimated as less than 11%.This type of error is associated with the imperfection of the matching algorithm, e.g. when rotation is disregarded in equation ( 4).Error of type 2 (also called the total error of the flow estimation method) calculated for all the measurements performed was defined as:

Conclusion
The algorithm described makes it possible to automatically collect the data used for the coronary blood flow measurement based on the analysis of a diluted indicator; though some occurrences of mismatch happen.The assessed error of the method did not exceed 11%.The algorithm of automatic data collection for the coronary flow estimation significantly accelerates the collection of the data without deteriorating the accuracy of the results obtained.Our method of automatic tracing of the measurement area not only reduces the manual effort of the operator involved in positioning the cross-sectional line in the frame sequence, but also provides a simple method of obtaining a densitometric curve with the time resolution equal to that of the timeline of the sequence investigated.Most of the other approaches to the measurements of coronary flow are based on the analysis of densitometric curves calculated from the frames referring to a particular phase of heart cycle and are later interpolated with a preset fitting curve.The presented algorithm makes it also possible to correct mismatching and incorporation of movement trajectories obtained to analyze the sequence range in which there is no contrast medium in the artery segment investigated (invisible in the frame).This permits perfusion estimation of the part of the myocardium which surrounds the artery segment.For the artery (LCA and RCA) analyzed segments not taking the rotations into consideration has no influence on the accuracy of the collected data due to the relatively small (<45°) rotation movement of the cross-section investigated during the heart cycle and the features of algorithms generating straight lines in the digitized space.It seems that this relatively simple algorithm can be accepted even in routine clinical tests due to a short time (few minutes) of the sequence analysis.Some other ideas related to the solutions of the problem of effective artery segment tracing can be mentioned here: • Inclusion of some searching strategy into the matching algorithm to reduce the number of positions of the template for comparison with a selected fragment [Orkisz, 1999].

Appendix A Coronary flow estimation using densitometric measurements
The calculation of the indicator concentration from a radiographic image is based on image brightness measurements.The differences in brightness in X-ray images are caused by spatially different attenuations of the X-ray beam by the tissue, as described by the Lambert-Beer law: where I o is the intensity of the monochromatic input beam, I 1 is the intensity of the attenuated beam at a distance l [cm], and s is the linear X-ray attenuation coefficient [cm -1 ].In angiographic images, morphometric and haemodynamic measurements are performed mainly on some parts of the image in the region of uptake of the contrast medium, i.e. in arteries or the myocardium filled with an X-ray indicator, whose attenuation coefficient s is significantly higher than that of the blood (the blood attenuation coefficient is neglected).The value of the distance l can be expressed as: For the region of interest (ROI) consisting of P ROI pixels, V ROIi , equation (A4) becomes a summation calculated for each pixel with coordinates m,n ⊂ ROI:

Fig. 1 .
Fig. 1.Sequence of coronarographic images (every 6 th image) from negative coronarographic film with visible contrast medium passing

Fig. 2 .Fig. 3 .
Fig. 2. Fragments of images of the left coronary artery (LCA) from a coronarographic sequence during injection, filling and washing out of a contrast medium

Fig. 4 .
Fig. 4. Determining the template and investigation windows The search for the area in frame I(t+1) which is most similar to the chosen area R in frame I(t) involves moving the reference area R on the magnified area S and computing the coefficient of mutual correlation (0 ≤ k ≤ k max , 0 ≤ l ≤ l max ): Fig. 5. Procedure for determining semiautomatic cross-sectional lines.

Fig. 8 .Fig. 9 .
Fig. 8. Example of mismatch of template R with window S in image 59

Fig. 11 .
Fig. 11.Examples of the reference values for estimation of the error caused by automatic cross-sectional line positioning (description in the text)

Fig. 12 .
Fig. 12.The changes in the slope angle for the manually positioned cross-sectional line in the T c period: a) the slope angle, b) the line length / (number of pixels -1), c) a diagram of the influence of the cross-sectional line slope with the axis of an artery fragment on the line length in pixels.

Fig. 13 .
Fig. 13.Densitometric curves A(t) and their low pass version for the same cross-sectional line of the catheter calculated for two background line interpolation algorithms: a) method I, b) method II

Fig. 14 .
Fig. 14.The position of the artery cross-sectional line on the reference image: a) correct position (visual estimation), b) incorrect position (intentional) Fig. 15.Ranges of values of the area ∑ ) t ( A for automatic correct, automatic incorrect and manual positioning of the cross-sectional line relationship between the X-ray intensities I and I 0 and the corresponding image brightness D and D 0 (image value without the indicator on the same site for the same exposure parameters , i.e. a background value) respectively, the value V pi of the indicator in the volume along the distance l, viewed as an area of one pixel a [cm 2 ], can be expressed as: