- Email: [email protected]
Isotropic four-point interpolation
COMPUTER GRAPHICS AND IMAGE PROCESSING
11, 192-196 (1979)
NOTE Isotropic Four-Point Interpolation JAMES E. MIDGLEY Physics Department, University of Texas at Dallas, Richardson, Texas 75080 Received June 8, 1979 An efficient algorithm is given for constructing a smooth curve through an arbitrary sequence of data points. The method is local, using only four points at a time ; isotropic, giving a result independent of the orientation of the curve relative to the underlying coordinate system; and applicable to curves in an arbitrary number of dimensions. It makes the direction of the curve at each data point that of a circle through that point and the two adjacent data points. 1. INTRODUCTION As long as there have been tables of data there has been interest in methods of interpolation. Most such tables enumerate discrete values of a single-valued function of the independent variable(s), and in such a case piecewise polynomial approximation is appropriate and generally used. The spline functions, which preserve continuity of the first few derivatives, are the ultimate of such approximations and their theory is highly developed. Such techniques cannot be used, however, for interpolating closed curves in two dimensions or a n y other multivalued function. Recently Rutkowski E l i proposed a m e t h o d for filling the gaps in curves with missing segments utilizing cubic polynomials. With appropriate modifications and extension this technique can be used for interpolation. We will here describe the algorithm which seems to produce the best interpolation with only m o d e r a t e calculation, developing it entirely in vector n o t a t i o n so t h a t it can be applied in a n y number of dimensions. 2. DERIVATION Given a sequence of points r~, i = 0, 1 2, . . . , n, to be interpolated, define the displacements di = ri -- 1"/_1. (1) If we define a parameter t which varies from 0 at ri to 1 at r~+l, then the straight line connecting these points would be written as r(t) = rl -f- td~+l. 192 0146-664X/79/100192-05502.00/0 Copyright ~ 1979 by Academic Press, Inc. All rights of reproduction in any form reserved.
(2)
193
ISOTROPIC F O U R - P O I N T INTERPOLATION
I n a local four-point interpolation algorithm the deviation f r o m t h a t line o v e r the interval can only depend on the three vectors d~, d~+l, a n d d~+2. Assume it is a linear combination of these vectors with coefficients which are polynomials in t. These m u s t vanish at t = 0 and t = 1, a n d the derivative of the coefficient of d~+2 m u s t vanish at t = 0 as m u s t the derivative of the coefficient of di at t = 1, in order to m a k e possible the continuity of the t a n g e n t vector at the d a t a points. T h e lowest degree polynomial which can satisfy three homogeneous b o u n d a r y conditions is a cubic, so the simplest four-point interpolation formula will be of the form r(t) = ri -4- td~÷l -4- f ~ t ( 1 - - t):d~ -4- t(1 -- t)(g~ -- G#)d~+l - h i t 2 ( 1 - t)d~+~.. (3) T h e continuity of the t a n g e n t vector places one additional condition on the four constants appearing in this formula: dr - - (1) = (1
--
gi_i
G~_~)di -[-
"3t-
hi-idi+l,
dt dr
- - (0) -= fidi A- (1 A- g~)d~+l dt
so t h a t
fi
1 - - g i - 1 -4- G i _ l -
hi-,
(4)
1 + gi
Clearly the t a n g e n t at a d a t a point m u s t be a linear c o m b i n a t i o n of the displacem e n t vectors on each side of the point, a n d as the length of a vector goes to zero the contribution of the other vector m u s t correspondingly go to zero, suggesting a relationship of the form =
•
(5)
l "-[- gi
C o m p u t a t i o n a l simplicity would dictate the choice n = 2, and f o r t u n a t e l y this is the choice which makes the t a n g e n t m a t c h a circle t h r o u g h t h a t a n d the two a d j a c e n t points. T h e reader m a y confirm this using the fact t h a t the radius vector f r o m the center of such a circle to the point ri is (p A- q)d~ R =
--
2(p
( __d, ) 2, q
p(1 -4- q)d~+l , --
where
p .
q2)
.
-di+l
.
.
di" di+, .
(6)
di+l 2
Using these two conditions leaves two constants still to be determined: r(t)
= r~ +
Lt(1 -
t)~d~ +
t
at -
2t~ +
(1 -
t)V~
.~
- - t ( 1 -- t)h~ \---/d,+l ] d,+l -- h # 2 ( 1 - - t)di+,.
(7)
194
JAMES E. MIDGLEY
Rutkowski's shape completion algorithm also has two constants which must be determined arbitrarily. He characterizes them as ~1 and ~s, which in the present notation would have the values /21 :
~ (0)
= L d i I d i -t- di+ll/di+l,
---- hidi+sldi+s -t- di+iI/di+l.
(8)
In the symmetrical case (di+2 = d , Idi+2 -t- di+ll = Idi -t- di+ll), when a circular arc is the proper completion, he notes that if pl = P2 the polynomial completion will pass through the midpoint of the circular arc if 2di+1 ~1 -
, 1
(9)
.-t- cos 01
where 81 is the angle between the tangent at r~ and the vector di+1. Using the circle rule for the tangent gives di.di+l -1- di+l s cos 01 =
(10)
dildi -~- di+ll These reasonable conditions (for the symmetrical case) translate into the following values for fi and hi: 2di+l s
fi = h i - -
(11) di'di+l -t- di s -t- dildi -t- di+ll
Things get a bit more arbitrary, though, when a general form is chosen, applicable to the nonsymmetrical case as well. It should be invariant to the reversal of the sequence of numbering, and it can be shown t h a t this will be true if and only if f~(d~, d~+l, di+s) = h~(--di+s, --d~+l, - d i ) .
(12)
E q u a t i o n (11) satisfies this condition only in the symmetrical case. Rutkowski chooses as his general relation 4di+1 72 1
=
lP s
(13)
=
2 -t- cos 01 -t- cos 0s which in the current context leads to the relations
fi = 4 / ( p + q + r(2 -t- (u + v)/w)), hi = 4 / ( u -t- v + w(2 + (p + q)/r)), where
p = diS/di+l 2,
q = di.di+l/di+l 2,
u = di+s2/di+l s,
v = di+s'di+l/di+l s,
r = (p2 + 2pq + p)i, w=
(u s + 2 u v + u ) ~ .
(14)
ISOTROPIC FOUR-POINT INTERPOLATION
195
These satisfy (12) and in the s y m m e t r i c a l case (di+2 = 2qdi+l - d i ) reduce to (11), b u t t h e y are not the simplest forms which do so. I n s t e a d we suggest the choice : fi = 2 / ( p - ~ - q + r),
h,i = 2 / ( u + v + w).
(15)
N o t only is there no simpler form which will satisfy the criteria, b u t of the m a n y a l t e r n a t i v e forms tried none gave a n y s m o o t h e r curves a n d m o s t gave obviously u n a c c e p t a b l e curves. R e a r r a n g i n g (7) for efficient calculation, the final interpolation formula is as follows : r(t) = r~ + At + Bt 2 + Ct a, (16) A = f i d i -~- fipdi+l, B = -2fidi
+ (3 -
2f,.p -
hiu)di+l -- hid,+2,
C = f~di ~- ( - - 2 + f~p + h~u)d~+~ + hidi+2. 3. END POINT INTERPOLATION Since do is undefined, (16) cannot be used to interpolate between r0 a n d rl without a further assumption. I t would seem best o v e r this interval to approxim a t e the circle which passes t h r o u g h the points r0, r,, a n d r2, a n d t h a t is just w h a t (16) will do if we define dl"d2~
do = 2 - - - - - d l - - d2. \ dl 2 /
(17)
Likewise (16) can be used over the final interval if we define (d.-d._~ d.+,
= 2
---
\
d, 2
/
d. -- d._1.
(lS)
4. PRACTICAL IMPLEMENTATION W h e n the purpose of the interpolation is to produce points for a vector p l o t t e r which will be close enough together to produce the a p p e a r a n c e of a s m o o t h curve, the absolute spacing needed will depend on the resolution of the plotter. U n d e r some circumstances it m a y be desirable to m a k e all the elemental vectors of about the same length As, in which case the i n c r e m e n t in t should be As At --
As =
{dr/dtl
(19) [A + 2Bt + 3Ct 2]
Under other circumstances it m a y be desirable to limit the angle between a d j a c e n t vectors to some m a x i m u m A0, in which case the i n c r e m e n t in t should be {A q- 2Bt -4- 3Ct~l At _<
A0.
2{2B + 6Ctl
(20)
196
JAMES E. MIDGLEY
FIG. 1. Smooth interpolation of a curve specified by only 15 points. The simplest algorithm, Eqs. (15) and (16), yields the solid curve; the more complicated algorithm, (14) and (16), yields the dotted curve. Usually, however, a sufficiently smooth curve is generated most economically b y the geometric mean of these two criteria:
)0
At =
•
(21)
rB + 3Ctl This is evaluated only once, at t = 0 or t = 1 (whichever makes it smaller), and used over the whole of the interval. Figure 1 illustrates the use of the m e t h o d to generate a rather complex curve specified b y only 15 data points, indicated b y tick marks on the curve. T h e solid curve was generated using (15) for f and h, while the dotted curve utilized (14); they are remarkable for their s i m i l a r i t y - - o t h e r choices gave rather different curves. The dotted curve also illustrates the varying vector size from the use of (21) ; the value of AOAOused there was 0.03. REFERENCE 1. W. S. Rutkowski, Shape completion, Computer Graphics and Image Processing 9, 1979, 89-101.
11, 192-196 (1979)
NOTE Isotropic Four-Point Interpolation JAMES E. MIDGLEY Physics Department, University of Texas at Dallas, Richardson, Texas 75080 Received June 8, 1979 An efficient algorithm is given for constructing a smooth curve through an arbitrary sequence of data points. The method is local, using only four points at a time ; isotropic, giving a result independent of the orientation of the curve relative to the underlying coordinate system; and applicable to curves in an arbitrary number of dimensions. It makes the direction of the curve at each data point that of a circle through that point and the two adjacent data points. 1. INTRODUCTION As long as there have been tables of data there has been interest in methods of interpolation. Most such tables enumerate discrete values of a single-valued function of the independent variable(s), and in such a case piecewise polynomial approximation is appropriate and generally used. The spline functions, which preserve continuity of the first few derivatives, are the ultimate of such approximations and their theory is highly developed. Such techniques cannot be used, however, for interpolating closed curves in two dimensions or a n y other multivalued function. Recently Rutkowski E l i proposed a m e t h o d for filling the gaps in curves with missing segments utilizing cubic polynomials. With appropriate modifications and extension this technique can be used for interpolation. We will here describe the algorithm which seems to produce the best interpolation with only m o d e r a t e calculation, developing it entirely in vector n o t a t i o n so t h a t it can be applied in a n y number of dimensions. 2. DERIVATION Given a sequence of points r~, i = 0, 1 2, . . . , n, to be interpolated, define the displacements di = ri -- 1"/_1. (1) If we define a parameter t which varies from 0 at ri to 1 at r~+l, then the straight line connecting these points would be written as r(t) = rl -f- td~+l. 192 0146-664X/79/100192-05502.00/0 Copyright ~ 1979 by Academic Press, Inc. All rights of reproduction in any form reserved.
(2)
193
ISOTROPIC F O U R - P O I N T INTERPOLATION
I n a local four-point interpolation algorithm the deviation f r o m t h a t line o v e r the interval can only depend on the three vectors d~, d~+l, a n d d~+2. Assume it is a linear combination of these vectors with coefficients which are polynomials in t. These m u s t vanish at t = 0 and t = 1, a n d the derivative of the coefficient of d~+2 m u s t vanish at t = 0 as m u s t the derivative of the coefficient of di at t = 1, in order to m a k e possible the continuity of the t a n g e n t vector at the d a t a points. T h e lowest degree polynomial which can satisfy three homogeneous b o u n d a r y conditions is a cubic, so the simplest four-point interpolation formula will be of the form r(t) = ri -4- td~÷l -4- f ~ t ( 1 - - t):d~ -4- t(1 -- t)(g~ -- G#)d~+l - h i t 2 ( 1 - t)d~+~.. (3) T h e continuity of the t a n g e n t vector places one additional condition on the four constants appearing in this formula: dr - - (1) = (1
--
gi_i
G~_~)di -[-
"3t-
hi-idi+l,
dt dr
- - (0) -= fidi A- (1 A- g~)d~+l dt
so t h a t
fi
1 - - g i - 1 -4- G i _ l -
hi-,
(4)
1 + gi
Clearly the t a n g e n t at a d a t a point m u s t be a linear c o m b i n a t i o n of the displacem e n t vectors on each side of the point, a n d as the length of a vector goes to zero the contribution of the other vector m u s t correspondingly go to zero, suggesting a relationship of the form =
•
(5)
l "-[- gi
C o m p u t a t i o n a l simplicity would dictate the choice n = 2, and f o r t u n a t e l y this is the choice which makes the t a n g e n t m a t c h a circle t h r o u g h t h a t a n d the two a d j a c e n t points. T h e reader m a y confirm this using the fact t h a t the radius vector f r o m the center of such a circle to the point ri is (p A- q)d~ R =
--
2(p
( __d, ) 2, q
p(1 -4- q)d~+l , --
where
p .
q2)
.
-di+l
.
.
di" di+, .
(6)
di+l 2
Using these two conditions leaves two constants still to be determined: r(t)
= r~ +
Lt(1 -
t)~d~ +
t
at -
2t~ +
(1 -
t)V~
.~
- - t ( 1 -- t)h~ \---/d,+l ] d,+l -- h # 2 ( 1 - - t)di+,.
(7)
194
JAMES E. MIDGLEY
Rutkowski's shape completion algorithm also has two constants which must be determined arbitrarily. He characterizes them as ~1 and ~s, which in the present notation would have the values /21 :
~ (0)
= L d i I d i -t- di+ll/di+l,
---- hidi+sldi+s -t- di+iI/di+l.
(8)
In the symmetrical case (di+2 = d , Idi+2 -t- di+ll = Idi -t- di+ll), when a circular arc is the proper completion, he notes that if pl = P2 the polynomial completion will pass through the midpoint of the circular arc if 2di+1 ~1 -
, 1
(9)
.-t- cos 01
where 81 is the angle between the tangent at r~ and the vector di+1. Using the circle rule for the tangent gives di.di+l -1- di+l s cos 01 =
(10)
dildi -~- di+ll These reasonable conditions (for the symmetrical case) translate into the following values for fi and hi: 2di+l s
fi = h i - -
(11) di'di+l -t- di s -t- dildi -t- di+ll
Things get a bit more arbitrary, though, when a general form is chosen, applicable to the nonsymmetrical case as well. It should be invariant to the reversal of the sequence of numbering, and it can be shown t h a t this will be true if and only if f~(d~, d~+l, di+s) = h~(--di+s, --d~+l, - d i ) .
(12)
E q u a t i o n (11) satisfies this condition only in the symmetrical case. Rutkowski chooses as his general relation 4di+1 72 1
=
lP s
(13)
=
2 -t- cos 01 -t- cos 0s which in the current context leads to the relations
fi = 4 / ( p + q + r(2 -t- (u + v)/w)), hi = 4 / ( u -t- v + w(2 + (p + q)/r)), where
p = diS/di+l 2,
q = di.di+l/di+l 2,
u = di+s2/di+l s,
v = di+s'di+l/di+l s,
r = (p2 + 2pq + p)i, w=
(u s + 2 u v + u ) ~ .
(14)
ISOTROPIC FOUR-POINT INTERPOLATION
195
These satisfy (12) and in the s y m m e t r i c a l case (di+2 = 2qdi+l - d i ) reduce to (11), b u t t h e y are not the simplest forms which do so. I n s t e a d we suggest the choice : fi = 2 / ( p - ~ - q + r),
h,i = 2 / ( u + v + w).
(15)
N o t only is there no simpler form which will satisfy the criteria, b u t of the m a n y a l t e r n a t i v e forms tried none gave a n y s m o o t h e r curves a n d m o s t gave obviously u n a c c e p t a b l e curves. R e a r r a n g i n g (7) for efficient calculation, the final interpolation formula is as follows : r(t) = r~ + At + Bt 2 + Ct a, (16) A = f i d i -~- fipdi+l, B = -2fidi
+ (3 -
2f,.p -
hiu)di+l -- hid,+2,
C = f~di ~- ( - - 2 + f~p + h~u)d~+~ + hidi+2. 3. END POINT INTERPOLATION Since do is undefined, (16) cannot be used to interpolate between r0 a n d rl without a further assumption. I t would seem best o v e r this interval to approxim a t e the circle which passes t h r o u g h the points r0, r,, a n d r2, a n d t h a t is just w h a t (16) will do if we define dl"d2~
do = 2 - - - - - d l - - d2. \ dl 2 /
(17)
Likewise (16) can be used over the final interval if we define (d.-d._~ d.+,
= 2
---
\
d, 2
/
d. -- d._1.
(lS)
4. PRACTICAL IMPLEMENTATION W h e n the purpose of the interpolation is to produce points for a vector p l o t t e r which will be close enough together to produce the a p p e a r a n c e of a s m o o t h curve, the absolute spacing needed will depend on the resolution of the plotter. U n d e r some circumstances it m a y be desirable to m a k e all the elemental vectors of about the same length As, in which case the i n c r e m e n t in t should be As At --
As =
{dr/dtl
(19) [A + 2Bt + 3Ct 2]
Under other circumstances it m a y be desirable to limit the angle between a d j a c e n t vectors to some m a x i m u m A0, in which case the i n c r e m e n t in t should be {A q- 2Bt -4- 3Ct~l At _<
A0.
2{2B + 6Ctl
(20)
196
JAMES E. MIDGLEY
FIG. 1. Smooth interpolation of a curve specified by only 15 points. The simplest algorithm, Eqs. (15) and (16), yields the solid curve; the more complicated algorithm, (14) and (16), yields the dotted curve. Usually, however, a sufficiently smooth curve is generated most economically b y the geometric mean of these two criteria:
)0
At =
•
(21)
rB + 3Ctl This is evaluated only once, at t = 0 or t = 1 (whichever makes it smaller), and used over the whole of the interval. Figure 1 illustrates the use of the m e t h o d to generate a rather complex curve specified b y only 15 data points, indicated b y tick marks on the curve. T h e solid curve was generated using (15) for f and h, while the dotted curve utilized (14); they are remarkable for their s i m i l a r i t y - - o t h e r choices gave rather different curves. The dotted curve also illustrates the varying vector size from the use of (21) ; the value of AOAOused there was 0.03. REFERENCE 1. W. S. Rutkowski, Shape completion, Computer Graphics and Image Processing 9, 1979, 89-101.