- Email: [email protected]
Digital curvature estimation for left ventricular shape analysis
Digital curvature estimation for left ventricular shape analysis Maurizio Baroni and Giuseppe
Curvature estimation from noisy digital curves is addressed for shape analysis of angiographic left ventricle (LV) images. A simple and fast algorithm, based on Fourier series approximation, changes the harmonic number and filter window along the closed (or partially open) contours to optimize smoothness and reconstruction errors. The curvatures estimated on a set of ventricle-shaped models are significantly closer to their analytical values than those assessed by four other methods. Preliminary results on LV data are presented. Keywords: left ventricular shape, curvature estimation, adaptive filtering
INTRODUCTION An important goal in image processing is shape description and recognition. As a result of the imaging process, a 3D object is mapped onto a bidimensional projection. This image is digitized with a square sampling grid and segmented into a bilevel picture or a boundary is extracted from it. Both region and boundary representations define the ‘shape’ of the object on a background, and in most situations still carry the main information needed for recognition purposes. Even with careful acquisition and segmentation procedures, various sources of errors are possible; moreover, both projection and sampling operators are non-invertible: therefore, digital shape analysis is not a simple problem, even if we restrict ourselves to planar shapes. Many papers have been written on this subject, and excellent reviews can be found’,2. Shape analysis
Barletta”
methods can be classified in series (global) or piecewise (local) representations, internal (area-based) or external (boundary-based). These latter can further be classified with regard to the boundary parameterization (rectangular, radial, angular or curvature parametrization). For example, the method of moments is global and internal; Fourier transform is a global method which can be internal (if applied in 2D) or external, using any parametrization of the boundary. Examples of local methods are Medial Axis Transformation (internal) and polygonal approximations of the boundary. Our attention is concentrated on boundary representations of simply connected curves, because we are interested in the study of left ventricle (LV) shape as it appears in standard cardiac imaging (contrast angiography, ecography). Endocardial shape, in particular, plays an important role in ventricular function’, and its analysis provides a promising tool for clinical and pathophysiological investigations. For this purpose, the curvatrure representation is chosen because it satisfies all the following criteria:
Dipartimento di Tngegneria Elettronica, via S Marta 3, 50130 Firenze, Italy ‘Institute of Internal Medicine, Cardiology Unit, University of Florence, Florence, Italy
1. Curvature is invariant under rotation and translation of the curve, and can easily be normalized with respect to scale changes. as it is possible to 2. It is a unique representation, match the curvature along a curve against that of homogeneous populations of curves. after suitable normalizations. able to carrying informa3. It is a local representation tion at varying resolutions, and to provide a structural description in terms of primitives and relations. 4. It can represent both closed and open curves, even in the case of missing data. to small changes in part of 5. It is relatively insensitive the curve and, theoretically, no arbitrary choice is needed for its computation.
Paper received: 5 April 1991; revked papers received: 1991
However,
0262-8856/92/007485vol IO no 7september
1992
2 September
10
0
the curvature
1992 Butterworth-Heinemann
computation
involves
first and
Ltd 485
second order differentiation of noisy digital data, so smoothing is necessary to obtain significant values. Polynomial fitting and standard spline interpolation4 are highly sensitive to noise as the equivalent filtering bandwidth is suboptimal. Regularization theory with spline functions” or Fourier series’ is characterized by high performance at the cost of long computation times. As an alternative, lowpass filtering may be optimized by estimating the signal bandwidth’. In this paper a new fast algorithm for digital curvature computation is proposed and compared with another four methods. Three of them are based on the orientation-length function which is differentiated by using either differences taken over a window, or FIR filtering, with a fixed number of points or a fixed sampling step. The fourth method consists of Gaussian filtering of contour coordinates. Our algorithm is a windowed Fourier series approximation of contours, in which the number of harmonics and filter-window are locally chosen to minimize the reconstruction errors, and to maximize the smoothness of the curve. The comparison was carried out on analytically derived curvatures of a set of ventricle-shaped models, using both quantized software-generated curves (here called shape models) and manually digitized contours of the same shapes (test contours). Moreover, preliminary results on LV data are also presented.
METHODS Curvature
and noise
The curvature of a planar curve at a point P is defined as the reciprocal of the radius of the osculating circle (the circle tangent to the curve at P with the highest degree of contact): K=
l/R
(I)
From differential geometry, the curvature can be expressed as the derivative of the unitary tangent vector with respect to arc length 1, measured along the curve from a starting point. Therefore, if we let Q(f) be the slope angle between the local tangent line and some arbitrary x-axis, the (scalar) curvature, k(f), is given by: K(1) = d@ldl If the curve is expressed can write: d@ -=-.dl
d@
dx
dx
dl
It is known
(2) in explicit
form, y = y(x),
that
dl= (dx2 + dy2)“2 = dx[(l + (dyldx)2]“2, dyldx, then: d@ -= dx
one
d2yldx2 1+ (dyldx)2
and
tg@ =
It follows
that: d2yldx2
K(x) = [l + (dyldx)2]“2 Since y(x) is a multivalued function, a parametric form x =x(t) and y = y(t) is more suitable if the parameter t varies monotonically along the curve. The simplest choice is t = 1, although t might be equal to the angle defined by radius vectors projected between the curve and its centroid. This choice, however, leads again to multivalued functions for non-convex shapes, unless one resorts to some artifice’. Now, letting ’ and ’ be the first and second derivatives with respect to t (i.e. x’ = dxldt, x” = d;l d;, etc.), one can write: dy dy ~=~‘~=x’
dt
y’
and d2y -= dx2
Finally,
equation
d(y’/x’) p= dx
x’y”-y’x” X
I2
.-
1 X’
(3) becomes
x’y”_y’x”
K(t) =
(4)
(x12+yt2)3/2
The curvature sign depends on vector direction: for closed curves it is conventionally taken as positive for the convex segments and negative for the concave ones. In discrete terms, when a collection of points is available, rather than a continuous function, equation (1) is implemented either using the cocircularity concept’ or with a simple geometric construction based on triplets of pixels4. The angles in equation (2) are locally estimated, and the derivatives in equations (2,3, 4) may be replaced by first differences. However the discrete curvature computation faces the problem of noise, given by both quantization and statistical errors. Here, the signals underlying discrete noisy data are assumed to be differentiable, and are therefore certainly smooth and hence practically band-limited as their spatial frequency spectra will be concentrated at low frequencies. In general, given a signal with limited bandwidth (B), sampled at equispaced steps (llfs), affected by additive stationary white noise (~3, the uncertainty in the computation of the mth order derivative’” is given by: &CT;.
B(2m+‘)l[fs* (2m + l)]
(5)
Therefore, the sampling frequency must be chosen not only with respect to signal bandwidth, but also in relation to the noise level and derivative order. In the one-dimensional case, the samples are also quantized with M bits. This is equivalent to an additive white noise” with variance given by: a; = 2-2M/12
(6)
In two dimensions, the sampling requirements to preserve shape are more restrictive, and if d is the linear dimension of the finest details, the Nyquist image and vision computing
condition
Shape models and test contours
(llfs) d (d/2) becomes4:
(llfs) d (d/W2)
Since the curve x(t), y(t) is sampled have:
on a square grid we
yi = int b(i/fJ
x, = int [x(i/fy) + 0.51,
+ OS]
and successive slope angles can differ only by a multiple of 45”. These abrupt changes in direction are amplified by differentiation and give rise to a still coarser quantization in curvature. Even the convexity property of the curve may be locally violated. This means that a certain order of smoothness must be imposed to digital curves. However the scale at which the smoothing is performed strongly influences the result: with over-smoothing small scale details are lost, while under-smoothing leaves undesirable noise fluctuations. Another problem is related to the non-uniform sampling of the length: in fact, the distance between two consecutive pixels is 1 along the sampling grid directions and V? in diagonal directions. For this reason, the digital contour must be resampled by linear interpolation with spacing equal to PIN, where N is the number of equispaced pixels, and the perimeter P is estimated along the original curve’2,‘3.
Shape models
The LV silhouttes seen in standard cardiac imaging are smooth curves with high and low, positive and negative More particularly, the LV contours curvatures. obtained in contrast angiography with RAO projection are characterized by a highly variable valvular zone and two low curvature regions separated by a high curvature apical zone. For the purpose of obtaining a set of ventricle-shaped closed curves, the parametric equations of the ellipse were generalized to allow accumulation of high order terms, just like a truncated Fourier series: x(t) = x0 + Ji,,A ,,I sin (nt + F,,,). y(t) = yo+ &,Avnsin
(nt+
F,.,,), n=l,2,3,1~[0,2rr],
(7a)
where (X,,, Y,) is the central point, the number of terms was limited to 3, and the other constants were chosen by means of an interactive graphical procedure. Table 1 lists the constant values chosen for 39 models, and some of them are showed in Figure 1. Curvature can be analytically computed according to equation (4) by differentiating equations (7a), i.e. for x(t): x’ (t) = 2;,, nA i,, cos (nt + F,,) , x”(t) = -Ii,,
n2 A,,, sin (nt + F,,,),
(7b)
Table 1. Model coeffkients for
Constant values A,,
A,.,
Fx,
Fw
AC;_7AYZ Friz FYJ
+60
+50
+40
+90
+20
0
+70
+ 80
+40
+90
+20
0
+70
+20
-t 90
+20
+80
+20
Ax?
A,,
F,.I Fy.7
0
0
+20
0
+50
+70
0
0
+20
0
+50
0
+70
0
0
+20
0
+50
0
+70
0
0
+20
0
+75
+70
0
0
+20
0
+50
39 ventricle-shaped curves
i + 80
4 5
I
f60
I
IO
+I10
II
+80
I
I
+40
I5 16
+x0
+70
+40
I
I +100
20 21
I 26 27 I 32 33
I
39
+ 80
+70
+40
-t 100
+20
+5
+30 +70
+70
+40
+90
+20
0
+70
0
0
+10
0
+50
+50
+40
+110
+20
0
+70
0
0
+5
0
+70
I
+95 +60
t i +80
WI IO no 7september
1992
487
The orientation-length function Q(I) can be approximated by taking either cosine or tangent measures on successive points, as reviewed by Teh and Chint4. Here equation (2) was rewritten as:
(8) where {Xi, YJi =:1 . . . N} are the coordinates of the digitized contour, and w B 1 is an integer parameter which defines the region of support (window) for the computation of differences. Choosing a w> 1 window has the effect of improving the direction quantization in a way that is similar to high order chain codes”.The angular quantization step can be written as: Figure I. Four sets of ventricle-shaded curves (3rd, Sth, 6th and 7th set in Table It Therefore, these shape models provide a gold standard to compare the methods for curvature estimation described in the next sections. They were used in two ways: 0
0
the coordinates generated by equations (7a) were sampled at steps of 2”, rounded to integer values, and directly given as input to the analysis programs; computer displays of shape models were printed on paper and manually digitized by two experienced observers using a 0.1 mm resolution graphic tablet, according to the usual procedure for digitizing actual LV data. The discrete coordinates of these test contours were resampied at 1 mm equispaced steps, and analysed by the same programs.
The former shape models generate a set of curves affected only by quantization noise (SNR, ~25 dB), while the latter test contours represent the worst case in which the statistical errors due to manual tracing are also present. Moreover, the points sampled using equation (7a) are equispaced in angle not in arc length, so that they were resampied at equidistant steps by an iterative procedure to allow a point-by-point comparison between the analytical curvatures and those estimated from test contours. Finally, since the number of points varies with contour size, the comparison was carried out by linearly interpolating ail the curvature values at N equidistant points. This size normalization implicitly assumes that equal length abnormalities are more significant in smaller ventricles. A value of N = 90 was chosen because the analytical curvature sampled at a lower number of points underscores the sharper peaks seen along our shape models.
Window-based method Now the methods
introduced. 488
used for curvature
estimation
are
@i+w-
and
Q,i_,,=2?T/ZM
M=w+2
(9)
Of course, the window size must be large compared with the spatial noise frequency so that it was chosen by the following procedure. Since the curvature function along a closed curve is periodic, Fourier analysis can be properly applied. The amplitude spectra of the curvatures, ana~yticaiiy computed on the shape models with those previously described, were compared obtained on the test contours using equation (8) and varying w in the range [2, 201. The best signal-to-noise ratio was achieved with a window size of w = 12, using 1 mm resolution.
Filtering the orientation function In the previous window method the curvature at each point depends on the orientation value in only two distant points, neglecting ail the others. This does not ensure the elimination of spurious oscillations, so a further smoothing of curvature values is needed. A better use of the available data is obtained by averaging all the orientation values over the window, or filtering them to reduce high frequency noise. For this purpose it is necessary to define a cumulative orientation-length function which is continuous along the contour. In fact, since the arctangent function returns angles between -7r and +n-, the slope angle function Q(j) has artificial discontinuities, i.e. for a closed curve Q(O) -Q(L) = 27r. At these points an offset of fn= was added, and the arc length IE [O, L], was replaced by the parameter te[O, 27r]ih. The resulting angular function: B(t) = @(tL/27)
- Q(O) + t
(10)
is periodic, and can be expanded in discrete Fourier series. Its power spectrum j O(f) 1’ was used to select the filter bandwidth. If we assume a white noise, and that f,22B/c (c = 0.8), the noise power No can be assessed by averaging the spectrum between c. fJ2 and f,/2. Now the cutoff frequency fc where the signai-tonoise-ratio reaches a threshold (KT> 1) can be chosen: i@(fc)?=
Kr+No
(11)
A FIR filter was then designed to reduce the frequency image and vision computing
Having chosen N = 90 for our ventricle-shaped tours and estimated NC = 12 by means of equation we have f,., = 12/90 = 0.1333.
con(1 l),
Filtering contour data
-50
1
0
0.1
I
0.2
0.3
0.4
0.5
Norml~r~dfwqumcy
Figure 2. Frequency response of digital filter used to differentiate slope angle data: f,, = 0.1333, CY= 3.39 with f >fC and, at the same time, to all the others, according to equation (2). response of a differentiator taken at an of samples” is given by
components differentiate The impulse odd number ho = 0.
hk = [2z-kf,, cos(2nkf,.,) - sin (2rkfC,)]/7rk2
where the cutoff frequency was normalized sampling frequency to provide a more general
(12) by the design:
(13) The filter samples must be truncated with a well shaped window to obtain a finite number of coefficients kE[--/,, &I. We used the Kaiser window” with Kh = 5 and shape factor (Y= 3.39, which offers a good compromise between sharp transition band and residual passband oscillations (see Figure 2). This filter was applied to the angle data estimated with equations (10) and (8) using w = 1 and two different sampling schemes. In the first, the test contours were sampled with a fixed frequency so a different number of points is obtained for different curves. In the second, the number N of samples is taken as fixed, and the spatial frequency is variable. In the latter case, the curvature values do not depend on size, while in the former, as well as in all the other methods discussed here, it is necessary to normalize the shape representation for size variations by multiplying the curvatures by the perimetral length. For comparison purpose, the curvatures computed with the fixed N were normalized by an average perimeter. Another difference is in the value of fCn. This is given by equation (13) once fixedf, (1 mm-‘) and estimated fC by equation (11) (fCn ~0.05). On the other hand, it is convenient to choose the fixed value of N as small as possible to reduce the quantization step (9) until it reaches the minimum value given by the sampling theorem. By writing fC = N,. f. as a multiple of the fundamental frequency f0 = 1lL = f,lN, equation (13) becomes: L.,, = NJN ~~1110 no 7 september 1992
(14)
Another common method for differentiating noisy data involves Gaussian filtering, according to Marr’s theory. Using the rectangular parametrization of the contours the x(t) and y(t) functions can be smoothed independently by convolving them with a Gaussian kernel of width CJ, G,(t) = exp(-t’/2a’)/a(2~)‘“. To compute curvature the smoothed x5(t) function must be differentiated twice, and letting * be the convolution integral, x,y’(t) = [x(t) :kG,,(t)]’ = x(t) * G,‘(t). x,“(t) =x(t) * G,“(r), where G,,‘(t) = [-tV~~(27r)“1)] G,,“(t) = [W’(27r)“‘]
exp( -tt’/2rr?) exp(-t’/2u’)(
1 - f’/a*)
The same holds for ys(t). The Gaussian function realizes the best compromise between frequency bandwidth and spatial window, according to the indetermination principle. Moreover, it can be easily truncated at a finite number of coefficients where it falls below a very low threshold. Finally, it has the unique property of not introducing extra zero-crossings of the second derivative’. Therefore, since the zerocrossings indicate the dominant points along the curve (the inflection points, in our case), the Gaussian function was used for multiscale filteringlx. In this work, however, the width value v was chosen proportionally to the number N of points, depending upon noise level and shape details.
Improving
Fourier series approximation
Starting again from a rectangular parametrization of closed contours, the x(t) and y(t) functions are periodic and can be independently expanded into a Fourier series. A truncated expansion can be used to reconstruct a smoothed and/or interpolated version of the contour and, at the same time, to compute curvature according to equations (4) and (7b). The trigonometric expansion of x(t) and y(t) can be expressed in a matrix form, and geometrically interpreted as a sum of ellipses rotating with a speed proportional to their harmonic number”. Similarly, the coordinate functions can be considered as real and imaginary parts of a complex discrete function z(n): z(n) =x(n)
+ jy(n) = Ck=-Nc,Nc hkck e’k(?T’N)n
(15)
where the complex Fourier coefficients are computed via FFT with N = 256 in the usual form: ck
=
(l/N)
&+,,Z,e-~“(2”‘N)k
The weighting factors hk are introduced reduce the spurious oscillations given by any of a trigonometric series. For this purpose, window was used with a shape factor
(16) in (15) to truncation the Kaiser cy chosen 489
proportionally to NC, starting from 3.39. (This filter is the dual form of that used for the orientation function.) As pointed out by Persoon and King-SunzO, any truncated expansion (15) converges to a closed curve, but the reconstructed points are not always equispaced in arc length, so they must be resampled during the same program loop implementing equations (15) and (4). Our concern is now to select the harmonic number NC. Using the amplitude spectra (16) averaged over all test contours does not ensure the best performance in each case. In fact, rewriting equation (ll), it follows that Ic(?Nc)/ 2~ K,N, is satisfied for KT = 20 at N,.=+12, for Kr=6 at NC=+24 and for Kr=2 at NC= +36. (T a k’mg into account the fact that the spectrum of the complex function z(n) is not symmetric, unlike the real case where ck * = c-k, we have considered an average of both the positive and negative harmonic number.) If a low value is chosen for NC (in our context N,.d24), the reconstruction errors e, defined as the Euclidean distances between the original and reconstructed points, are much greater for some test contours than for the smoother ones. In other words, high curvatures are underestimated (see Figure 3a). On the other hand, by increasing NC the errors are reduced on the average, but at the cost of leaving spurious oscillations (see Figure 3b). Further filtering, such as increasing the Kaiser window parameter LY, has an effect similar to that of decreasing NC. A possible approach consists of using a high number of harmonics only for those contours in which a low NC value gives a maximum error above some threshold. Nevertheless, the choice of this threshold may be critical. Another approach stems from the following empirical observation. In the reconstructed test contours the most unrealistic oscillations are apparent in smoother regions, where higher harmonics are redundant (see Figure 3b). Therefore, it seems more reasonable to vary NC locally along each contour (see Figure 3c), instead of choosing NC from curve to curve. A problem now arises to avoid any discontinuity where the value of NC is changed (subsequently referred as transition points). As a first solution, the following simple algorithm was tested (FOURry method): 1. compute FFI using N = 256 points resampled at equidistant steps along the contour perimeter (Equation 16); 256 points using equation (15) with NC = 2. reconstruct 12, calculating the euclidean errors e,, their peak value emax, mean value pu,, and standard deviation cr,? are labelled as 3. the points with e,,
a
-4od b
C Figure 3. (a) Normalized curvatures estimated using N, = 12 and (Y = 3.39 against those analytically computed for test contour no. 2; (b) as 3a using N, = 36 and (Y = 9; (c) results by FOURxy method using same data as 3a, where e,
Both steps (4) and (5) accomplish a sort of morphological filter: removing isolated good points is equivalent to an erosion operation, and precedes the dilation operation to prevent macroscopic discontinuities where the image and vision computing
smoothed contour comes across the original one. However, such discontinuties may be apparent in any point of the transition regions. So it is better to restrict these transition points only between runs of good and bad points. To this aim, the sixth step in the earlier algorithm was removed, and the curvature values at each transition point, as well as in its two neighbouring points, were linearly interpolated. In conclusion, it is worth noting that this method, though based on discrete Fourier series, is not only suited for closed curves. In fact, the other methods based on filtering require some extrapolation of open curves, until half the filter window is beyond the endpoints. Whenever a similar extrapolation could be accomplished for closing the contour, Fourier series can be employed to compute curvature also for open curves. In the angiographic LV contours the aortic plane is usually closed by a straight segment. This may introduce one (or two) corners not existing in the endocardial wall. By extrapolating the contour beyond the valvular margins and then closing it with a straight segment, curvature can be properly measured along the entire LV silhouette. Among the available techniques we used the followextrapolation: x(n + i) = 2x(n) ing symmetrical x(n-i), where i=l, 5 for n=Nand i=-I, -5 for II = 1. The same holds for y.
RESULTS Validation
of methods
Several tests were undertaken to validate the five methods previously described, using both shape models and test contours. Test programs were written in Pascal and run on a 80386 computer (25 MHz, 80387 coprocessor) with execution times ranging from 2s for the window method to 15s for the FOURxy method. First, curvature invariance with rotation, translation and magnification was verified. In particular, the curvature values computed before and after 60” and 90 rotation were found to be identical up to the fifth decimal digit. It is well known that curvature depends on size, the curvature of a circle with radius R being twice the curvature of a circle with radius 2R. So the curvature values, expressed in mm-‘, were multiplied by the perimetral length. This normalized curvature is dimensionless, and for any circle is conventionally equal to 27~. After scaling the test contours by factors of 0.5 and 2, the changes in normalized curvature were less than 1%. Then the accuracy of the curvature estimated by the five methods from discrete noisy rectangular coordinates were investigated with respect to their analytical values. The comparison was accomplished by linear regression and correlation analysis. Table 2 shows the results of this comparison for the set of 39 ventricle-shaped models listed in Table 1. Similar results were achieved for circular and elliptic curves. Obviously, our algorithm (FOURxy) gives curvature values identical to the analytical ones when using
Table 2. Curvature estimated on shape models by five methods against analytical values. (WIN12 = Window method with w = 12); FIRP-FP or FIRP-FS = Filtering of the orientation function using N = 90 (Fixed Points) or f, = 1 mm-’ (Fixed Sampling); FOURXY = Fourier approximation with N, = 12, LY= 3.39 where e, <2pu,; N, = 36, (Y= 9 elsewhere; GAUSS = Gaussian filtering with u = 6 Method
Correlation
Intercept
Slope
SEE
WIN12 FIRP-FS FIRP-FP FOURXY GAUSS
0.5!8 0.627 0.629 0.998 O.Y46
5.68 1.Y6 2.02 0.28 I.51
0.36 0.20 0.21 O.Y4 0.68
x.22 3.41 3.61 0.82 3.1’)
NC = 3. However, even if the harmonic number is not matched to that of the generating model (equation (7a)), the results are better (R>0.99, standard error of estimate - S.E.E. - ~1) than those of all the other methods. In addition, noise oscillations are also minimal: in particular, the curvatures estimated by Gaussian filtering with CT= 2 are very close to the true values (R = 0.992, slope = 0.98) but fluctuate about them. This can be seen by comparing the derivatives of the estimated and true curvature-length functions. Table 3 shows the results of the comparison for the set of 39 test contours. After manual digitization the errors increase. Also the spatial frequency components of the rectangular coordinate signals increase, so that the reconstruction errors are acceptable (order of magnitude 1 pixel) only for N,.z 12. An assessment of the observer variability is given in Figure 4. Following the t-test for paired data, the mean interobserver difference was found to be about 1%) whereas the mean deviation of the FOURxy method from the analytical data was about 10%. According to correlation analysis, FOURxy showed R > 0.90. On the other hand, Gaussian filtering with g E [2, 61 gives R = 0.70, with a low reproducibility for ~<4 due to residual oscillations. The window-based method was the worst (R <0.5 with w = 12): moreover, the resulting curvature-length function exhibited high frequency noise. The method based on filtering of the orientation Table 3. Curvature estimated on test contours by five methods against analytical values. (WIN12 = Window method with w = 12; FIRP-FP or FIRP-FS = Filtering of the orientation function using N = 90 (Fixed Points) or f, = 1 mm-’ (Fixed Sampling); FOURXY = Fourier approximation with N, = 12, cy = 3.39 where e, < 2,~~; N, = 36, LY= 9 elsewhere; GAUSS = Gaussian filtering with u = 6)
WIN12 FIRP-FS FIRP-FP FOURXY GAUSS
OBSI 0BS2 OBSI OBS2 OBSI OBS2 OBSI OBS2 OBSl OBS2
Correlation
Intercept
Slope
SEE
0.290 0.470 0.314 0.315 0.820 0.821 0.919 0.907 0.704 0.723
5.45 4.67 4.80 4.81 2.20 2.1’) I.25 1.32 2.53 2.57
0.22 0.37 0.08 0.08 0.82 0.81 0.88 0.85 0.63 0.62
22.x 22.0 7.33 7.31 6.32 6.29 4. I8 4.39 7.04 6.60
80
60
I
40 N t
Em
d
0
-L 1 ,
10
-20
0
20 Observer
I
I
40
60
60
70 80 mfer. wall
90
80
function gave different results for fixed N (R = 0.31) with respect to taking a fixed sampling step (R = 0.82), but the latter method showed a worse performance in smoothness and variability. In conclusion, the FOURxy method outperforms the others, in terms of both reconstruction errors and noise oscillations, in estimating curvature of ventricle-shaped curves. Moreover, it is relatively insensitive to significant changes in the parameter values. The error threshold eT was changed from pu,+cP to 2,uu,, and also a fixed value was used (er = 2), without significant differences in the estimated curvature. The Kaiser window parameter CYhas no critical value, as well as the morphological filter parameter NT, if limited to few units. Some minor differences were found by changing the low harmonic number (i.e. from 12 to 18), since some test contours have very high curvature variations. The same differences were found by using a lower value for er. On the contrary, no significant difference was found varying the high harmonic number in the range [30, 481. to contrast ventriculography
Our algorithm was applied to the X-ray LV cineangiograms in RAO projection. The LV contours of 16 patients who were found to have normal cardiac anatomy and function were digitized frame-to-frame from 35 mm films, using a moviola and a graphic tablet. Similarly, the LV contours of about 40 patients with aneurysm following anterior infarct were digitized at the end of systolic and diastolic phases. The endsystolic curvature at 90 equidistant points along the contours, averaged on normal and pathological groups, are shown in Figure 5. (The LV contours are traced in the clockwise direction, starting from the anterior aortic margin and ending at the mitral valve fornix, on the contrary of the closed curves generated by equation (7a). However, the convention about the curvature sign is held easily.) Inspection of Figure 5 suggests some typical curvature patterns. The intraclass variation of curvature is relatively small. The correlation with the surgery outcome might confirm the 492
50 Apex
I
Figure 4. Interobserver variability: curvature estimated with FOURxy method from 39 test contours digitized by two observers
Application
40
Figure 5. Curvature at 90 equidistant points along the end-systolic LV contours, averaged on control (solid squares) and aneurysm (solid circles) groups
-20
-40
20 30 Anterwall
prognostic value of shape analysis, besides its diagnostic performance. Other pathological patterns in aortic insufficiency, stable angina and myocardial infarct are under investigation.
DISCUSSION Two philosophical approaches are possible for discrete curvature estimation as well as, in general, for derivative assessment from noisy data: digital filtering and model fitting. In case of very noisy data, matched filtering provides the best performance on condition that a priori knowledge of signal shape is available. By assessing the frequency spectrum and exploiting some assumptions on noise and signal characteristics, a filter bandwidth can be selected. Instead of a unique cutoff frequency, the scale-space approach’* uses a range of scales: at large values of Gaussian widths only significant structures are retained, but complex algorithms are needed to track the structures to small Gaussian widths to place them accurately. Another difficulty arises from blurring of neighbouring large-scale structures. Without neglecting the importance of the multiscale approach in vision problems, our goal is to recover a continuous representation from a discrete noisy one. In fact, LV papillary muscles and wall trabecular structures are excluded from contours since they are not visible in angiography, nor relevant. It follows that the contours under examination are highly smooth, and cannot easily be described in terms of simple curvature primitives2’. Amon model fitting techniques, autoregressive models’. 9* and Tikhonov’s regularization theory were also used. This latter method can be formulated in our terms as follows: find a function z(t,, cklk = 1, NC) that minimizes (l/N)ZnIz(t,, cJ-z,~~+A&~” z(t,, ck)dtm. The first term in this functional form is a mean squared measure of the reconstruction error; the second term is a measure of the noise oscillations; and A is a real non-negative smoothing parameter. From the methods of calculus of variations, the solution is a set of spline functions of order 2m. The value of h is critical, and may be determined through the maximum likelihood estimatior$ or by using the generalized image and vision computing
cross-validation’. However, this optimal performance is paid for by incurring long computation times and
curvature has been used as an objective system to assess wall motion2h.
reference
highly complex algorithms. Our results show that a close approximation to the true curvature values can be achieved through a rather faster algorithm, fitting trigonometric polynomials and exploiting simple heuristics to minimize both reconstruction errors and noise oscillations. Such a method was developed for test contours generated by other trigonometric polynomials. However, the frequency characteristics of these test contours are highly modified by the digitization procedure (N,z 12) with respect of the theoretical spectra of their generating shape models (N,.= 3). Thereby, the method is validated for a more general class of planar shapes, taking into account that the parameter values of the FOURxy algorithm, are ideally suited to smooth curves. In the presence of sharp corners between straight segments, for example, it will be necessary to increase the high value of NC. Whenever the resulting errors cannot be accepted, such discontinuities must be detected and processed independently. Another advantage is that the method provides Fourier descriptors which are themselves useful for shape coding and classification. Finally, even if discrete Fourier series is suited only to periodic functions such as closed contours, a simple procedure is proposed to compute curvature also for partially open or occluded curves. As to the comparison among methods, only the Gaussian filtering is able to approach the performance of the FOURxy method. The other methods based on the orientation-length function, though widely used in cardiology2”*24, systematically under-estimate the highest curvature values and overestimate the curvature valleys because of suboptimal smoothing. No regression analysis can correct this nonlinear relation versus analytical curvatures. It would be necessary to improve the angle estimation, perhaps using an adaptive window14. The FOURwy method, varying the harmonic number along the contour, is indeed an adaptive filtering. Once the curvature-length function is estimated, a number of interpretation schemes are possible. Statistical analysis has been widely used for global curve matching in the parameter space by fitting a model and using a training set. On the other hand, subjective evaluation of planar curves seems to occur in space domain. It is well known that information on the shape of a curve is concentrated at dominant points such as local maxima of curvature (Attneave), local minima (the codons representation), inflection points (zerocrossing of curvature), as well as on their combinalion (curvature primitives). Accordingly, a number of curve partitioning (or growing) algorithms were developed’.‘. The choice of the interpretation scheme depends on the application. For this reason, a further improvement in shape understanding may be expected by exploiting knowledge-based systems. In our application, statistical and artificial intelligence techniques are combined to produce symbolic descriptions in terms of primitives and relations. After the learning phase, fuzzy logic and forward-backward inference are employed in a diagnostic expert system 25. More particularly, since LV shape changes smoothly during the cardiac cycle, \lol IO no 7septemher
1992
ACKNOWLEDGEMENTS This work was supported Education (MPI).
by the
Italian
Ministry
of
REFERENCES
4 5
6
7
8
9
10
11
12
13
14
15
16
17
18
Special Issue on Shape Analysis. IEEE Trans. PAMI, Vol 8 No 1 (January 1986) Marshall, S ‘Review of shape coding techniques’, Image & Vision Comput., Vol7 (1989) pp 281-294 Hutchins, G M, Bulkley, B H, Moore, G W, Piasio, M A and Lohr, F T ‘Shape of the human cardiac ventricles’, Am. .I. Cardiology, Vo141 (1978) p 646 Pavlidis, T Algorithms for Graphics and Image Processing, Springer-Verlag, Berlin (1981) Woltring, H J ‘A FORTRAN package for generalized cross-validatory spline smoothing and differentiation’, Adv. Eng. Softw., Vol8 (1986) pp 104113 Anderssen, R S and Bloomfield, P ‘Numerical differentiation procedures for non-exact data’, Numer. Math., Vol 22 (1974) pp 157-182 Craven, P and Wahba, G ‘Smoothing noisy data with spline functions - Estimating the correct degree of smoothing by the method of generalized cross-validation’, Num. Math.. No 31 (1979) pp 377-403 Dubois, S R and Glanz, F H ‘An autoregressive model approach to two-dimensional shape classification’, IEEE Trans. PAMI, Vol 8 (1986) pp 5565 Parent, P and Zucker, S W ‘Trace inference, curvature consistency, and curve detection’, IEEE Trans. PA MI, Vol 11 (1989) pp 823-839 Lanshammar, H ‘On precision limits for derivatives numerically calculated from noisy data’, Med. Biol. Eng. & Comput., Vol 15 (1982) pp 459-470 Cappellini, V, Costantinides, A G and Emiliani, P Digital Filters and their applications, Academic Press, New York (1978) Dorsf, L and Smeulders, A ‘Length estimators for digitized contours’, Comput. Vision, Graph. & Image Process., Vol 40 (1987) pp 31 l-333 Vernon, D ‘Two-dimensional object recognition using partial contours’, Image & Vision Comput., Vol 7 No 1 (1987) pp 21-27 Teh, C-H and Chin, R T ‘On the detection of dominant points on digital curves’, IEEE Trans. PAMI, Vol 11 (1989) pp 859-872 Freeman, H and Saghri, J A ‘Analysis of the precision of generalized chain code for the representation of planar curves’, IEEE Trans. PAMI., Vol 3 (1981) pp 533-539 Zhan, C T and Roskies, R Z ‘Fourier descriptors for plane closed curves’, IEEE Trans. Comput., Vol 21 (1972) pp 269-281 Kaiser, J F ‘Digital filters’, in F F Kuo and J F Kaiser (eds), System Analysis by Digital Computer, Wiley, New York (1966) Witkins, A ‘Scale-space filtering’, Proc. 7th Int. 493
Joint
Conf
Artif.
Intell.,
Kalsruhe,
Germany
(1983) p 1019 19 Kuhl, F P and Giardina, C R ‘Elliptic Fourier features of a closed contour’, Comput. Graph. & Image Process., Vol 18 (1982) pp 236-258 20 Persoon, E and King-Sun, F ‘Shape discrimination using Fourier Descriptors’, IEEE Trans. PAMZ, Vol 8 (1986) pp 388-397 21 Asada, H and Brady, M ‘Curvature primal sketch’, IEEE Trans. PAMI, Vol 8 (1986) pp 2-14 22 Chellappa, R and Bagdazian, R ‘Fourier coding of image boundaries’, IEEE Trans. PAMZ, Vol 6 (1984) pp 102-105 23 Mancini, G J, LeFree, M and Vogel, R A ‘Curvature analysis of normal ventriculograms: fundamental framework for the assessment of shape changes in man’, Proc. Comp. in Card. (1985) pp 141-144 24 Duncan, J S, Smeulders, A, Lee, F and Zaret, B L ‘Measurement of end diastolic deformity using bending energy’, Proc. Comp. in Card., Washington, DC, (1988) p 277-280 2.5 Baroni, M, Barletta, G and Fantini, F ‘LV wall motion: from quantitative analysis to knowledgebased understanding’, Proc. Comp. in Card., Jerusalem, Israel (1989) pp 483-86 26 Baroni, M, Barletta, G, Fantini, A, Toso, A and Fantini, F ‘Assessing LV wall motion by frame-toframe curvature matching and optical flow estimation’, Proc. Comput. in Card., Venice, Italy (1991)
APPENDIX
A: LIST OF SYMBOLS
R = radius of curvature. K(t), Ki = continuous or discrete curvature. i, k, n = integer indexes. t = parameter.
I = arc length. L = perimetral length. @ = slope angle, orientation-length function. 0 = cumulative orientation-length function. w = number of points in half a window. u = Gaussian function width. uz = (total) noise variance. (T: = quantization noise variance. m = derivative order. M = number of bits. B = signal bandwidth. F, = sampling frequency. f. = fundamental frequency. fc = cutoff frequency. fc,, = normalized cutoff frequency. N, = cutoff harmonic number. hk = filter samples. Kh = number of filter samples. G,(t) = Gaussian function. (Y= Kaiser window parameter. A = regularization parameter. 1. .I = modulus of complex function. N, = noise power. ck = complex Fourier coefficients. z, = complex contour data. c, KT, NT, eT= COWSantS. A,i, Ayi, Fxi, Fyi = amplitudes and phases of the ith sinusoidal component for rectangular coordinates of ventricle-shaped curves. e,, ,u<, Us = reconstruction errors and their mean and s.d. values.
x, y = rectangular coordinates of contours. N = number of contour points.
494
image and vision computing
Curvature estimation from noisy digital curves is addressed for shape analysis of angiographic left ventricle (LV) images. A simple and fast algorithm, based on Fourier series approximation, changes the harmonic number and filter window along the closed (or partially open) contours to optimize smoothness and reconstruction errors. The curvatures estimated on a set of ventricle-shaped models are significantly closer to their analytical values than those assessed by four other methods. Preliminary results on LV data are presented. Keywords: left ventricular shape, curvature estimation, adaptive filtering
INTRODUCTION An important goal in image processing is shape description and recognition. As a result of the imaging process, a 3D object is mapped onto a bidimensional projection. This image is digitized with a square sampling grid and segmented into a bilevel picture or a boundary is extracted from it. Both region and boundary representations define the ‘shape’ of the object on a background, and in most situations still carry the main information needed for recognition purposes. Even with careful acquisition and segmentation procedures, various sources of errors are possible; moreover, both projection and sampling operators are non-invertible: therefore, digital shape analysis is not a simple problem, even if we restrict ourselves to planar shapes. Many papers have been written on this subject, and excellent reviews can be found’,2. Shape analysis
Barletta”
methods can be classified in series (global) or piecewise (local) representations, internal (area-based) or external (boundary-based). These latter can further be classified with regard to the boundary parameterization (rectangular, radial, angular or curvature parametrization). For example, the method of moments is global and internal; Fourier transform is a global method which can be internal (if applied in 2D) or external, using any parametrization of the boundary. Examples of local methods are Medial Axis Transformation (internal) and polygonal approximations of the boundary. Our attention is concentrated on boundary representations of simply connected curves, because we are interested in the study of left ventricle (LV) shape as it appears in standard cardiac imaging (contrast angiography, ecography). Endocardial shape, in particular, plays an important role in ventricular function’, and its analysis provides a promising tool for clinical and pathophysiological investigations. For this purpose, the curvatrure representation is chosen because it satisfies all the following criteria:
Dipartimento di Tngegneria Elettronica, via S Marta 3, 50130 Firenze, Italy ‘Institute of Internal Medicine, Cardiology Unit, University of Florence, Florence, Italy
1. Curvature is invariant under rotation and translation of the curve, and can easily be normalized with respect to scale changes. as it is possible to 2. It is a unique representation, match the curvature along a curve against that of homogeneous populations of curves. after suitable normalizations. able to carrying informa3. It is a local representation tion at varying resolutions, and to provide a structural description in terms of primitives and relations. 4. It can represent both closed and open curves, even in the case of missing data. to small changes in part of 5. It is relatively insensitive the curve and, theoretically, no arbitrary choice is needed for its computation.
Paper received: 5 April 1991; revked papers received: 1991
However,
0262-8856/92/007485vol IO no 7september
1992
2 September
10
0
the curvature
1992 Butterworth-Heinemann
computation
involves
first and
Ltd 485
second order differentiation of noisy digital data, so smoothing is necessary to obtain significant values. Polynomial fitting and standard spline interpolation4 are highly sensitive to noise as the equivalent filtering bandwidth is suboptimal. Regularization theory with spline functions” or Fourier series’ is characterized by high performance at the cost of long computation times. As an alternative, lowpass filtering may be optimized by estimating the signal bandwidth’. In this paper a new fast algorithm for digital curvature computation is proposed and compared with another four methods. Three of them are based on the orientation-length function which is differentiated by using either differences taken over a window, or FIR filtering, with a fixed number of points or a fixed sampling step. The fourth method consists of Gaussian filtering of contour coordinates. Our algorithm is a windowed Fourier series approximation of contours, in which the number of harmonics and filter-window are locally chosen to minimize the reconstruction errors, and to maximize the smoothness of the curve. The comparison was carried out on analytically derived curvatures of a set of ventricle-shaped models, using both quantized software-generated curves (here called shape models) and manually digitized contours of the same shapes (test contours). Moreover, preliminary results on LV data are also presented.
METHODS Curvature
and noise
The curvature of a planar curve at a point P is defined as the reciprocal of the radius of the osculating circle (the circle tangent to the curve at P with the highest degree of contact): K=
l/R
(I)
From differential geometry, the curvature can be expressed as the derivative of the unitary tangent vector with respect to arc length 1, measured along the curve from a starting point. Therefore, if we let Q(f) be the slope angle between the local tangent line and some arbitrary x-axis, the (scalar) curvature, k(f), is given by: K(1) = d@ldl If the curve is expressed can write: d@ -=-.dl
d@
dx
dx
dl
It is known
(2) in explicit
form, y = y(x),
that
dl= (dx2 + dy2)“2 = dx[(l + (dyldx)2]“2, dyldx, then: d@ -= dx
one
d2yldx2 1+ (dyldx)2
and
tg@ =
It follows
that: d2yldx2
K(x) = [l + (dyldx)2]“2 Since y(x) is a multivalued function, a parametric form x =x(t) and y = y(t) is more suitable if the parameter t varies monotonically along the curve. The simplest choice is t = 1, although t might be equal to the angle defined by radius vectors projected between the curve and its centroid. This choice, however, leads again to multivalued functions for non-convex shapes, unless one resorts to some artifice’. Now, letting ’ and ’ be the first and second derivatives with respect to t (i.e. x’ = dxldt, x” = d;l d;, etc.), one can write: dy dy ~=~‘~=x’
dt
y’
and d2y -= dx2
Finally,
equation
d(y’/x’) p= dx
x’y”-y’x” X
I2
.-
1 X’
(3) becomes
x’y”_y’x”
K(t) =
(4)
(x12+yt2)3/2
The curvature sign depends on vector direction: for closed curves it is conventionally taken as positive for the convex segments and negative for the concave ones. In discrete terms, when a collection of points is available, rather than a continuous function, equation (1) is implemented either using the cocircularity concept’ or with a simple geometric construction based on triplets of pixels4. The angles in equation (2) are locally estimated, and the derivatives in equations (2,3, 4) may be replaced by first differences. However the discrete curvature computation faces the problem of noise, given by both quantization and statistical errors. Here, the signals underlying discrete noisy data are assumed to be differentiable, and are therefore certainly smooth and hence practically band-limited as their spatial frequency spectra will be concentrated at low frequencies. In general, given a signal with limited bandwidth (B), sampled at equispaced steps (llfs), affected by additive stationary white noise (~3, the uncertainty in the computation of the mth order derivative’” is given by: &CT;.
B(2m+‘)l[fs* (2m + l)]
(5)
Therefore, the sampling frequency must be chosen not only with respect to signal bandwidth, but also in relation to the noise level and derivative order. In the one-dimensional case, the samples are also quantized with M bits. This is equivalent to an additive white noise” with variance given by: a; = 2-2M/12
(6)
In two dimensions, the sampling requirements to preserve shape are more restrictive, and if d is the linear dimension of the finest details, the Nyquist image and vision computing
condition
Shape models and test contours
(llfs) d (d/2) becomes4:
(llfs) d (d/W2)
Since the curve x(t), y(t) is sampled have:
on a square grid we
yi = int b(i/fJ
x, = int [x(i/fy) + 0.51,
+ OS]
and successive slope angles can differ only by a multiple of 45”. These abrupt changes in direction are amplified by differentiation and give rise to a still coarser quantization in curvature. Even the convexity property of the curve may be locally violated. This means that a certain order of smoothness must be imposed to digital curves. However the scale at which the smoothing is performed strongly influences the result: with over-smoothing small scale details are lost, while under-smoothing leaves undesirable noise fluctuations. Another problem is related to the non-uniform sampling of the length: in fact, the distance between two consecutive pixels is 1 along the sampling grid directions and V? in diagonal directions. For this reason, the digital contour must be resampled by linear interpolation with spacing equal to PIN, where N is the number of equispaced pixels, and the perimeter P is estimated along the original curve’2,‘3.
Shape models
The LV silhouttes seen in standard cardiac imaging are smooth curves with high and low, positive and negative More particularly, the LV contours curvatures. obtained in contrast angiography with RAO projection are characterized by a highly variable valvular zone and two low curvature regions separated by a high curvature apical zone. For the purpose of obtaining a set of ventricle-shaped closed curves, the parametric equations of the ellipse were generalized to allow accumulation of high order terms, just like a truncated Fourier series: x(t) = x0 + Ji,,A ,,I sin (nt + F,,,). y(t) = yo+ &,Avnsin
(nt+
F,.,,), n=l,2,3,1~[0,2rr],
(7a)
where (X,,, Y,) is the central point, the number of terms was limited to 3, and the other constants were chosen by means of an interactive graphical procedure. Table 1 lists the constant values chosen for 39 models, and some of them are showed in Figure 1. Curvature can be analytically computed according to equation (4) by differentiating equations (7a), i.e. for x(t): x’ (t) = 2;,, nA i,, cos (nt + F,,) , x”(t) = -Ii,,
n2 A,,, sin (nt + F,,,),
(7b)
Table 1. Model coeffkients for
Constant values A,,
A,.,
Fx,
Fw
AC;_7AYZ Friz FYJ
+60
+50
+40
+90
+20
0
+70
+ 80
+40
+90
+20
0
+70
+20
-t 90
+20
+80
+20
Ax?
A,,
F,.I Fy.7
0
0
+20
0
+50
+70
0
0
+20
0
+50
0
+70
0
0
+20
0
+50
0
+70
0
0
+20
0
+75
+70
0
0
+20
0
+50
39 ventricle-shaped curves
i + 80
4 5
I
f60
I
IO
+I10
II
+80
I
I
+40
I5 16
+x0
+70
+40
I
I +100
20 21
I 26 27 I 32 33
I
39
+ 80
+70
+40
-t 100
+20
+5
+30 +70
+70
+40
+90
+20
0
+70
0
0
+10
0
+50
+50
+40
+110
+20
0
+70
0
0
+5
0
+70
I
+95 +60
t i +80
WI IO no 7september
1992
487
The orientation-length function Q(I) can be approximated by taking either cosine or tangent measures on successive points, as reviewed by Teh and Chint4. Here equation (2) was rewritten as:
(8) where {Xi, YJi =:1 . . . N} are the coordinates of the digitized contour, and w B 1 is an integer parameter which defines the region of support (window) for the computation of differences. Choosing a w> 1 window has the effect of improving the direction quantization in a way that is similar to high order chain codes”.The angular quantization step can be written as: Figure I. Four sets of ventricle-shaded curves (3rd, Sth, 6th and 7th set in Table It Therefore, these shape models provide a gold standard to compare the methods for curvature estimation described in the next sections. They were used in two ways: 0
0
the coordinates generated by equations (7a) were sampled at steps of 2”, rounded to integer values, and directly given as input to the analysis programs; computer displays of shape models were printed on paper and manually digitized by two experienced observers using a 0.1 mm resolution graphic tablet, according to the usual procedure for digitizing actual LV data. The discrete coordinates of these test contours were resampied at 1 mm equispaced steps, and analysed by the same programs.
The former shape models generate a set of curves affected only by quantization noise (SNR, ~25 dB), while the latter test contours represent the worst case in which the statistical errors due to manual tracing are also present. Moreover, the points sampled using equation (7a) are equispaced in angle not in arc length, so that they were resampied at equidistant steps by an iterative procedure to allow a point-by-point comparison between the analytical curvatures and those estimated from test contours. Finally, since the number of points varies with contour size, the comparison was carried out by linearly interpolating ail the curvature values at N equidistant points. This size normalization implicitly assumes that equal length abnormalities are more significant in smaller ventricles. A value of N = 90 was chosen because the analytical curvature sampled at a lower number of points underscores the sharper peaks seen along our shape models.
Window-based method Now the methods
introduced. 488
used for curvature
estimation
are
@i+w-
and
Q,i_,,=2?T/ZM
M=w+2
(9)
Of course, the window size must be large compared with the spatial noise frequency so that it was chosen by the following procedure. Since the curvature function along a closed curve is periodic, Fourier analysis can be properly applied. The amplitude spectra of the curvatures, ana~yticaiiy computed on the shape models with those previously described, were compared obtained on the test contours using equation (8) and varying w in the range [2, 201. The best signal-to-noise ratio was achieved with a window size of w = 12, using 1 mm resolution.
Filtering the orientation function In the previous window method the curvature at each point depends on the orientation value in only two distant points, neglecting ail the others. This does not ensure the elimination of spurious oscillations, so a further smoothing of curvature values is needed. A better use of the available data is obtained by averaging all the orientation values over the window, or filtering them to reduce high frequency noise. For this purpose it is necessary to define a cumulative orientation-length function which is continuous along the contour. In fact, since the arctangent function returns angles between -7r and +n-, the slope angle function Q(j) has artificial discontinuities, i.e. for a closed curve Q(O) -Q(L) = 27r. At these points an offset of fn= was added, and the arc length IE [O, L], was replaced by the parameter te[O, 27r]ih. The resulting angular function: B(t) = @(tL/27)
- Q(O) + t
(10)
is periodic, and can be expanded in discrete Fourier series. Its power spectrum j O(f) 1’ was used to select the filter bandwidth. If we assume a white noise, and that f,22B/c (c = 0.8), the noise power No can be assessed by averaging the spectrum between c. fJ2 and f,/2. Now the cutoff frequency fc where the signai-tonoise-ratio reaches a threshold (KT> 1) can be chosen: i@(fc)?=
Kr+No
(11)
A FIR filter was then designed to reduce the frequency image and vision computing
Having chosen N = 90 for our ventricle-shaped tours and estimated NC = 12 by means of equation we have f,., = 12/90 = 0.1333.
con(1 l),
Filtering contour data
-50
1
0
0.1
I
0.2
0.3
0.4
0.5
Norml~r~dfwqumcy
Figure 2. Frequency response of digital filter used to differentiate slope angle data: f,, = 0.1333, CY= 3.39 with f >fC and, at the same time, to all the others, according to equation (2). response of a differentiator taken at an of samples” is given by
components differentiate The impulse odd number ho = 0.
hk = [2z-kf,, cos(2nkf,.,) - sin (2rkfC,)]/7rk2
where the cutoff frequency was normalized sampling frequency to provide a more general
(12) by the design:
(13) The filter samples must be truncated with a well shaped window to obtain a finite number of coefficients kE[--/,, &I. We used the Kaiser window” with Kh = 5 and shape factor (Y= 3.39, which offers a good compromise between sharp transition band and residual passband oscillations (see Figure 2). This filter was applied to the angle data estimated with equations (10) and (8) using w = 1 and two different sampling schemes. In the first, the test contours were sampled with a fixed frequency so a different number of points is obtained for different curves. In the second, the number N of samples is taken as fixed, and the spatial frequency is variable. In the latter case, the curvature values do not depend on size, while in the former, as well as in all the other methods discussed here, it is necessary to normalize the shape representation for size variations by multiplying the curvatures by the perimetral length. For comparison purpose, the curvatures computed with the fixed N were normalized by an average perimeter. Another difference is in the value of fCn. This is given by equation (13) once fixedf, (1 mm-‘) and estimated fC by equation (11) (fCn ~0.05). On the other hand, it is convenient to choose the fixed value of N as small as possible to reduce the quantization step (9) until it reaches the minimum value given by the sampling theorem. By writing fC = N,. f. as a multiple of the fundamental frequency f0 = 1lL = f,lN, equation (13) becomes: L.,, = NJN ~~1110 no 7 september 1992
(14)
Another common method for differentiating noisy data involves Gaussian filtering, according to Marr’s theory. Using the rectangular parametrization of the contours the x(t) and y(t) functions can be smoothed independently by convolving them with a Gaussian kernel of width CJ, G,(t) = exp(-t’/2a’)/a(2~)‘“. To compute curvature the smoothed x5(t) function must be differentiated twice, and letting * be the convolution integral, x,y’(t) = [x(t) :kG,,(t)]’ = x(t) * G,‘(t). x,“(t) =x(t) * G,“(r), where G,,‘(t) = [-tV~~(27r)“1)] G,,“(t) = [W’(27r)“‘]
exp( -tt’/2rr?) exp(-t’/2u’)(
1 - f’/a*)
The same holds for ys(t). The Gaussian function realizes the best compromise between frequency bandwidth and spatial window, according to the indetermination principle. Moreover, it can be easily truncated at a finite number of coefficients where it falls below a very low threshold. Finally, it has the unique property of not introducing extra zero-crossings of the second derivative’. Therefore, since the zerocrossings indicate the dominant points along the curve (the inflection points, in our case), the Gaussian function was used for multiscale filteringlx. In this work, however, the width value v was chosen proportionally to the number N of points, depending upon noise level and shape details.
Improving
Fourier series approximation
Starting again from a rectangular parametrization of closed contours, the x(t) and y(t) functions are periodic and can be independently expanded into a Fourier series. A truncated expansion can be used to reconstruct a smoothed and/or interpolated version of the contour and, at the same time, to compute curvature according to equations (4) and (7b). The trigonometric expansion of x(t) and y(t) can be expressed in a matrix form, and geometrically interpreted as a sum of ellipses rotating with a speed proportional to their harmonic number”. Similarly, the coordinate functions can be considered as real and imaginary parts of a complex discrete function z(n): z(n) =x(n)
+ jy(n) = Ck=-Nc,Nc hkck e’k(?T’N)n
(15)
where the complex Fourier coefficients are computed via FFT with N = 256 in the usual form: ck
=
(l/N)
&+,,Z,e-~“(2”‘N)k
The weighting factors hk are introduced reduce the spurious oscillations given by any of a trigonometric series. For this purpose, window was used with a shape factor
(16) in (15) to truncation the Kaiser cy chosen 489
proportionally to NC, starting from 3.39. (This filter is the dual form of that used for the orientation function.) As pointed out by Persoon and King-SunzO, any truncated expansion (15) converges to a closed curve, but the reconstructed points are not always equispaced in arc length, so they must be resampled during the same program loop implementing equations (15) and (4). Our concern is now to select the harmonic number NC. Using the amplitude spectra (16) averaged over all test contours does not ensure the best performance in each case. In fact, rewriting equation (ll), it follows that Ic(?Nc)/ 2~ K,N, is satisfied for KT = 20 at N,.=+12, for Kr=6 at NC=+24 and for Kr=2 at NC= +36. (T a k’mg into account the fact that the spectrum of the complex function z(n) is not symmetric, unlike the real case where ck * = c-k, we have considered an average of both the positive and negative harmonic number.) If a low value is chosen for NC (in our context N,.d24), the reconstruction errors e, defined as the Euclidean distances between the original and reconstructed points, are much greater for some test contours than for the smoother ones. In other words, high curvatures are underestimated (see Figure 3a). On the other hand, by increasing NC the errors are reduced on the average, but at the cost of leaving spurious oscillations (see Figure 3b). Further filtering, such as increasing the Kaiser window parameter LY, has an effect similar to that of decreasing NC. A possible approach consists of using a high number of harmonics only for those contours in which a low NC value gives a maximum error above some threshold. Nevertheless, the choice of this threshold may be critical. Another approach stems from the following empirical observation. In the reconstructed test contours the most unrealistic oscillations are apparent in smoother regions, where higher harmonics are redundant (see Figure 3b). Therefore, it seems more reasonable to vary NC locally along each contour (see Figure 3c), instead of choosing NC from curve to curve. A problem now arises to avoid any discontinuity where the value of NC is changed (subsequently referred as transition points). As a first solution, the following simple algorithm was tested (FOURry method): 1. compute FFI using N = 256 points resampled at equidistant steps along the contour perimeter (Equation 16); 256 points using equation (15) with NC = 2. reconstruct 12, calculating the euclidean errors e,, their peak value emax, mean value pu,, and standard deviation cr,? are labelled as 3. the points with e,,
a
-4od b
C Figure 3. (a) Normalized curvatures estimated using N, = 12 and (Y = 3.39 against those analytically computed for test contour no. 2; (b) as 3a using N, = 36 and (Y = 9; (c) results by FOURxy method using same data as 3a, where e,
Both steps (4) and (5) accomplish a sort of morphological filter: removing isolated good points is equivalent to an erosion operation, and precedes the dilation operation to prevent macroscopic discontinuities where the image and vision computing
smoothed contour comes across the original one. However, such discontinuties may be apparent in any point of the transition regions. So it is better to restrict these transition points only between runs of good and bad points. To this aim, the sixth step in the earlier algorithm was removed, and the curvature values at each transition point, as well as in its two neighbouring points, were linearly interpolated. In conclusion, it is worth noting that this method, though based on discrete Fourier series, is not only suited for closed curves. In fact, the other methods based on filtering require some extrapolation of open curves, until half the filter window is beyond the endpoints. Whenever a similar extrapolation could be accomplished for closing the contour, Fourier series can be employed to compute curvature also for open curves. In the angiographic LV contours the aortic plane is usually closed by a straight segment. This may introduce one (or two) corners not existing in the endocardial wall. By extrapolating the contour beyond the valvular margins and then closing it with a straight segment, curvature can be properly measured along the entire LV silhouette. Among the available techniques we used the followextrapolation: x(n + i) = 2x(n) ing symmetrical x(n-i), where i=l, 5 for n=Nand i=-I, -5 for II = 1. The same holds for y.
RESULTS Validation
of methods
Several tests were undertaken to validate the five methods previously described, using both shape models and test contours. Test programs were written in Pascal and run on a 80386 computer (25 MHz, 80387 coprocessor) with execution times ranging from 2s for the window method to 15s for the FOURxy method. First, curvature invariance with rotation, translation and magnification was verified. In particular, the curvature values computed before and after 60” and 90 rotation were found to be identical up to the fifth decimal digit. It is well known that curvature depends on size, the curvature of a circle with radius R being twice the curvature of a circle with radius 2R. So the curvature values, expressed in mm-‘, were multiplied by the perimetral length. This normalized curvature is dimensionless, and for any circle is conventionally equal to 27~. After scaling the test contours by factors of 0.5 and 2, the changes in normalized curvature were less than 1%. Then the accuracy of the curvature estimated by the five methods from discrete noisy rectangular coordinates were investigated with respect to their analytical values. The comparison was accomplished by linear regression and correlation analysis. Table 2 shows the results of this comparison for the set of 39 ventricle-shaped models listed in Table 1. Similar results were achieved for circular and elliptic curves. Obviously, our algorithm (FOURxy) gives curvature values identical to the analytical ones when using
Table 2. Curvature estimated on shape models by five methods against analytical values. (WIN12 = Window method with w = 12); FIRP-FP or FIRP-FS = Filtering of the orientation function using N = 90 (Fixed Points) or f, = 1 mm-’ (Fixed Sampling); FOURXY = Fourier approximation with N, = 12, LY= 3.39 where e, <2pu,; N, = 36, (Y= 9 elsewhere; GAUSS = Gaussian filtering with u = 6 Method
Correlation
Intercept
Slope
SEE
WIN12 FIRP-FS FIRP-FP FOURXY GAUSS
0.5!8 0.627 0.629 0.998 O.Y46
5.68 1.Y6 2.02 0.28 I.51
0.36 0.20 0.21 O.Y4 0.68
x.22 3.41 3.61 0.82 3.1’)
NC = 3. However, even if the harmonic number is not matched to that of the generating model (equation (7a)), the results are better (R>0.99, standard error of estimate - S.E.E. - ~1) than those of all the other methods. In addition, noise oscillations are also minimal: in particular, the curvatures estimated by Gaussian filtering with CT= 2 are very close to the true values (R = 0.992, slope = 0.98) but fluctuate about them. This can be seen by comparing the derivatives of the estimated and true curvature-length functions. Table 3 shows the results of the comparison for the set of 39 test contours. After manual digitization the errors increase. Also the spatial frequency components of the rectangular coordinate signals increase, so that the reconstruction errors are acceptable (order of magnitude 1 pixel) only for N,.z 12. An assessment of the observer variability is given in Figure 4. Following the t-test for paired data, the mean interobserver difference was found to be about 1%) whereas the mean deviation of the FOURxy method from the analytical data was about 10%. According to correlation analysis, FOURxy showed R > 0.90. On the other hand, Gaussian filtering with g E [2, 61 gives R = 0.70, with a low reproducibility for ~<4 due to residual oscillations. The window-based method was the worst (R <0.5 with w = 12): moreover, the resulting curvature-length function exhibited high frequency noise. The method based on filtering of the orientation Table 3. Curvature estimated on test contours by five methods against analytical values. (WIN12 = Window method with w = 12; FIRP-FP or FIRP-FS = Filtering of the orientation function using N = 90 (Fixed Points) or f, = 1 mm-’ (Fixed Sampling); FOURXY = Fourier approximation with N, = 12, cy = 3.39 where e, < 2,~~; N, = 36, LY= 9 elsewhere; GAUSS = Gaussian filtering with u = 6)
WIN12 FIRP-FS FIRP-FP FOURXY GAUSS
OBSI 0BS2 OBSI OBS2 OBSI OBS2 OBSI OBS2 OBSl OBS2
Correlation
Intercept
Slope
SEE
0.290 0.470 0.314 0.315 0.820 0.821 0.919 0.907 0.704 0.723
5.45 4.67 4.80 4.81 2.20 2.1’) I.25 1.32 2.53 2.57
0.22 0.37 0.08 0.08 0.82 0.81 0.88 0.85 0.63 0.62
22.x 22.0 7.33 7.31 6.32 6.29 4. I8 4.39 7.04 6.60
80
60
I
40 N t
Em
d
0
-L 1 ,
10
-20
0
20 Observer
I
I
40
60
60
70 80 mfer. wall
90
80
function gave different results for fixed N (R = 0.31) with respect to taking a fixed sampling step (R = 0.82), but the latter method showed a worse performance in smoothness and variability. In conclusion, the FOURxy method outperforms the others, in terms of both reconstruction errors and noise oscillations, in estimating curvature of ventricle-shaped curves. Moreover, it is relatively insensitive to significant changes in the parameter values. The error threshold eT was changed from pu,+cP to 2,uu,, and also a fixed value was used (er = 2), without significant differences in the estimated curvature. The Kaiser window parameter CYhas no critical value, as well as the morphological filter parameter NT, if limited to few units. Some minor differences were found by changing the low harmonic number (i.e. from 12 to 18), since some test contours have very high curvature variations. The same differences were found by using a lower value for er. On the contrary, no significant difference was found varying the high harmonic number in the range [30, 481. to contrast ventriculography
Our algorithm was applied to the X-ray LV cineangiograms in RAO projection. The LV contours of 16 patients who were found to have normal cardiac anatomy and function were digitized frame-to-frame from 35 mm films, using a moviola and a graphic tablet. Similarly, the LV contours of about 40 patients with aneurysm following anterior infarct were digitized at the end of systolic and diastolic phases. The endsystolic curvature at 90 equidistant points along the contours, averaged on normal and pathological groups, are shown in Figure 5. (The LV contours are traced in the clockwise direction, starting from the anterior aortic margin and ending at the mitral valve fornix, on the contrary of the closed curves generated by equation (7a). However, the convention about the curvature sign is held easily.) Inspection of Figure 5 suggests some typical curvature patterns. The intraclass variation of curvature is relatively small. The correlation with the surgery outcome might confirm the 492
50 Apex
I
Figure 4. Interobserver variability: curvature estimated with FOURxy method from 39 test contours digitized by two observers
Application
40
Figure 5. Curvature at 90 equidistant points along the end-systolic LV contours, averaged on control (solid squares) and aneurysm (solid circles) groups
-20
-40
20 30 Anterwall
prognostic value of shape analysis, besides its diagnostic performance. Other pathological patterns in aortic insufficiency, stable angina and myocardial infarct are under investigation.
DISCUSSION Two philosophical approaches are possible for discrete curvature estimation as well as, in general, for derivative assessment from noisy data: digital filtering and model fitting. In case of very noisy data, matched filtering provides the best performance on condition that a priori knowledge of signal shape is available. By assessing the frequency spectrum and exploiting some assumptions on noise and signal characteristics, a filter bandwidth can be selected. Instead of a unique cutoff frequency, the scale-space approach’* uses a range of scales: at large values of Gaussian widths only significant structures are retained, but complex algorithms are needed to track the structures to small Gaussian widths to place them accurately. Another difficulty arises from blurring of neighbouring large-scale structures. Without neglecting the importance of the multiscale approach in vision problems, our goal is to recover a continuous representation from a discrete noisy one. In fact, LV papillary muscles and wall trabecular structures are excluded from contours since they are not visible in angiography, nor relevant. It follows that the contours under examination are highly smooth, and cannot easily be described in terms of simple curvature primitives2’. Amon model fitting techniques, autoregressive models’. 9* and Tikhonov’s regularization theory were also used. This latter method can be formulated in our terms as follows: find a function z(t,, cklk = 1, NC) that minimizes (l/N)ZnIz(t,, cJ-z,~~+A&~” z(t,, ck)dtm. The first term in this functional form is a mean squared measure of the reconstruction error; the second term is a measure of the noise oscillations; and A is a real non-negative smoothing parameter. From the methods of calculus of variations, the solution is a set of spline functions of order 2m. The value of h is critical, and may be determined through the maximum likelihood estimatior$ or by using the generalized image and vision computing
cross-validation’. However, this optimal performance is paid for by incurring long computation times and
curvature has been used as an objective system to assess wall motion2h.
reference
highly complex algorithms. Our results show that a close approximation to the true curvature values can be achieved through a rather faster algorithm, fitting trigonometric polynomials and exploiting simple heuristics to minimize both reconstruction errors and noise oscillations. Such a method was developed for test contours generated by other trigonometric polynomials. However, the frequency characteristics of these test contours are highly modified by the digitization procedure (N,z 12) with respect of the theoretical spectra of their generating shape models (N,.= 3). Thereby, the method is validated for a more general class of planar shapes, taking into account that the parameter values of the FOURxy algorithm, are ideally suited to smooth curves. In the presence of sharp corners between straight segments, for example, it will be necessary to increase the high value of NC. Whenever the resulting errors cannot be accepted, such discontinuities must be detected and processed independently. Another advantage is that the method provides Fourier descriptors which are themselves useful for shape coding and classification. Finally, even if discrete Fourier series is suited only to periodic functions such as closed contours, a simple procedure is proposed to compute curvature also for partially open or occluded curves. As to the comparison among methods, only the Gaussian filtering is able to approach the performance of the FOURxy method. The other methods based on the orientation-length function, though widely used in cardiology2”*24, systematically under-estimate the highest curvature values and overestimate the curvature valleys because of suboptimal smoothing. No regression analysis can correct this nonlinear relation versus analytical curvatures. It would be necessary to improve the angle estimation, perhaps using an adaptive window14. The FOURwy method, varying the harmonic number along the contour, is indeed an adaptive filtering. Once the curvature-length function is estimated, a number of interpretation schemes are possible. Statistical analysis has been widely used for global curve matching in the parameter space by fitting a model and using a training set. On the other hand, subjective evaluation of planar curves seems to occur in space domain. It is well known that information on the shape of a curve is concentrated at dominant points such as local maxima of curvature (Attneave), local minima (the codons representation), inflection points (zerocrossing of curvature), as well as on their combinalion (curvature primitives). Accordingly, a number of curve partitioning (or growing) algorithms were developed’.‘. The choice of the interpretation scheme depends on the application. For this reason, a further improvement in shape understanding may be expected by exploiting knowledge-based systems. In our application, statistical and artificial intelligence techniques are combined to produce symbolic descriptions in terms of primitives and relations. After the learning phase, fuzzy logic and forward-backward inference are employed in a diagnostic expert system 25. More particularly, since LV shape changes smoothly during the cardiac cycle, \lol IO no 7septemher
1992
ACKNOWLEDGEMENTS This work was supported Education (MPI).
by the
Italian
Ministry
of
REFERENCES
4 5
6
7
8
9
10
11
12
13
14
15
16
17
18
Special Issue on Shape Analysis. IEEE Trans. PAMI, Vol 8 No 1 (January 1986) Marshall, S ‘Review of shape coding techniques’, Image & Vision Comput., Vol7 (1989) pp 281-294 Hutchins, G M, Bulkley, B H, Moore, G W, Piasio, M A and Lohr, F T ‘Shape of the human cardiac ventricles’, Am. .I. Cardiology, Vo141 (1978) p 646 Pavlidis, T Algorithms for Graphics and Image Processing, Springer-Verlag, Berlin (1981) Woltring, H J ‘A FORTRAN package for generalized cross-validatory spline smoothing and differentiation’, Adv. Eng. Softw., Vol8 (1986) pp 104113 Anderssen, R S and Bloomfield, P ‘Numerical differentiation procedures for non-exact data’, Numer. Math., Vol 22 (1974) pp 157-182 Craven, P and Wahba, G ‘Smoothing noisy data with spline functions - Estimating the correct degree of smoothing by the method of generalized cross-validation’, Num. Math.. No 31 (1979) pp 377-403 Dubois, S R and Glanz, F H ‘An autoregressive model approach to two-dimensional shape classification’, IEEE Trans. PAMI, Vol 8 (1986) pp 5565 Parent, P and Zucker, S W ‘Trace inference, curvature consistency, and curve detection’, IEEE Trans. PA MI, Vol 11 (1989) pp 823-839 Lanshammar, H ‘On precision limits for derivatives numerically calculated from noisy data’, Med. Biol. Eng. & Comput., Vol 15 (1982) pp 459-470 Cappellini, V, Costantinides, A G and Emiliani, P Digital Filters and their applications, Academic Press, New York (1978) Dorsf, L and Smeulders, A ‘Length estimators for digitized contours’, Comput. Vision, Graph. & Image Process., Vol 40 (1987) pp 31 l-333 Vernon, D ‘Two-dimensional object recognition using partial contours’, Image & Vision Comput., Vol 7 No 1 (1987) pp 21-27 Teh, C-H and Chin, R T ‘On the detection of dominant points on digital curves’, IEEE Trans. PAMI, Vol 11 (1989) pp 859-872 Freeman, H and Saghri, J A ‘Analysis of the precision of generalized chain code for the representation of planar curves’, IEEE Trans. PAMI., Vol 3 (1981) pp 533-539 Zhan, C T and Roskies, R Z ‘Fourier descriptors for plane closed curves’, IEEE Trans. Comput., Vol 21 (1972) pp 269-281 Kaiser, J F ‘Digital filters’, in F F Kuo and J F Kaiser (eds), System Analysis by Digital Computer, Wiley, New York (1966) Witkins, A ‘Scale-space filtering’, Proc. 7th Int. 493
Joint
Conf
Artif.
Intell.,
Kalsruhe,
Germany
(1983) p 1019 19 Kuhl, F P and Giardina, C R ‘Elliptic Fourier features of a closed contour’, Comput. Graph. & Image Process., Vol 18 (1982) pp 236-258 20 Persoon, E and King-Sun, F ‘Shape discrimination using Fourier Descriptors’, IEEE Trans. PAMZ, Vol 8 (1986) pp 388-397 21 Asada, H and Brady, M ‘Curvature primal sketch’, IEEE Trans. PAMI, Vol 8 (1986) pp 2-14 22 Chellappa, R and Bagdazian, R ‘Fourier coding of image boundaries’, IEEE Trans. PAMZ, Vol 6 (1984) pp 102-105 23 Mancini, G J, LeFree, M and Vogel, R A ‘Curvature analysis of normal ventriculograms: fundamental framework for the assessment of shape changes in man’, Proc. Comp. in Card. (1985) pp 141-144 24 Duncan, J S, Smeulders, A, Lee, F and Zaret, B L ‘Measurement of end diastolic deformity using bending energy’, Proc. Comp. in Card., Washington, DC, (1988) p 277-280 2.5 Baroni, M, Barletta, G and Fantini, F ‘LV wall motion: from quantitative analysis to knowledgebased understanding’, Proc. Comp. in Card., Jerusalem, Israel (1989) pp 483-86 26 Baroni, M, Barletta, G, Fantini, A, Toso, A and Fantini, F ‘Assessing LV wall motion by frame-toframe curvature matching and optical flow estimation’, Proc. Comput. in Card., Venice, Italy (1991)
APPENDIX
A: LIST OF SYMBOLS
R = radius of curvature. K(t), Ki = continuous or discrete curvature. i, k, n = integer indexes. t = parameter.
I = arc length. L = perimetral length. @ = slope angle, orientation-length function. 0 = cumulative orientation-length function. w = number of points in half a window. u = Gaussian function width. uz = (total) noise variance. (T: = quantization noise variance. m = derivative order. M = number of bits. B = signal bandwidth. F, = sampling frequency. f. = fundamental frequency. fc = cutoff frequency. fc,, = normalized cutoff frequency. N, = cutoff harmonic number. hk = filter samples. Kh = number of filter samples. G,(t) = Gaussian function. (Y= Kaiser window parameter. A = regularization parameter. 1. .I = modulus of complex function. N, = noise power. ck = complex Fourier coefficients. z, = complex contour data. c, KT, NT, eT= COWSantS. A,i, Ayi, Fxi, Fyi = amplitudes and phases of the ith sinusoidal component for rectangular coordinates of ventricle-shaped curves. e,, ,u<, Us = reconstruction errors and their mean and s.d. values.
x, y = rectangular coordinates of contours. N = number of contour points.
494
image and vision computing