Interpolated Background Subtraction Method For Coronary Stenosis Quantification

The determination of percent stenosis of coronary arteries is an important task in medicine. In this paper we discuss three different algorithms which can be used in conjunction with videodensitometry to measure this quantity. These algorithms may be used in subtracting the background under a vessel segment thus eliminating the need for a preinjection mask. Mathematical details of the algorithms and experimental results are presented.


INTRODUCTION
Recent advances in the digital x-ray video imaging systems has made the visualization and measurement of small amounts of radiopaque contrast agents within the vasculature a reality.The components of a typical digital radiographic system are the x-ray generator, x-ray tube, intensifier (II), television (TV) camera, video digitizer, display and recording devices.Two general classes of measurements one can perform on the digital video imaging are: (a) geometrical^ and (b) videodensitometric quantification^,6 Here, we will only be concerned with the latter.Videodensitometric measurements * are based on the linear relationship between the intensity of the video signal and the thickness of the tissue transversed by x-rays causing that signal.Many physical problems, some of which are intrinsic to such imaging systems and others due to the nature of interaction of x-rays with an object, usually cause a deviation from this desired linear relationship4 .Thus, it is quite essential to determine the response of the imaging chain prior to making measurements on images obtained with such systems.In this paper, we will focus on the videodensitometric measurement of coronary stenosis.The most common technique used for this purpose is the "mask subtraction" 3. in this technique, a mask image of the patient is taken prior to the injection of iodinated contrast material into the coronary vessels.Then, a pixel by pixel subtraction of the logarithmically amplified mask image from the post-injection image results in an image which "ideally" has all the noniodinated parts of the object cancelled.The net iodine signal intensity is linearly related to the thickness of the vessel containing it.Thus, a correct measurement requires, among other things, a total cancellation of the noniodinated overlaying structures.If the patient moves in between the time mask and post-injection images are acquired, then the difference images may be degraded due to motion artifacts.The measurement of coronary stenosis using such images would give rise to different degrees of error depending on the magnitude and location of motion artifacts.Since the physiological motions which occur in a patient are not simple translational motion, one can not remedy the situation by simply shifting the mask to obtain an acceptable registration with the iodinated image.Attempts to alleviate this problem by means of rubber sheet algorithms have proven to be difficult and at best, approximate.In this paper, we discuss several different algorithms to estimate the background under the iodinated vessel.Since the background is derived from the iodinated images, one does not have to use any preinjection mask images.Section 2 covers the methods.The results are presented in Section 3. The conclusions are given in Section 4.

Basic Formalism
The expression for a two dimensional video image Q(x,y) obtained with a TV based x-ray imaging system is given by Kruger2, where a is an amplification factor and y a constant.It should be emphasized that for a dynamic system such as the cardiovascular system, images are also an explicit function of time.In this paper the time dependence will be kept implicit.The function I is related to the transmitted photon fluence by, I(x,y) = a N(x,y) and N is defined by, N(x,y) = P(x,y) + S(x,y) In equation ( 3), P and S are the primary and scattered x -ray distribution components, respectively.The primary x -ray component P is related to the incident one by the following equation.
where No is the incident x -ray flux at location (x,y) and p. the spatially position dependent linear attenuation coefficient along the x -ray beam direction which was assumed to be along the z axis.Note that in writing eq.( 4), we ignored the polychromaticity of the x -ray beam.The approximation of taking the attenuation coefficient at an effective energy <E> appears to be a good approximation when one attempts to determine coronary stenosis by videodensitometry.The scatter component S is related to the incident flux by a complex relationship.In most physical situations, the scatter component S does not contain any high spatial frequencies, i.e. it is a slowly varying function of spatial coordinates (x,y) over the dimensions of the vessel.In equations ( 2) to (4), we also assumed that other physical problems such as veiling glare which are inherent to image intensifier based systems are not significant for the measurement of coronary stenosis even though they cause an overall reduction in image contrast.
Since the measurement of coronary stenosis involves the injection of iodine into the vessels, we will now deal with this situation.The primary signal after the contrast material is injected into the vessels is given by, PI(x,Y) = P(x,y) exp [-f µI(x,Y,z: <E>) dz] (5) where µl is the iodine linear attenuation coefficient measured with respect to blood which was displaced by the contrast agent.From now on, we will implicitly assume that the spatial extent of (x,y) in all the equations is limited to within the boundaries of the vessel where a videodensitometric measurement will be performed.Using the new primary signal PI in eq. ( 3) and taking the logarithm of the video signal Ql results in, where a is a constant and N defined by eq. ( 3) is the background signal under the iodinated vessel.In writing eq. ( 6), it was also assumed that the signal from the iodinated vessel is much smaller than the background signal under it.We will now rewrite eq. ( 6) as, ln QI(x,y) = a ln N(x,y) -R I µI (x,y,z: <E>) dz (7) where we defined R = P(x,y)/N(x,y).The constant R is independent of (x,y) in a first order approximation.The last term in eq. ( 7) can further be simplified by assuming that the contrast material is homogeneously mixed within the vessel.This results in, l' (x,Y) = R µI (< E >) I dz where a is an amplification factor and y a constant.It should be emphasized that for a dynamic system such as the cardiovascular system, images are also an explicit function of time.In this paper the time dependence will be kept implicit.
The function I is related to the transmitted photon fluence by, and N is defined by, N(x,y) = P(x,y) + S(x,y) In equation ( 3), P and S are the primary and scattered x-ray distribution components, respectively.The primary x-ray component P is related to the incident one by the following equation.
where No is the incident x-ray flux at location (x,y) and |i the spatially position dependent linear attenuation coefficient along the x-ray beam direction which was assumed to be along the z axis.Note that in writing eq.( 4), we ignored the polychromaticity of the x-ray beam.The approximation of taking the attenuation coefficient at an effective energy <E> appears to be a good approximation when one attempts to determine coronary stenosis by videodensitometry.The scatter component S is related to the incident flux by a complex relationship.In most physical situations, the scatter component S does not contain any high spatial frequencies, i.e. it is a slowly varying function of spatial coordinates (x,y) over the dimensions of the vessel.In equations ( 2) to (4), we also assumed that other physical problems such as veiling glare which are inherent to image intensifier based systems are not significant for the measurement of coronary stenosis even though they cause an overall reduction in image contrast.
Since the measurement of coronary stenosis involves the injection of iodine into the vessels, we will now deal with this situation.The primary signal after the contrast material is injected into the vessels is given by, Pl(x,y) = P(x,y) exp [ -J w(x,y,z: <E>) dz] (5 where jij is the iodine linear attenuation coefficient measured with respect to blood which was displaced by the contrast agent.From now on, we will implicitly assume that the spatial extent of (x,y) in all the equations is limited to within the boundaries of the vessel where a videodensitometric measurement will be performed.Using the new primary signal PI in eq. ( 3) and taking the logarithm of the video signal Qi results in, where a is a constant and N defined by eq. ( 3) is the background signal under the iodinated vessel.In writing eq. ( 6), it was also assumed that the signal from the iodinated vessel is much smaller than the background signal under it.We will now rewrite eq. ( 6) as, In Qi(x,y) = a In N(x,y)p J m (x,y,z: <E>) dz (7) where we defined P = P(x,y)/N(x,y).The constant P is independent of (x,y) in a first order approximation.The last term in eq. ( 7) can further be simplified by assuming that the contrast material is homogeneously mixed within the vessel.This results in, ,y) = p W (<E>) Jdz ( 8) It should be pointed out that in eq. ( 8), the integral over z is equal to the thickness of vessel L(x,y) in the direction of the x -ray beam at a location (x,y).The first term in eq. ( 7) results from the background under the vessel.A videodensitometric determination of vessel volume requires the subtraction of the background term from eq. ( 7) and integration of the resultant signal over (x,y) within the boundaries of the vessel segment under consideration.In the mask subtraction technique, the background term in eq. ( 7) is eliminated by choosing an image of the region of interest prior to the arrival of iodine.However, since the mask and post-injection images are obtained at different times, there is a great probability that they may not be in perfect registration even if one uses ECG gating.
Let us assume that we have a vessel segment where we wish to measure the cross section videodensitometrically, this is illustrated in Fig. 1.The cross section to be determined is perpendicular to the length of the vessel which was taken to be in the direction of A -D.If the cross section of the vessel does not change appreciably along its length, it proves to be advantageous to measure the vessel volume bounded by A'B'C'D' and dividing it by the height h to obtain an average cross sectional area.The logarithmically amplified video signal from eq. ( 7) can be integrated along the y direction to obtain, 1 (x) _ -J ln QI (x',y') dy' = -a J ln N(x,y) dy + R µl (< E >) J L(x',y) dy' (9) The function r) which is defined in eq. ( 9) is illustrated in Fig. 2. In order to calculate the net video signal coming from the iodine filled vessel (shaded area), we have to subtract the background area within A' <-x' <_ B'.Let us denote the background signal along the x' direction by BG(x').This corresponds to the first term on the right hand side of eq. ( 9).The computation of the cross sectional area is achieved by evaluating B' The double integral term in eq. ( 10) is the volume of the vessel bounded within the area A'B'C'D'.Dividing equation (10) by the length of the vessel h yields 692 / SPIE Vol.914 Medica / Imaging ll (1988) where we moved the attenuation coefficient pi outside the integral sign.It should be pointed out that in eq. ( 8), the integral over z is equal to the thickness of vessel L(x,y) in the direction of the x-ray beam at a location (x,y).The first term in eq. ( 7) results from the background under the vessel.A videodensitometric determination of vessel volume requires the subtraction of the background term from eq. ( 7) and integration of the resultant signal over (x,y) within the boundaries of the vessel segment under consideration.In the mask subtraction technique, the background term in eq. ( 7) is eliminated by choosing an image of the region of interest prior to the arrival of iodine.However, since the mask and post-injection images are obtained at different times, there is a great probability that they may not be in perfect registration even if one uses ECG gating.
Let us assume that we have a vessel segment where we wish to measure the cross section videodensitometrically, this is illustrated in Fig. 1.The cross section to be determined is perpendicular to the length of the vessel which was taken to be in the direction of A-D.If the cross section of the vessel does not change appreciably along its length, it proves to be advantageous to measure the vessel volume bounded by A'B'C'D' and dividing it by the height h to obtain an average cross sectional area.The logarithmically amplified video signal from eq. ( 7) can be integrated along the y' direction to obtain, TI (x') m -J In Qi (x',y') dy' =a J In N(x,y') dy + p m (< E >) J L(x',y') dy' (9) The function TI which is defined in eq. ( 9) is illustrated in Fig. 2. In order to calculate the net video signal coming from the iodine filled vessel (shaded area), we have to subtract the background area within A' < x' < B'.Let us denote the background signal along the x' direction by BG(x').This corresponds to the first term on the right hand side of eq. ( 9).The computation of the cross sectional area is achieved by evaluating B ' The double integral term in eq. ( 10 where A s and 0 p are related to the cross sectional areas of the stenotic and patent vessels, respectively [c.f.eq. ( 11).In the following sections, we will introduce various algorithms for the estimation of background BG(x) under the vessel, i.e. for A'<_x'S B'.

Background Interpolation Algorithms
The problem to be solved is the estimation of the background BG(x) within A' <_ x' <_ B' from a knowledge of its values for x'< A' and x' > B', see Fig. 2.

Linear Interpolation (LIN) Method
In this subsection, we will assume the background to be linear under the vessel, i.e. for A' <_ x' <_ B'.The background is estimated by computing its average values from two adjacent regions outside the vessel.These regions will be taken to be: The average values are given by A' BGL = J BG(x') dx'/(A'-A) A B BGR =.r BG(x') dx'/(B-B')

B'
The estimated background BGEi is defined by The videodensitometric cross sectional area is computed by where <A> is the average cross sectional area of the vessel along the segment under consideration.Percent area stenosis is computed by, where A s and A p are related to the cross sectional areas of the stenotic and patent vessels, respectively [c.f.eq. ( 11).In the following sections, we will introduce various algorithms for the estimation of background BG(x') under the vessel, i.e. for A' < x1 < B'.

Background Interpolation Algorithms
The problem to be solved is the estimation of the background BG(x') within A1 < xf < B' from a knowledge of its values for x'< A! and x1 > Bf, see Fig. 2.

Linear Interpolation (LDsH Method
In this subsection, we will assume the background to be linear under the vessel, i.e. for A1 < x1 < B'.The background is estimated by computing its average values from two adjacent regions outside the vessel.These regions will be taken to be: which is similar to eq. (11).

Least Mean Squared Error (LMSE) Method
Let us assume that the background can be approximated as a polynomial given by M BGE2(x) = E ak (x')k k=0 The mean error between the actual background BG(x') and the estimated one is minimized by computing, (aQ/aak) =0 (19) where In equation ( 20), the summation over j is performed for values of x' given in eq. ( 13), i.e. regions adjacent to but outside the vessel boundaries.A further simplification takes place if x' =0 is chosen to be halfway between A' and B'.In this case, all the odd powers from the summation in eq. ( 18) drop out due to asymmetry and we are left with a smaller array to deal with.Computation of the coefficients ak is achieved by using eqs.(18 -20) and using a matrix inversion to obtain ak.The net iodine signal is given by, B' The procedure can be explained with the aid of Figure 2. First, we perform a separate mean least squares fit to the background in the two regions using: Once we determine ak and ck according to the procedure given in Section 2.2.2, we then know both of these background functions.The cubic spline process involves fitting a third order polynomial BGE3(x) within A' S x' _< B' in such a way that the following conditions hold:

Least Mean Squared Error fLMSE^ Method
Let us assume that the background can be approximated as a polynomial given by M BGE2(x') = I aktx'^ (18) k=0 The mean error between the actual background BG(x') and the estimated one is minimized by computing, OQ/3ak ) = 0 (19) where

20) j=l
In equation ( 20), the summation over j is performed for values of x* given in eq. ( 13), i.e. regions adjacent to but outside the vessel boundaries.A further simplification takes place if x =0 is chosen to be halfway between A1 and B!.In this case, all the odd powers from the summation in eq. ( 18) drop out due to asymmetry and we are left with a smaller array to deal with.Computation of the coefficients ak is achieved by using eqs.(18-20) and using a matrix inversion to obtain ak.The net iodine signal is given by, B 1 A2 = J[Tl(x')-BGE2(x')]dxyh (21) A1

Cubic Spline Interpolation (CSD Method
The procedure can be explained with the aid of Figure 2. First, we perform a separate mean least squares fit to the background in the two regions using: Once we determine ak and ck according to the procedure given in Section 2.2.2, we then know both of these background functions.The cubic spline process involves fitting a third order polynomial BGE3(x*) within A' < x' < B f in such a way that the following conditions hold: we can determine Ai for 0 <_ i <_ 3 using the conditions given in eq. ( 23) since we already know bL and bR from eq. ( 22).
The net iodine signal is determined by A3 evaluated as B' 3. RESULTS

System Linearity (25)
In order to test the linearity of the system in measuring volumes, we performed a series of experiments in which a balloon was filled with different amounts of 10% Renografin -76.The volume of the balloon was determined videodensitometrically after logarithmic mask subtraction by integrating the video signal within the boundaries of the balloon.The results obtained by changing the balloon volume are shown in Fig. 3.The vertical axis is the videodensitometrically measured volumes in cubic centimeters.The actual volumes of the balloon are along the horizontal axis.The normalization factor for the actual volumes along the y -axis was obtained by using the straight line which was computed using the linear least squares fitting.The images were obtained at 70 kVp with 4 -mm aluminum filtration and an x -ray grid was placed on the II.The linearity of the data is quite impressive.
we can determine Ai for 0 < i < 3 using the conditions given in eq. ( 23) since we already know bL and bR from eq. ( 22).The net iodine signal is determined by A3 evaluated as In order to test the linearity of the system in measuring volumes, we performed a series of experiments in which a balloon was filled with different amounts of 10% Renografin-76.The volume of the balloon was determined videodensitometrically after logarithmic mask subtraction by integrating the video signal within the boundaries of the balloon.The results obtained by changing the balloon volume are shown in Fig. 3.The vertical axis is the videodensitometrically measured volumes in cubic centimeters.The actual volumes of the balloon are along the horizontal axis.The normalization factor for the actual volumes along the y-axis was obtained by using the straight line which was computed using the linear least squares fitting.The images were obtained at 70 kVp with 4-mm aluminum filtration and an x-ray grid was placed on the II.The linearity of the data is quite impressive.

Phantom Studies
In these studies, a phantom consisting of a Plexiglas block containing holes was used.An illustration of the phantom is seen in Figure 4.The diameter of the patent segment was 4 -mm.The diameter of the stenotic segment was 0.51, 0.71, 0.99, 1.51 and 1.99 (all in mm) in five different blocks.The holes were filled with 50% Cysto -Conray II (81 mg/ml organically bound iodine).The x -ray tube was operated at 75 kVp with a 4 -mm aluminum added filter.In the first case, five Plexiglas blocks were placed behind 10 cm thickness of plexiglas to provide a uniform background.A second study was performed in which the same blocks were placed within a chest phantom one by one to simulate a more realistic nonuniform background.The x -ray images were digitized into a 512X512 format using an 8 -bit analog -to-digital converter (ADC) and stored in a video RAM for image analysis.

Linear Background Interpolation (LIN)
Five images with the uniform background were analyzed using the linear background interpolation algorithm given in Section 2.2.1.A similar calculation was repeated on the second set of five images which included the chest phantom as the background.The net video signal was computed using eq.( 17) for the patent and stenotic regions for all images.The percent area stenosis was calculated using eq.( 12).The results comparing the actual and videodensitometrically measured percent area stenosis are plotted in Fig. 5.The straight lines shown in the figure were obtained by a linear least squares fit.The data displayed by open squares was obtained with the uniform background phantom.The data arising from the chest phantom is shown with solid triangles.We observe that with the linear background interpolation method (LIN) the videodensitometric measurement of percent stenosis for both phantoms correlates well with the actual values.The LMSE algorithm described in Section 2.2.2 was used in analyzing the same images as in the previous subsection.The order of the polynomial in eq. ( 18) was taken to be 4.The graphs of the measured videodensitometric percent stenosis for the two sets of images are shown in Fig. 6.The correlation of the videodensitometric data with the actual percent stenosis measurements seems quite satisfactory.

Phantom Studies
In these studies, a phantom consisting of a Plexiglas block containing holes was used.An illustration of the phantom is seen in Figure 4.The diameter of the patent segment was 4-mm.The diameter of the stenotic segment was 0.51,0.71,0.99, 1.51 and 1.99 (all in mm) in five different blocks.The holes were filled with 50% Cysto-Conray II (81 mg/ml organically bound iodine).The x-ray tube was operated at 75 kVp with a 4-mm aluminum added filter.In the first case, five Plexiglas blocks were placed behind 10 cm thickness of plexiglas to provide a uniform background.A second study was performed in which the same blocks were placed within a chest phantom one by one to simulate a more realistic nonuniform background.The x-ray images were digitized into a 512X512 format using an 8-bit analog-to-digital converter (ADC) and stored in a video RAM for image analysis.

Linear Background Interpolation (LIK)
Five images with the uniform background were analyzed using the linear background interpolation algorithm given in Section 2.2.1.A similar calculation was repeated on the second set of five images which included the chest phantom as the background.The net video signal was computed using eq.( 17) for the patent and stenotic regions for all images.The percent area stenosis was calculated using eq.( 12).The results comparing the actual and videodensitometrically measured percent area stenosis are plotted in Fig. 5.The straight lines shown in the figure were obtained by a linear least squares fit.The data displayed by open squares was obtained with the uniform background phantom.The data arising from the chest phantom is shown with solid triangles.We observe that with the linear background interpolation method (LIN) the videodensitometric measurement of percent stenosis for both phantoms correlates well with the actual values.The LMSE algorithm described in Section 2.2.2 was used in analyzing the same images as in the previous subsection.The order of the polynomial in eq. ( 18) was taken to be 4.The graphs of the measured videodensitometric percent stenosis for the two sets of images are shown in Fig. 6.The correlation of the videodensitometric data with the actual percent stenosis measurements seems quite satisfactory.In this case, we applied the CSI algorithm described in Section 2.2.3 to compute the percent stenosis using the same images as discussed before.A plot of videodensitometric CSI measurements versus the actual ones is shown in Fig. 7.The correlation of measured data with the actual ones again is quite well.In this case, we applied the CSI algorithm described in Section 2.2.3 to compute the percent stenosis using the same images as discussed before.A plot of videodensitometric CSI measurements versus the actual ones is shown in Fig. 7.The correlation of measured data with the actual ones again is quite well.the vessel varies rapidly within the dimension of the vessel, linear background subtraction may be inferior to the LMSE or CSI techniques.We have also seen from the presented results that the physical problems such as veiling glare, beam hardening and x -ray scatter do not affect the measurement of percent stenosis.All the computations presented here can be performed quite rapidly in a video image processor.The main advantage of the three techniques discussed in this paper is that one does not have to perform subtraction of any mask obtained prior to contrast material injection thus avoiding any possible misregistration artifacts.the vessel varies rapidly within the dimension of the vessel, linear background subtraction may be inferior to the LMSE or CSI techniques.We have also seen from the presented results that the physical problems such as veiling glare, beam hardening and x-ray scatter do not affect the measurement of percent stenosis.All the computations presented here can be performed quite rapidly in a video image processor.The main advantage of the three techniques discussed in this paper is that one does not have to perform subtraction of any mask obtained prior to contrast material injection thus avoiding any possible misregistration artifacts.

Figure 1 .Figure 2 .
Figure 1.Selection of a region -of-interest (ROI) bounded by ABCD along a vessel bounded by A'B'C'D'.

Figure 1 .
Figure 1.Selection of a region-of-interest (ROI) bounded by ABCD along a vessel bounded by A'B'C'D*.

Figure 2 .
Figure 2. Integrated videodensity profile T] as a function of distance x' perpendicular to the vessel edges.

Figure 3 .
Figure 3. Videodensitometrically measured balloon volume vs. actual balloon volume, all in cc's.The best fit line is also shown.

Figure 6 .
Figure 7. CSI interpolation Figure 7. CSI interpolation supported in part by PHS Grant No. HL31440 awarded by the National Heart, Lung and Blood Institute, DHHS.The authors also acknowledge a research grant awarded by Philips Medical Systems, Inc.