- Email: [email protected]
Fitting conic sections to scattered data
COILIPL~TI~:It
OIt.Ll?‘HICS.\ND IKLGIC PItOCIGSSIKG
Fitting
Conic
9, 56-71 (1979)
Sections
to Scattered
Data
FRED L. BOOKSTEIN Center for Human Growth and Biostatistics, The University
Development of Michigan,
and Departments of Statistics Ann Arbor, Michigan @109
and
Received January 6, 1978; revised March 28, 1978 The problem of fitting conic sections to scattered data has arisen in several applied literatures. The quadratic form Ax2 + Bxy + Cy2 + Dx + Eu + F that is minimized in mean-square is proportional to the ratio of two squared distances along rays through the center of a conic. Considerations of invariance under translation, rotation, and scaling of the data configuration lead to a straightforward method of estimation somewhat different from earlier suggestions. The method permits an extension to conic splines around extended digitized curves, expediting a smooth reconstruction of their curvature. Some examples are presented indicating how the technique might be applied in morphometrics. 1. INTILOI~UCTION
In the course of investigations into the shape of plane curves, several authors, attempting to fit conic se&ions to scatters, have recently and independently discovered algorithms that arc reducible to conventional linear computations. All minimize the sum-of-squares of a form Q(x, y) = AJ~ + Rxy + Cy2 + Dx + Ey + F over a scatter in the (z, y)-plane. Paton [l-2] sets A2 + B2 + C2 + D2 + E2 + F2 = 1 as the constraint on this minimization, whereupon the solution vector (A, B, C, D, E, F) is a principal component of the cross-product matrix of the variables x2, zy, y2, 2, y, 1. Biggerstaff [S], Albano [4], and Cooper and Yalabik [5] set F = 1 as the constraint, whereupon there arises a set of normal equations for A, B, C, D, E in terms of the same cross-product matrix. Gnanadesikan [S] set’s A2 + B2 + C2 + D2 + E2 = 1, yielding for extrcmum the last principal component of the variance-covariancc matrix of x2, my, y2, 2, y. All these authors loosely characterize the optimal conic as “best-fitting” of all the tonics of the plane, but the nature of the fit each is optimizing, whose computation is so surprisingly straightforward, has not been elucidated. WC seek an estimation rule which is general, simple to compute, and invariant. Generality requires a single algorithm by which we fit both ellipses and hyperbolas, with circles and parabolas included as constrained special cases. For instance, we cannot generalize the fitting of parabolas usual in quadratic regression, in which error-of-fit, is measured parallel to the axis of the parabola, because the
0146~GG4X/79/010056-16$02.00/O Copyright All rights
0 1970 by Academic Press, Inc. of reproduction in any form reserred.
Iti\~ariatwc is to lx \vitli rcq)c~c~t~ to tllc> c~cjuiform h~~roup of transi’ormatioits oi’ t Ii(b l~~u~litloati platir : rotations, t~ratislatiotis, and c~liattgc~sof scalp. if tlio wortlito an ~~l~~m~wloi’ tiutctr of t’lici data aw trattsformc~d, point 1)~ point, wwrding this group, tlwti tlw kst-fitting coriic~ to the ti~\v wnttc>r must 1)~sc~~at*tIy t tics ws~ill of t.lic same transformation applicstl to tlica cwtiich I\-lticli ~2s hwt-fitting to tlirl original scatter. ‘l?i~~wtiniatiott rules proposc~d lwr(~ mitiimizw a cluadratic* form in tlw utiktio\vn cotticd l)aranic$t~ri sut)jwt to ati w:wt wrlstraiiit on tlw val\tct of anothc~r cluadratic form. It \\ill 1~~sIw\vtt t 1~~1Iwth tlicw forms arcs invariants und(ar the action of th(a c~c~ttit’ornlgrortl), and lwtrw so is tlw tic4 wtimstion proccdurc~. To insist on this invariant is to alickr our intttitivc notion of curvesfitt,ittg “to a dcpondwit varisl,l(s.” \\.ltcsti al~wiss:~ and ordinatc~ wpwwtit 2111 arbit,rary coordinate q~stc~m, or arbitrary oriotttatiott. plapcbd mpor~ an undwlyitt~; gc~omc~trical ol)jxt, swh us tlicl outliirc~ 01’ a skull, t It(x axcbs wparatrbl!. liavo no gc~omc%ric mclaning at all; tlw~. arc’ \\ Itoll~~ (‘Oliilil(‘tisitr;Lt(’ ant1 may lw rotattd frwl~~. WC mrtst gclt)tlw sam( j gwnic~tric~ rcxsttlt iri all (2~s. 2. (:I~:OMI~;Tl:I(’
IN’~I~:I:l’l~l~:T.21’IO?j
01; Q(.r, !/)
1~~1a giwn wnt,ral conic wctioti 11k~vc~ c~cluatiott c)b, y) = .l~:’ + lj.r,y + C’,$ + 11.r + 1Q + E’ = 0. I,(+ (xl,, yO) 1)~~an arbitrary Jwint. ot’ tltct plaw, and s(‘~ (xl, ~1) qua1 to the point on 11~~1 conic on the ray from thtl wntw through (.I’,,, !I,~). Cltoosc: a IWW coordinat.cl systclm, d(wotcld 1)~ primw, n-hic*h is wntercd at t ll(t cotitc>r of tlw conic and aligrictd with its priwipal aws. J>~finc~tlicl form Q’(.r’, !/‘:I on tSlrcbrw\\- coordinates so that Q’(/ , !/‘j = ()(.I*, y) OII tlic l)lanc (A, y’). In tlti!; IWK coordinate systrm, Q’(.r’ y’) = :l’.P + 13’jg’2+ 1.” for somchA’, B’, F’. In the primc\d systcsm! the> wntc~r is thcb origin. (,(I’, 0’). By c~olliwarit\~ OI (,r,,, g,,), (x1, ,ql), and tho ccwtor, (bxpr(wed in the, prim(bd cwordinat.c> s;!stc~m. UC’ Ilaw d.‘~‘!/‘I = do ‘~‘0 or .r’,) = .r’] t!/‘,, !/‘I). Th(~ti
-
/“‘[rJ/‘,, ‘!/‘,)”
I j,
i II
since (c’~, g’lj is on the conic C)‘(J’ , y’) = 0. But the term (g’Cv’g’l) is just tllcbratio of distances writer-to-point and c,c:itter-to-c~otti(~don p the ray JVCarc using -tlw
5s
FlU3:D
L.
BOOKSTEIN
Ellipse
Parobofa FIG. 1. The geometry of error-of-fit conic. The error is always proportional
Hyperbola around a conic section, to [(d + d$/#] -
as measured 1.
by the
equation
of the
ratio of d + dI to cl in Fig. 1; for projection onto any axis preserves distance along a line. Hence &(x0, 1~01= [I (d + dJ2/C121 - 1, and as the equiform transformations do not alter distance ratios they preserve this proportionality (though they may alter F’). We have [(cl + d$ - cP]/cP = &(cZI + 2d)/d2. When the data are quite close to a conic, cl1 is small, dl $ 2d is approximately equal to 2d, and Q(z,,, y,,) is approximately proportional to dl/d, the distance from the conic along a line through the conic’s center, measured in units of distance from the center to the conic in that direction. The fitting of a parabola is a limiting case, exactly transitional between ellipse and hyperbola. As the center of an ellipse moves off toward infinity while its major axis and the curvature of one end are held constant, the rays from its center become parallel, and the difference of squared distances which is the numerator of Q becomes linear in distance. The quantity being minimized in the fitting then approaches the usual deviation, Euclidean distance from the curve measured parallel to the axis.
For a circle, I propose using the instance of least mean-squared Q. Let there be given points (~1, yl), . . ., (xn, yn) and let the sample linear regression of x2 + y2 on x, y, 1 be x2 + y2 - 2ax + 2by + c, where a, b, c are the least-squares
4.
IN\‘.\I:L.\N’I’
E‘OI:
‘l’lll~:
(;I,:NI~:I:.\I,
C’ozIt
‘
1~01~an ahitrary conic, tlic gxwc~ralizatiott of Ilit> cwor-of-fit i.r - 0’):: + (y - b)’ - A? is clwrl\tlitr valnc~ (J(,.i’, y) = .I tl hy miriinnizing CGrLB,., (i/d. + I)!/ + (.jz sul)j~.~~tto tlw cwnstraitll (1” + /I” = 1. For conics, w: might tq’, IJ~ clirwt niimict~y, the, for~rl -1’ + 1:’ + C”. I’IIfortutIatc>ly, the liypcrbola Lh.y = 1, with .t” + 13’ + f’” = 1, can lx, rotated into th form .J.‘)- p2 = 1, with A’ + 11” + C2 = 2, and so tit,, first would automal imlly fit t\viw as l~adly if tltc data w(‘rc rotat cd 15”; hut. tltcs sclccted normalization tnust be invariant, utidw rotation. Ititroduc~ti~~w lo tlw algel~raic tliwry of conic85 tell us that the forms A + (’ ant1 11’ --- .4.1C aw invariant iitickr the I3tclideati Sgrou~),l’hc~ onl\- r)ositivc~-tl(~finit(~itiv:tri;tlit t Ilat (%ti l~c’ formcbtl from tlrwc cpin~~t i-
tics is (fl + 6)” + (B2 - UC)/2 = A” + B2/2 + C2. We may set a unique scale for Q, and thus proceed to actual minimization of ZQ2, b*y setting this norm equal to some conventional constant value. I suggrst the value 2, so that the equation of a circle is in the usual form. Given data points (51, yl), . . . , (z~, yn) scattcwd about a conic, I propose to estimate the conic of best fit by minimizing C;[G,(zi, yJ12 subject to the constraint VDV’ = 2, where V = (A, B, C, D, E, F) is the vector of coefficients to be estimat.ed, and D is the matrix diag (1, 4, 1, 0, 0, 0). The minimand may be written VSV’, where S is the scatter matrix about 0 of the 6-vector (z?, xiyi, yi2, z;, y;, 1) as a function of the data points (x;, vi). The invariance of this procedure is proved straightforwardly. If the coordinates of the data points are rotated or translated, the constraint A2 + B2/2 + C2 is unchanged, likewise the term F’ of Eq. (I), which is gotten by undoing the transformation; and the distance ratio, the other multiplicand of (l), is clearly unchanged. Then the fitting is invariant under translation and rotation. If the coordinates are resealed by k, (2, y) + (x’, y’) = (kx, Icy), let V, be the vector of coefficients of the resulting best-fitting conic by this algorithm. The conic As2 + . . . + F = 0, after resealing and renormalizing, becomes Ad2 + Bz’y’ + Cyf2 + kDx’ + kEy’ + k2F = 0. The resealing has replaced the matrix S by SI, = DkSDb with Dk = diag (k2, k2, k2, Ic, lc, 1). Note that DkDDk = k4D. Extremization of Vk(DkSDk)Vk’ for V,DVk constant is equivalent t’o extremizing it for VkDkDDkVk’ constant. Since Dk is nonsingular, the extremum is V/Dk, where V is the extremum before resealing. But V/Dk gives the coefficients of the old conic in the new system, for k2V/D~ = (A, B, C, kD, kE, k2F). Hence this method of fitting is invariant under equiform transformations. As the authors cited in Section 1 all use D and E in the norms for their fitting, none of their solutions are invariant under the Euclidean group. Those which involve F, further, are not even invariant under translation : For by setting F = 1 it becomes impossible to fit tonics through the origin (F = 0) at all. These authors must therefore forcibly standardize their data before conic-fitt’ing in any of various ad hoc ways. Translation to center-of-mass and rotation to principal axes of scatter are the most common. This approach seems unsatisfactory. Dat’a must be fit in a manner invariant under t’he Euclidean group, not arbit’rarily constrained with respect to it. 5. ESTIMATION
We wish to minimize
(VI 1V,), both components
FOR
THE
GENERAL
CONIC
VSV’ subject to VDV’
= constant. of length 3, and let S bc partitioned
Partition V into correspondingly :
SC(+ I )* x12
21
VSV’ = VlSllV’l
+
x22
2v1s,,v’2
+
v2s22v2.
We must minimize this subject to VIDIV’ 1 = constant, where DI = diag (1, 8, l)*
‘I‘llis
i9mputution
upon
1 tw
cdimat
ing,
I\-hiw
M
anil \vi’
=
matrix niinus (‘.gL, ktraint Sctt \\:lys,
a wt
is a matrix,
cdrcmllm 1-D 1”
remains
parami+~r5
scc~k:
-
M
iliwotiss it, must h may ing s~)niil
of
iwnstraint Hao
l’hc
vcdor
I’M
=
0
rlwil [S,
tlli> -I)
ortlinaty-
\vib plains
tlc~noti~
tlw
lc.(i]
instruct,5
smallid
iilvc~r,q~. invcw~.
\vitll ‘l’lic~
of riqwt
to
mat,ris D
in IIS
T-S\”
\\ Irtww
linear of
togi+twr
i~igi~nvcdor D
arbitrary
Ci-wctor
\\-ritt,csrl
Swtion T- of
is
gi~ni~raliaid an
\\.lic~n to
,q is
(M’S--‘M)-M’S any
haw
I- is
t81wn
wwdant, (I
t ractal)l~~
If
the in this
T’M
formula
subjcd
to
c~igi~nvallw
S.
Hew
wit
form
largwt, S mud
ma!.
constraints
i~oiBkii~nts
the
for
tlw
of
tllc,
supcwcript, clcfiniti~.
and
thi>
CY)II-
1x3 rdunda~~t.
linear ilriitit
const~raints usdul.
crodw IIc~w
:iw
this wn1i’
avuilal,li~
diyycw
I,iwsil)ilitiw
:
of
friwlom
0, thi.
constraints
1~’ positive
1)~s singular
arc’ =
in
various
. FIG.
2.
Two tonics abutting
at a. knot where the tangent crosses the curve.
(5) 2Aso + By, + D = 0, Bxo + 2Cy0 + E = 0. The conic has its center at (x0, ~0). (6) Ar? + Br~yl + Cy,Z + D.z:~+ Eyl + P = 0. The conic goes through (Xl,
Yl).
(7) @Ax:1 + BYI + 0) cos cr + (B.xI + 2Cy1 + &) sin o( = 0. If the point (x1, yl) is on the conic, this is the criterion that the tangent there make angle LY with the positive x axis. (8) A (- 21 sin LY+ R cos2a) + B (-05y1 sin2 a! + 0.521 cos2 o( + R cos a sin a) + C(y, cos c11+ R sin2 a) + D( -0.5 sin a) + E(0.5 cos a) = 0. The conic through (x1, yl) with tangent making an angle cr with the positive x axis is by this constraint required to have radius of curvature there equal to R. 7. CONIC
SPLINING
There is a rapidly growing literature (cf. Poirier [9], Wold [lo]) on polynomial spline regression. In this technique, the expected value of a dependent variable is fitted by different functions-usually polynomials-over distinct ranges of an independent variable. At ‘(knots,” where these ranges abut, the predictions are required to agree in value and in their first few derivatives. After the knots are set, all the polynomials are computed at once, by minimizing the tot’al squared error-of-fit about all of them in one clever dummy-variable regression. The method becomes much more complex when error-of-fit is not always measured along a fixed axis. The constraints which force the curves to line up become nonlinear in the coefficients being estimated, and linear methods fail (cf. Thomas [ll] and references therein). In particular, one cannot spline a closed curve linearly by polynomials. But just as parabolas (and higher-order polynomial graphs) generalize straight-line regression, so do tonics (and higherorder algebraic curves) generalize st’raight-line curve-fitting in general orientation. It-is then of interest to explore the possibility of fitting multiple tonics to separate parts of a curve simultaneously in a single linear computation generalizing the previous discussion for fitting tonics singly. , Suppose we have a data set in the plane which looks like two tonics, in the manner of Fig. 2, impinging upon each other at some visible nexus P. It would be nice to fit two conic segments simultaneously, one to the left of P, one to the right, which minimize the net error-of-fit. We would like the tonics to connect at P without any corner, that is, to have the same tangent there. The inclination
ITlwN(:
(10X1(’
S1~:~‘1’10NS ‘1‘0 s(‘.~‘1”1’I~:l: I,:I) I ).\‘I’.1
(it:
is not givcw ill IIu: flata, Iio\vc~vc~r; \\‘I: n11bd f~sliumtc it, I oo, as that position of 1lit: langc~nt lint> at I’ for whkh t ht\ bc~sttwo conic?; through I’ \vith that tangcbnt, 011~: I’or f~af*lk sid(x, togcstIiw huvc~ lfmt nlcan-squarfd fvorof-fit. To fbxf>cutc this plan \vf’ must link thfl tlormalizations of thf: tmm conic+. The bwt way to do this is by insisting that, a givctn fxrror-of-fit near I’ 1~ tlio Sam0 I\-hichww caonic WCconsidw it on. To a linclar approximation, thcl error for a point such a:: R rclatiw to a conic (r)(s, IJ) = 0 through 1’ is c~lual to t,hc>gradicant of the walsr fif4d CI)(.r!y) timw thcl distant from I’ to 22 nornml to tlifb fs~irvf’. ‘I‘hcw~ \vill bf> (qua1 for t,licBt\\-0 f~oriks impinging at I’ if’ and only- if the, gratlicnts of t,hc> t ~1-0 scalar fields arc ftflual at I’. But to spccif~. that thwl two gradifsnts l)fb c~111alis jwt a [lair of’ linfbar COIIst raints! for thf: gradients are liwar forms in th(l cwffic4c~nts of the cork. Spc’ificdl\- : if thf> two conies arc’ (Jzi.r. y) = -1 ,.? + B,J!/ + C1!/ll + 11,~ + E,y + P = 0. cI’. s): thcxu thcl t ii-0 con&tc~nc~!~crit(~ria :WI I’ = 1, 2: and if P have t hc coordinate
0I t,liis mul)ual tangfwi,
“:l 11’+ up + D, -- (2.1 .‘t’ + B,s + II,) = 0, 3c,s + n,r + I:‘, -
iX’,s + II,,, + 7:‘,) = 0.
IVe must, also insist, that both caonics actually pa.hs through thr point (I,: s) at ivhich KC:arc specifying t,hAr tangents, and this is t\vo mow constraints :
..w + R,,,~+ c,2 + rjzr+ E:‘,.~ + I?,= 0,
I’ = I,?.
\\Tc can estimate both conks, all 12 cof+Ikicnts, simultanc~ously by minimizing their total squared error-of-fit subject to t hfi normalization -4,’ -+ HIZi’2 + C1’ = 2 and th(h four constraints above. This tact’ic is easily gcncralized for thft fitting of conic spliws to cxxtcndcd chain?, cwn closed curves. Let thaw be I arcs, and on the it,h arc data points (I,,, vi,). ,j = 1. . . ., f/,, 7, = 1, . . ., I. Let the arcs abut at knots (.fi: !I,):, i = 1, . . . , I, and sot (II. 1, yf+l) = (./.I, yl). Let) Si he thf: G X Ci scatter matrix about 0 for tlift points (Xi,‘, sX,j,!Jij,uij 2, T,j, !/aj, 1) dflrivctl from the data points (s,~, ?/lj) of thfn it,h arc. Define S, (61) X (GI), equal to blockdiag (A’,, . . . , SI). Consider a conic* spline which is equal to AL.? + R,.fy + C,y’ + D,z + E,y + I;‘i on the ith arc, i=l I, and let. Al, 1 = .-11, . . . . PI .m1= PI. The error sum-of-squarf~s al)oilt the itll’ki for the data of the I’th arc is (.A;, 13,, C,, I),, I:‘,, F,)AS’,(-I,. Z3,, C’,. II,. I!‘;, Pi)‘. Dcfinc 1’ = (AI. . . . , k’lJ A,. . . . , I”,-1, rlr, . . . , F,), a ((il)-vwtor. l’hf~~ the tot$al sum-of-squares for wror-of-fit about all tlicx arcs is 1 JI-‘. This may I)(, minimized subject to a set of homogcncous linear caowtraint:: 1-M = 0 by the formula of Rao quoted above. In t5his application, S must lw nonsingular ; thcsrl each S, must, be nonsingular, which implies that 71, > 5, all i, and that thr point:: of all arcs must lie in “general position” (not all upon on< ronip).
G4
FlmIl
L.
BOOKSTEIN
and also subject to a normalization. As the tangency constraints explicitly link normalizations of one conic t,o the next, we can only assign one normalizing constraint for any subarc chain of continuously turning tangent. In the usual case where this subarc is the whole spline, as in Figs. 3, 5, 6, and 7, this means that any constraint Ai + Bi2/2 + Ci2 = constant d&ermines the normalization of all the arcs all the way around. 8. WORKED
EXAMPLE
The vaguely biological outline, Fig. 3a, derives from a freehand sketch on graph paper. I chose three “landmarks” : No. 1, the “corner” in lower center; No. 2, the extreme right point; and T\‘o. 3, the ext’reme top point. The word “landmark” is used here, in preference to the “knot” of the splining literature, t’o suggest determination by anatomical criteria rather than asympt.otic optimization tactics. If this were a skull X ray, Ko. 1 might be the audit,ory mcafus, No. 2 the bridge of the nose, and No. 3 the bregma. Points between the landmarks have been digitized crudely by cyc. The spline arcs are computed by a program executing the following algorithm. (a) The data are read in arc by arc and the scatter matrix S accumulated and stored. A vector of 18 unknown coefficienOs V = (Al, . . . ) F1! . . . ) Aa, . . , F3) is allocated, bearing the conic coefficients arc-wise. (b) The constraints are written out in full. Arcs 3 and 1 pass through point 1 of coordinates (0, -3). This requires
Ai.
+ Bi.0 + Ci.9 + Di.0 + Ei.(-3)
+ F, = 0,
i = 3, 1,
or
Similarly,
9Cs - 3E3+FJ = 0,
(1)
9C1- 3E1+ FI = 0.
(2)
arcs 1 and 2 pass through point 2 = (8.5, 0), requiring
and for point 3 = (-2,
At landmark are no further constraint and first conic and
72.2581+ 8.50, + F1 = 0,
(3)
72.25A, + 8.50, + Fz = 0;
(4)
5.7),
4Az - 11.4Bz+ 32.49cz - 202 + 5.732 + Fz = 0,
(5)
4A3 - 11.4B3 + 32.49Cs - 203+5.7Es+Fs
(‘3)
= 0.
1 we allow the arcs to int’ersect at an arbitrary angle, and there constraints to be had here. At landmark 2, we pick up a seventh an eighth from componentwise equality of the gradients of the the second: 17A1+
Dl-
17Az - Dz = 0,
8.5B1+ El - 8.5Bz - Ez = 0.
(7) (9
FITTTKC;
.
(:ONI(’
.
SWTlON~
TO ,S(‘.4TT1~:l:I~:I)
l).~T.\
. .
. . . .
a
FIG. 3. (a) Hypothetical data set. (b) 4 three-arc ronic spline approxim:ltinp knot,s at the points numbered in (a).
The final two constraints
-&I,
this scAtw, with
derive from point X similarI).:
+ 5.7Ba + D:! + 1A:j - 5.7B3 - 1)s = O>
11.4C2 - 2Bz + E’z - 11.4C3 + 2B3 - E, = 0.
(!I) (10)
(c) For any vector 17 of coefficients for any spline, the total error-of-fit is the scalar VSV’, where S is the cross-product matrix accumulated in (a). WP
GG
FRED
L.
BOOKSTEIN
nocd to niinimizo this for fixed invariant norm 111~+ Bi’/2 + Cl2 = VD V’, D = diag(l,$,l,O, . ..) 0), subject to constraints (1) through (10). Conversely, we can instead maximize VDV’ for fixed error VSV’ subject to the same 10 constraints. The theorem of Rao quoted previously applied verbatim here to provide a computing formula for the solution wc seek. A single eigenvector of a certain large matrix manifests the coefficients of all three tonics, in blocks of six, for its loadings. The spline corresponding to this eigenvector exactly passes through all the knot,s with well-defined tangents there and minimizes the total error-of-fit from the tonics in between. This curve is conic between every pair of consecutive landmarks, and its tangent turns continuously throughout, except where we have expressly omitted t,he constraint that it do so. (d) These arcs are drawn out between appropriate endpoints (by short segmental approximations) and then superimposed over the data in Fig. 3b. The spline fits the data quite well except on the far left, where the conic arc has only one curvature maximum whereas the data apparently have two (the northwest and southwest “corners”). This systematic deviation can, if we wish, be declared a “local feature” of the outline relative to the spline, as the nose is a local feature relative to the plane of the face ; it is not necessarily a flaw in the fitting. The splining-with-landmarks replaces the observed (digitized) boundary with smooth conventional arcs. Its tangent and normal, turning smoothly, provide a moving coordinate system to relate deviations from the “standard” in different regions all around the curve. The spline is more elegant and more easily read than a detailed tracing, as local features do not disrupt the global aspects of the shape, the ebb and flow of curvature around the form. Should a landmark, say the second, bc itself known only with error, one can free the tonics from passing exact’ly through it. The revised equations would require, instead of Eqs. (3), (4), (7), and (8), only that the landmark have a common distance from the two tonics and the same polar with respect to them. The appropriate constraint equation replaces (3) and (4) with the equality of their left-hand sides : 72.2581 + 8.501 + F1 - 72.2582 - 8.502 - Fz = 0;
Eqs. (7) and (8) persist unchanged. With this revision, the landmarks serve only a housekeeping function as knots, and do not serve any special role as data themselves. 9. RECTIFIABLE
CURVES
AND
MORPHOMETRIC
DATA
ANALYSIS
The technique just presented applies directly to a nagging problem in morphometrics: the function representation of outline. In morphometrics t,oday it is the custom to measure shape by methods which are indirect. Landmarks (reliably locatable points, “places to point to”) may be located in images, and distances and angles among them measured as variables ; or a center of coordinates may be set and Cartesian or polar coordinates for the curve extracted and condensed into moments or Fourier coefficients. In these procedures the underlying physical object, which I presume to be a biological outline, is never actually given a
FlTTIx(;
(‘oNI(:
sI’(TIoSS
‘l-0
sc.\‘rTI~;l:l~;l~
1 ).\‘I‘.\
. l
.
l
.
.
. .
.
. . . l
. .
. .
.
. .
Frc. -4. Typical :m-length does
l
.
.
data for ol7tlilrc-fitt,ill~: rmt mist.
.
frulc:tic)n:tl
rc~l:&mships
of Ilnkrl~rwn
strllcture.
I
geometric expression free of an impowd coordinate system. Then all these methods study a measurement vector, not the outline. h critical review of the scattered literature on this subject may be fourld in Bookstein [la, Chaps. II, III]. There is an alternative to be had from the general theory of plane curves. The curvature of a smooth plant curve is the rate of change of tangent direction with respect to arc-length, the reciprocal of the radius of the “osculating circle,” t’he circle making closest contact with the curve at a point. From curvature as an cbxplicit function of arc-length, measured from any origin on the curve, the whole curve can be reconstructed without error. Then a curve is ~-holly rcpresentcld 1))
FRED
L. BOOKSTEIN
FIG. 6. Two-arc conic spline to digitization of the outline of a pelecypod, Pisidium compressum. Data from Burch [17, p. 141. The open area at the top is a projection of the hinge obscuring the line of closure.
its curvature as a function of arc-length, the intrinsic representation (Stoker [13, Chap. II]). It does not matter at which point of bhe Cartesian plane the integration begins or in which direction the tangent line initially points; vary these and the intrinsic function draws out curves all perfectly congruent. Unfortunately, such a function is never given to us in the course of empirical research. We have, rather, selected points upon such a curve, in the manner of Fig. 4. Before we can proceed at all, these isolated points must be combined into a smooth curve somehow. It is clear that we cannot connect the given points with straight lines, for the resulting function has it’s tangent angle undefined at exactly the points, now vertices of a polygon, at which we need to sample it. In case the data were chain-coded, the situation is even worse, for every azimuth is one of eight values. Nor can we say that the data approximate to a “true” curve, that with greater and greater density of points sampled we would converge to a limiting “true” locus. At the limit, the tangent angle is not a continuous function of arc-length, for arc-length does not exist. The greater the detail in which one examines any real locus, natural or fabricated, the more irregular it appears, the longer the arc-length, and the denser its local extrema of curvature. See Steinhaus [14] and Mandclbrot [15] for a discussion of this remarkably general phenomenon, that arc-length generally increases exponentially with the scale of observation. This paradox is not peculiar to biometrics ; it complicates also the study of seacoasts, galaxies, turbulence, evolutionary histories, and soap. Then to utilize the intrinsic representation theorem in data analysis it is necessary to construct a rectifiable curve fitting the data before we proceed. WC could, for instance, pass a suitably smooth curve through every point of the digitized outline, thus enshrining all measurement error on a par with form itself, Such curves would tend to be a good deal longer than what WC intuitively feel is
i
70
FRED
L.
BOOKSTEIN
the appropriate “length” of the outline. We could contrariwise try to capture the form by a “best-fitting” specimen of some family of simple mathematical curves. These families, however, would have to have a number of parameters which is some multiple of the count of “pieces” of outline (arcs between landmarks) to be satisfactorily fit. They are most simply generated as splines such as WC have been discussing: linked series of curves each of which best fits the data in its vicinity and passes smoothly into its neighbors. Then the value of conic splining is that it supplies a substitute curve, easily fitted closely to the data, with a curvature computable simply from the coefficients of the arcs and with arc-length equal to a familiar elliptic integral offered in the various collections of scientific subroutines. Regular local features of shapeindentat,ions inside the spline or bumps outside-become new shapes which can be analyzed further. The spline captures the “general shape,” and provides an arc-length argument for the local features, in a wholly coordinate-free way. Measures of this general shape can then be derived from inspection of the splined intrinsic function and its variation over sets of shapes, without any ratios, distances, or other “shape variables” having been specified in advance. I forsee t1y-o statistical methods by which these functions might be analyzed. In one, a multivariate measurement vector may be constructed by sampling the intrinsic function directly at finitely many point’s, corresponding to either actual landmark locations or aliquots of arc-length. Then the full information about the curving outline is preserved until the last stages of the analysis, whereas it is lost early if one samples isolated points from the outline. I suggest a specific statistical algorithm in Bookstein [12, Chap. IV]. In the other, the spline coefficient vector V is used itself as the measurement vector. In an Z-arc spline with one “corner,” V exists in projective (21)-space. Means, group contrasts, and the like can be examined there by multivariate analyses upon the cross-product matrix of direction cosines of V. These techniques should apply whenever the derivative of curvature has few changes of sign over the outline of a form. Figures 5, 6, and 7 are examples of fits to suitable biological data. The first of these is a sea-cucumber outline that I digitized from a 3- by j-cm image in an old zoology text. It can be seen how closely a three-arc conic can be made to follow this creature of very odd shape, a shape, in fact, describable only by reference to the smooth spline. The next figure is of a clam shell. The fit of two elliptical arcs to this natural form is close enough that the intrinsic function should differentiate among growth stages of this organism and among congeneric adults, by shell outline alone (should anyone care to digitize more of the images). Figure 7 is of more general interest. Of particular irksomeness in craniometrics is the absence of identifiable landmarks over the vault of the human skull. It is very difficult to describe that form having only quantities (distances, angles, ratios) of rather large reach. Figure 7 shows that an arbitrary specimen of this form is quite closely fit by either two or three conic arcs splincd together-tht: optional third arc makes it possible to represent the little bump above the inion, at the far left. There is then no need for landmarks at all in this part of the skull, but the whole shape itself, as smoothly reconstructed by the spline, may be used
FITTING
CONIC’ SECTIONS
'1'0 S(‘ATTI~;Kl:l,
I)hTA
71
;ts a function-valued or coefficient-valued measure in analysis of gro\vth pat,tcwlF. racial differences, evolution of t’he hominid form, and t.he l&c.
1. K. A. Paton, Conic sections in chr~~mosome analysis, Puttern. ftecuy. 2, 1970, 30. 2. Ii. li. Pat,on, Conic sections in automatic chromosome analysis, Machine Znlelliyrrlce 5, 1970. $11. :;. II. II. Biggerstaff, Three variations iu dental arch form c4nmtcd by :I quadratic eyuatioil, J. Dent. Res. 51, 1972, 1599. 4. A. Albano, Representation of digitized contours in terms of conic arcs and straight-lirrc segments, Computer Graphics and Image Processing 3, 1974, 23. 5. I). I(. Cooper and N. Yalabik, On the computational cost of approximating and recognizing n,lise-perturbed straight lines and quadratic :u-cs in the plane, Il?EE Trans. Comprrters c-25, 1976, 1020. 6. 1:. ( ;nanadesikan, .Ilethotls for Stutisticd DC&I Analysis of J/ ultizuriute Observations, Wiley, New York, 1977. 7. K. Pearson, On lines and planes of closest tit to systems of point,s in space, I’hilos. dlag. Ser. /i 2, 1901, 559. S. C. Ii. Rao, Lineur Statistical Inference and Its Applications, 2nd ed., Wiley, New York, 197:;. 9. I). J. Poiricr, Piecewise regression using cubic splines, J. Amer. Slatisl. Assoc. 68, 1973, 515. 10. S. \\‘old, Bpline functions in data analysis, Technomet. 16, 1974, 1. 11. I). 1‘. Thomas, Pseudospline interpolation for space curves, dfath. Conrp. 30, 1976, 58. 12. F. I,. Bookstein, The Jleasurement of Biological Shape and Shape Change, Springer-Veriag. NPW York, 1978. 18. J. J. Stoker, Di’erential Geometry, Wiley, New York, 19ti9. 14. H. Steinhaus, Length, shape, and area, Colloq. Math. 3, 1951, 1. I ,5. B. AIandclbrot, Fractals: Form, Chance, Dimension, Freeman, San Francisco, 1977. 16. li. Storer, General Zoology, McGralv-Hill, New York, 1951. 17. J. Burch, Freshwafer Sphaeriacean Clnms of ,Yorth Americu, ;\I:dacological Publications, II:imbirrg, Mich., 197,5.
OIt.Ll?‘HICS.\ND IKLGIC PItOCIGSSIKG
Fitting
Conic
9, 56-71 (1979)
Sections
to Scattered
Data
FRED L. BOOKSTEIN Center for Human Growth and Biostatistics, The University
Development of Michigan,
and Departments of Statistics Ann Arbor, Michigan @109
and
Received January 6, 1978; revised March 28, 1978 The problem of fitting conic sections to scattered data has arisen in several applied literatures. The quadratic form Ax2 + Bxy + Cy2 + Dx + Eu + F that is minimized in mean-square is proportional to the ratio of two squared distances along rays through the center of a conic. Considerations of invariance under translation, rotation, and scaling of the data configuration lead to a straightforward method of estimation somewhat different from earlier suggestions. The method permits an extension to conic splines around extended digitized curves, expediting a smooth reconstruction of their curvature. Some examples are presented indicating how the technique might be applied in morphometrics. 1. INTILOI~UCTION
In the course of investigations into the shape of plane curves, several authors, attempting to fit conic se&ions to scatters, have recently and independently discovered algorithms that arc reducible to conventional linear computations. All minimize the sum-of-squares of a form Q(x, y) = AJ~ + Rxy + Cy2 + Dx + Ey + F over a scatter in the (z, y)-plane. Paton [l-2] sets A2 + B2 + C2 + D2 + E2 + F2 = 1 as the constraint on this minimization, whereupon the solution vector (A, B, C, D, E, F) is a principal component of the cross-product matrix of the variables x2, zy, y2, 2, y, 1. Biggerstaff [S], Albano [4], and Cooper and Yalabik [5] set F = 1 as the constraint, whereupon there arises a set of normal equations for A, B, C, D, E in terms of the same cross-product matrix. Gnanadesikan [S] set’s A2 + B2 + C2 + D2 + E2 = 1, yielding for extrcmum the last principal component of the variance-covariancc matrix of x2, my, y2, 2, y. All these authors loosely characterize the optimal conic as “best-fitting” of all the tonics of the plane, but the nature of the fit each is optimizing, whose computation is so surprisingly straightforward, has not been elucidated. WC seek an estimation rule which is general, simple to compute, and invariant. Generality requires a single algorithm by which we fit both ellipses and hyperbolas, with circles and parabolas included as constrained special cases. For instance, we cannot generalize the fitting of parabolas usual in quadratic regression, in which error-of-fit, is measured parallel to the axis of the parabola, because the
0146~GG4X/79/010056-16$02.00/O Copyright All rights
0 1970 by Academic Press, Inc. of reproduction in any form reserred.
Iti\~ariatwc is to lx \vitli rcq)c~c~t~ to tllc> c~cjuiform h~~roup of transi’ormatioits oi’ t Ii(b l~~u~litloati platir : rotations, t~ratislatiotis, and c~liattgc~sof scalp. if tlio wortlito an ~~l~~m~wloi’ tiutctr of t’lici data aw trattsformc~d, point 1)~ point, wwrding this group, tlwti tlw kst-fitting coriic~ to the ti~\v wnttc>r must 1)~sc~~at*tIy t tics ws~ill of t.lic same transformation applicstl to tlica cwtiich I\-lticli ~2s hwt-fitting to tlirl original scatter. ‘l?i~~wtiniatiott rules proposc~d lwr(~ mitiimizw a cluadratic* form in tlw utiktio\vn cotticd l)aranic$t~ri sut)jwt to ati w:wt wrlstraiiit on tlw val\tct of anothc~r cluadratic form. It \\ill 1~~sIw\vtt t 1~~1Iwth tlicw forms arcs invariants und(ar the action of th(a c~c~ttit’ornlgrortl), and lwtrw so is tlw tic4 wtimstion proccdurc~. To insist on this invariant is to alickr our intttitivc notion of curvesfitt,ittg “to a dcpondwit varisl,l(s.” \\.ltcsti al~wiss:~ and ordinatc~ wpwwtit 2111 arbit,rary coordinate q~stc~m, or arbitrary oriotttatiott. plapcbd mpor~ an undwlyitt~; gc~omc~trical ol)jxt, swh us tlicl outliirc~ 01’ a skull, t It(x axcbs wparatrbl!. liavo no gc~omc%ric mclaning at all; tlw~. arc’ \\ Itoll~~ (‘Oliilil(‘tisitr;Lt(’ ant1 may lw rotattd frwl~~. WC mrtst gclt)tlw sam( j gwnic~tric~ rcxsttlt iri all (2~s. 2. (:I~:OMI~;Tl:I(’
IN’~I~:I:l’l~l~:T.21’IO?j
01; Q(.r, !/)
1~~1a giwn wnt,ral conic wctioti 11k~vc~ c~cluatiott c)b, y) = .l~:’ + lj.r,y + C’,$ + 11.r + 1Q + E’ = 0. I,(+ (xl,, yO) 1)~~an arbitrary Jwint. ot’ tltct plaw, and s(‘~ (xl, ~1) qua1 to the point on 11~~1 conic on the ray from thtl wntw through (.I’,,, !I,~). Cltoosc: a IWW coordinat.cl systclm, d(wotcld 1)~ primw, n-hic*h is wntercd at t ll(t cotitc>r of tlw conic and aligrictd with its priwipal aws. J>~finc~tlicl form Q’(.r’, !/‘:I on tSlrcbrw\\- coordinates so that Q’(/ , !/‘j = ()(.I*, y) OII tlic l)lanc (A, y’). In tlti!; IWK coordinate systrm, Q’(.r’ y’) = :l’.P + 13’jg’2+ 1.” for somchA’, B’, F’. In the primc\d systcsm! the> wntc~r is thcb origin. (,(I’, 0’). By c~olliwarit\~ OI (,r,,, g,,), (x1, ,ql), and tho ccwtor, (bxpr(wed in the, prim(bd cwordinat.c> s;!stc~m. UC’ Ilaw d.‘~‘!/‘I = do ‘~‘0 or .r’,) = .r’] t!/‘,, !/‘I). Th(~ti
-
/“‘[rJ/‘,, ‘!/‘,)”
I j,
i II
since (c’~, g’lj is on the conic C)‘(J’ , y’) = 0. But the term (g’Cv’g’l) is just tllcbratio of distances writer-to-point and c,c:itter-to-c~otti(~don p the ray JVCarc using -tlw
5s
FlU3:D
L.
BOOKSTEIN
Ellipse
Parobofa FIG. 1. The geometry of error-of-fit conic. The error is always proportional
Hyperbola around a conic section, to [(d + d$/#] -
as measured 1.
by the
equation
of the
ratio of d + dI to cl in Fig. 1; for projection onto any axis preserves distance along a line. Hence &(x0, 1~01= [I (d + dJ2/C121 - 1, and as the equiform transformations do not alter distance ratios they preserve this proportionality (though they may alter F’). We have [(cl + d$ - cP]/cP = &(cZI + 2d)/d2. When the data are quite close to a conic, cl1 is small, dl $ 2d is approximately equal to 2d, and Q(z,,, y,,) is approximately proportional to dl/d, the distance from the conic along a line through the conic’s center, measured in units of distance from the center to the conic in that direction. The fitting of a parabola is a limiting case, exactly transitional between ellipse and hyperbola. As the center of an ellipse moves off toward infinity while its major axis and the curvature of one end are held constant, the rays from its center become parallel, and the difference of squared distances which is the numerator of Q becomes linear in distance. The quantity being minimized in the fitting then approaches the usual deviation, Euclidean distance from the curve measured parallel to the axis.
For a circle, I propose using the instance of least mean-squared Q. Let there be given points (~1, yl), . . ., (xn, yn) and let the sample linear regression of x2 + y2 on x, y, 1 be x2 + y2 - 2ax + 2by + c, where a, b, c are the least-squares
4.
IN\‘.\I:L.\N’I’
E‘OI:
‘l’lll~:
(;I,:NI~:I:.\I,
C’ozIt
‘
1~01~an ahitrary conic, tlic gxwc~ralizatiott of Ilit> cwor-of-fit i.r - 0’):: + (y - b)’ - A? is clwrl\tlitr valnc~ (J(,.i’, y) = .I tl hy miriinnizing CGrLB,., (i/d. + I)!/ + (.jz sul)j~.~~tto tlw cwnstraitll (1” + /I” = 1. For conics, w: might tq’, IJ~ clirwt niimict~y, the, for~rl -1’ + 1:’ + C”. I’IIfortutIatc>ly, the liypcrbola Lh.y = 1, with .t” + 13’ + f’” = 1, can lx, rotated into th form .J.‘)- p2 = 1, with A’ + 11” + C2 = 2, and so tit,, first would automal imlly fit t\viw as l~adly if tltc data w(‘rc rotat cd 15”; hut. tltcs sclccted normalization tnust be invariant, utidw rotation. Ititroduc~ti~~w lo tlw algel~raic tliwry of conic85 tell us that the forms A + (’ ant1 11’ --- .4.1C aw invariant iitickr the I3tclideati Sgrou~),l’hc~ onl\- r)ositivc~-tl(~finit(~itiv:tri;tlit t Ilat (%ti l~c’ formcbtl from tlrwc cpin~~t i-
tics is (fl + 6)” + (B2 - UC)/2 = A” + B2/2 + C2. We may set a unique scale for Q, and thus proceed to actual minimization of ZQ2, b*y setting this norm equal to some conventional constant value. I suggrst the value 2, so that the equation of a circle is in the usual form. Given data points (51, yl), . . . , (z~, yn) scattcwd about a conic, I propose to estimate the conic of best fit by minimizing C;[G,(zi, yJ12 subject to the constraint VDV’ = 2, where V = (A, B, C, D, E, F) is the vector of coefficients to be estimat.ed, and D is the matrix diag (1, 4, 1, 0, 0, 0). The minimand may be written VSV’, where S is the scatter matrix about 0 of the 6-vector (z?, xiyi, yi2, z;, y;, 1) as a function of the data points (x;, vi). The invariance of this procedure is proved straightforwardly. If the coordinates of the data points are rotated or translated, the constraint A2 + B2/2 + C2 is unchanged, likewise the term F’ of Eq. (I), which is gotten by undoing the transformation; and the distance ratio, the other multiplicand of (l), is clearly unchanged. Then the fitting is invariant under translation and rotation. If the coordinates are resealed by k, (2, y) + (x’, y’) = (kx, Icy), let V, be the vector of coefficients of the resulting best-fitting conic by this algorithm. The conic As2 + . . . + F = 0, after resealing and renormalizing, becomes Ad2 + Bz’y’ + Cyf2 + kDx’ + kEy’ + k2F = 0. The resealing has replaced the matrix S by SI, = DkSDb with Dk = diag (k2, k2, k2, Ic, lc, 1). Note that DkDDk = k4D. Extremization of Vk(DkSDk)Vk’ for V,DVk constant is equivalent t’o extremizing it for VkDkDDkVk’ constant. Since Dk is nonsingular, the extremum is V/Dk, where V is the extremum before resealing. But V/Dk gives the coefficients of the old conic in the new system, for k2V/D~ = (A, B, C, kD, kE, k2F). Hence this method of fitting is invariant under equiform transformations. As the authors cited in Section 1 all use D and E in the norms for their fitting, none of their solutions are invariant under the Euclidean group. Those which involve F, further, are not even invariant under translation : For by setting F = 1 it becomes impossible to fit tonics through the origin (F = 0) at all. These authors must therefore forcibly standardize their data before conic-fitt’ing in any of various ad hoc ways. Translation to center-of-mass and rotation to principal axes of scatter are the most common. This approach seems unsatisfactory. Dat’a must be fit in a manner invariant under t’he Euclidean group, not arbit’rarily constrained with respect to it. 5. ESTIMATION
We wish to minimize
(VI 1V,), both components
FOR
THE
GENERAL
CONIC
VSV’ subject to VDV’
= constant. of length 3, and let S bc partitioned
Partition V into correspondingly :
SC(+ I )* x12
21
VSV’ = VlSllV’l
+
x22
2v1s,,v’2
+
v2s22v2.
We must minimize this subject to VIDIV’ 1 = constant, where DI = diag (1, 8, l)*
‘I‘llis
i9mputution
upon
1 tw
cdimat
ing,
I\-hiw
M
anil \vi’
=
matrix niinus (‘.gL, ktraint Sctt \\:lys,
a wt
is a matrix,
cdrcmllm 1-D 1”
remains
parami+~r5
scc~k:
-
M
iliwotiss it, must h may ing s~)niil
of
iwnstraint Hao
l’hc
vcdor
I’M
=
0
rlwil [S,
tlli> -I)
ortlinaty-
\vib plains
tlc~noti~
tlw
lc.(i]
instruct,5
smallid
iilvc~r,q~. invcw~.
\vitll ‘l’lic~
of riqwt
to
mat,ris D
in IIS
T-S\”
\\ Irtww
linear of
togi+twr
i~igi~nvcdor D
arbitrary
Ci-wctor
\\-ritt,csrl
Swtion T- of
is
gi~ni~raliaid an
\\.lic~n to
,q is
(M’S--‘M)-M’S any
haw
I- is
t81wn
wwdant, (I
t ractal)l~~
If
the in this
T’M
formula
subjcd
to
c~igi~nvallw
S.
Hew
wit
form
largwt, S mud
ma!.
constraints
i~oiBkii~nts
the
for
tlw
of
tllc,
supcwcript, clcfiniti~.
and
thi>
CY)II-
1x3 rdunda~~t.
linear ilriitit
const~raints usdul.
crodw IIc~w
:iw
this wn1i’
avuilal,li~
diyycw
I,iwsil)ilitiw
:
of
friwlom
0, thi.
constraints
1~’ positive
1)~s singular
arc’ =
in
various
. FIG.
2.
Two tonics abutting
at a. knot where the tangent crosses the curve.
(5) 2Aso + By, + D = 0, Bxo + 2Cy0 + E = 0. The conic has its center at (x0, ~0). (6) Ar? + Br~yl + Cy,Z + D.z:~+ Eyl + P = 0. The conic goes through (Xl,
Yl).
(7) @Ax:1 + BYI + 0) cos cr + (B.xI + 2Cy1 + &) sin o( = 0. If the point (x1, yl) is on the conic, this is the criterion that the tangent there make angle LY with the positive x axis. (8) A (- 21 sin LY+ R cos2a) + B (-05y1 sin2 a! + 0.521 cos2 o( + R cos a sin a) + C(y, cos c11+ R sin2 a) + D( -0.5 sin a) + E(0.5 cos a) = 0. The conic through (x1, yl) with tangent making an angle cr with the positive x axis is by this constraint required to have radius of curvature there equal to R. 7. CONIC
SPLINING
There is a rapidly growing literature (cf. Poirier [9], Wold [lo]) on polynomial spline regression. In this technique, the expected value of a dependent variable is fitted by different functions-usually polynomials-over distinct ranges of an independent variable. At ‘(knots,” where these ranges abut, the predictions are required to agree in value and in their first few derivatives. After the knots are set, all the polynomials are computed at once, by minimizing the tot’al squared error-of-fit about all of them in one clever dummy-variable regression. The method becomes much more complex when error-of-fit is not always measured along a fixed axis. The constraints which force the curves to line up become nonlinear in the coefficients being estimated, and linear methods fail (cf. Thomas [ll] and references therein). In particular, one cannot spline a closed curve linearly by polynomials. But just as parabolas (and higher-order polynomial graphs) generalize straight-line regression, so do tonics (and higherorder algebraic curves) generalize st’raight-line curve-fitting in general orientation. It-is then of interest to explore the possibility of fitting multiple tonics to separate parts of a curve simultaneously in a single linear computation generalizing the previous discussion for fitting tonics singly. , Suppose we have a data set in the plane which looks like two tonics, in the manner of Fig. 2, impinging upon each other at some visible nexus P. It would be nice to fit two conic segments simultaneously, one to the left of P, one to the right, which minimize the net error-of-fit. We would like the tonics to connect at P without any corner, that is, to have the same tangent there. The inclination
ITlwN(:
(10X1(’
S1~:~‘1’10NS ‘1‘0 s(‘.~‘1”1’I~:l: I,:I) I ).\‘I’.1
(it:
is not givcw ill IIu: flata, Iio\vc~vc~r; \\‘I: n11bd f~sliumtc it, I oo, as that position of 1lit: langc~nt lint> at I’ for whkh t ht\ bc~sttwo conic?; through I’ \vith that tangcbnt, 011~: I’or f~af*lk sid(x, togcstIiw huvc~ lfmt nlcan-squarfd fvorof-fit. To fbxf>cutc this plan \vf’ must link thfl tlormalizations of thf: tmm conic+. The bwt way to do this is by insisting that, a givctn fxrror-of-fit near I’ 1~ tlio Sam0 I\-hichww caonic WCconsidw it on. To a linclar approximation, thcl error for a point such a:: R rclatiw to a conic (r)(s, IJ) = 0 through 1’ is c~lual to t,hc>gradicant of the walsr fif4d CI)(.r!y) timw thcl distant from I’ to 22 nornml to tlifb fs~irvf’. ‘I‘hcw~ \vill bf> (qua1 for t,licBt\\-0 f~oriks impinging at I’ if’ and only- if the, gratlicnts of t,hc> t ~1-0 scalar fields arc ftflual at I’. But to spccif~. that thwl two gradifsnts l)fb c~111alis jwt a [lair of’ linfbar COIIst raints! for thf: gradients are liwar forms in th(l cwffic4c~nts of the cork. Spc’ificdl\- : if thf> two conies arc’ (Jzi.r. y) = -1 ,.? + B,J!/ + C1!/ll + 11,~ + E,y + P = 0. cI’. s): thcxu thcl t ii-0 con&tc~nc~!~crit(~ria :WI I’ = 1, 2: and if P have t hc coordinate
0I t,liis mul)ual tangfwi,
“:l 11’+ up + D, -- (2.1 .‘t’ + B,s + II,) = 0, 3c,s + n,r + I:‘, -
iX’,s + II,,, + 7:‘,) = 0.
IVe must, also insist, that both caonics actually pa.hs through thr point (I,: s) at ivhich KC:arc specifying t,hAr tangents, and this is t\vo mow constraints :
..w + R,,,~+ c,2 + rjzr+ E:‘,.~ + I?,= 0,
I’ = I,?.
\\Tc can estimate both conks, all 12 cof+Ikicnts, simultanc~ously by minimizing their total squared error-of-fit subject to t hfi normalization -4,’ -+ HIZi’2 + C1’ = 2 and th(h four constraints above. This tact’ic is easily gcncralized for thft fitting of conic spliws to cxxtcndcd chain?, cwn closed curves. Let thaw be I arcs, and on the it,h arc data points (I,,, vi,). ,j = 1. . . ., f/,, 7, = 1, . . ., I. Let the arcs abut at knots (.fi: !I,):, i = 1, . . . , I, and sot (II. 1, yf+l) = (./.I, yl). Let) Si he thf: G X Ci scatter matrix about 0 for tlift points (Xi,‘, sX,j,!Jij,uij 2, T,j, !/aj, 1) dflrivctl from the data points (s,~, ?/lj) of thfn it,h arc. Define S, (61) X (GI), equal to blockdiag (A’,, . . . , SI). Consider a conic* spline which is equal to AL.? + R,.fy + C,y’ + D,z + E,y + I;‘i on the ith arc, i=l I, and let. Al, 1 = .-11, . . . . PI .m1= PI. The error sum-of-squarf~s al)oilt the itll’ki for the data of the I’th arc is (.A;, 13,, C,, I),, I:‘,, F,)AS’,(-I,. Z3,, C’,. II,. I!‘;, Pi)‘. Dcfinc 1’ = (AI. . . . , k’lJ A,. . . . , I”,-1, rlr, . . . , F,), a ((il)-vwtor. l’hf~~ the tot$al sum-of-squares for wror-of-fit about all tlicx arcs is 1 JI-‘. This may I)(, minimized subject to a set of homogcncous linear caowtraint:: 1-M = 0 by the formula of Rao quoted above. In t5his application, S must lw nonsingular ; thcsrl each S, must, be nonsingular, which implies that 71, > 5, all i, and that thr point:: of all arcs must lie in “general position” (not all upon on< ronip).
G4
FlmIl
L.
BOOKSTEIN
and also subject to a normalization. As the tangency constraints explicitly link normalizations of one conic t,o the next, we can only assign one normalizing constraint for any subarc chain of continuously turning tangent. In the usual case where this subarc is the whole spline, as in Figs. 3, 5, 6, and 7, this means that any constraint Ai + Bi2/2 + Ci2 = constant d&ermines the normalization of all the arcs all the way around. 8. WORKED
EXAMPLE
The vaguely biological outline, Fig. 3a, derives from a freehand sketch on graph paper. I chose three “landmarks” : No. 1, the “corner” in lower center; No. 2, the extreme right point; and T\‘o. 3, the ext’reme top point. The word “landmark” is used here, in preference to the “knot” of the splining literature, t’o suggest determination by anatomical criteria rather than asympt.otic optimization tactics. If this were a skull X ray, Ko. 1 might be the audit,ory mcafus, No. 2 the bridge of the nose, and No. 3 the bregma. Points between the landmarks have been digitized crudely by cyc. The spline arcs are computed by a program executing the following algorithm. (a) The data are read in arc by arc and the scatter matrix S accumulated and stored. A vector of 18 unknown coefficienOs V = (Al, . . . ) F1! . . . ) Aa, . . , F3) is allocated, bearing the conic coefficients arc-wise. (b) The constraints are written out in full. Arcs 3 and 1 pass through point 1 of coordinates (0, -3). This requires
Ai.
+ Bi.0 + Ci.9 + Di.0 + Ei.(-3)
+ F, = 0,
i = 3, 1,
or
Similarly,
9Cs - 3E3+FJ = 0,
(1)
9C1- 3E1+ FI = 0.
(2)
arcs 1 and 2 pass through point 2 = (8.5, 0), requiring
and for point 3 = (-2,
At landmark are no further constraint and first conic and
72.2581+ 8.50, + F1 = 0,
(3)
72.25A, + 8.50, + Fz = 0;
(4)
5.7),
4Az - 11.4Bz+ 32.49cz - 202 + 5.732 + Fz = 0,
(5)
4A3 - 11.4B3 + 32.49Cs - 203+5.7Es+Fs
(‘3)
= 0.
1 we allow the arcs to int’ersect at an arbitrary angle, and there constraints to be had here. At landmark 2, we pick up a seventh an eighth from componentwise equality of the gradients of the the second: 17A1+
Dl-
17Az - Dz = 0,
8.5B1+ El - 8.5Bz - Ez = 0.
(7) (9
FITTTKC;
.
(:ONI(’
.
SWTlON~
TO ,S(‘.4TT1~:l:I~:I)
l).~T.\
. .
. . . .
a
FIG. 3. (a) Hypothetical data set. (b) 4 three-arc ronic spline approxim:ltinp knot,s at the points numbered in (a).
The final two constraints
-&I,
this scAtw, with
derive from point X similarI).:
+ 5.7Ba + D:! + 1A:j - 5.7B3 - 1)s = O>
11.4C2 - 2Bz + E’z - 11.4C3 + 2B3 - E, = 0.
(!I) (10)
(c) For any vector 17 of coefficients for any spline, the total error-of-fit is the scalar VSV’, where S is the cross-product matrix accumulated in (a). WP
GG
FRED
L.
BOOKSTEIN
nocd to niinimizo this for fixed invariant norm 111~+ Bi’/2 + Cl2 = VD V’, D = diag(l,$,l,O, . ..) 0), subject to constraints (1) through (10). Conversely, we can instead maximize VDV’ for fixed error VSV’ subject to the same 10 constraints. The theorem of Rao quoted previously applied verbatim here to provide a computing formula for the solution wc seek. A single eigenvector of a certain large matrix manifests the coefficients of all three tonics, in blocks of six, for its loadings. The spline corresponding to this eigenvector exactly passes through all the knot,s with well-defined tangents there and minimizes the total error-of-fit from the tonics in between. This curve is conic between every pair of consecutive landmarks, and its tangent turns continuously throughout, except where we have expressly omitted t,he constraint that it do so. (d) These arcs are drawn out between appropriate endpoints (by short segmental approximations) and then superimposed over the data in Fig. 3b. The spline fits the data quite well except on the far left, where the conic arc has only one curvature maximum whereas the data apparently have two (the northwest and southwest “corners”). This systematic deviation can, if we wish, be declared a “local feature” of the outline relative to the spline, as the nose is a local feature relative to the plane of the face ; it is not necessarily a flaw in the fitting. The splining-with-landmarks replaces the observed (digitized) boundary with smooth conventional arcs. Its tangent and normal, turning smoothly, provide a moving coordinate system to relate deviations from the “standard” in different regions all around the curve. The spline is more elegant and more easily read than a detailed tracing, as local features do not disrupt the global aspects of the shape, the ebb and flow of curvature around the form. Should a landmark, say the second, bc itself known only with error, one can free the tonics from passing exact’ly through it. The revised equations would require, instead of Eqs. (3), (4), (7), and (8), only that the landmark have a common distance from the two tonics and the same polar with respect to them. The appropriate constraint equation replaces (3) and (4) with the equality of their left-hand sides : 72.2581 + 8.501 + F1 - 72.2582 - 8.502 - Fz = 0;
Eqs. (7) and (8) persist unchanged. With this revision, the landmarks serve only a housekeeping function as knots, and do not serve any special role as data themselves. 9. RECTIFIABLE
CURVES
AND
MORPHOMETRIC
DATA
ANALYSIS
The technique just presented applies directly to a nagging problem in morphometrics: the function representation of outline. In morphometrics t,oday it is the custom to measure shape by methods which are indirect. Landmarks (reliably locatable points, “places to point to”) may be located in images, and distances and angles among them measured as variables ; or a center of coordinates may be set and Cartesian or polar coordinates for the curve extracted and condensed into moments or Fourier coefficients. In these procedures the underlying physical object, which I presume to be a biological outline, is never actually given a
FlTTIx(;
(‘oNI(:
sI’(TIoSS
‘l-0
sc.\‘rTI~;l:l~;l~
1 ).\‘I‘.\
. l
.
l
.
.
. .
.
. . . l
. .
. .
.
. .
Frc. -4. Typical :m-length does
l
.
.
data for ol7tlilrc-fitt,ill~: rmt mist.
.
frulc:tic)n:tl
rc~l:&mships
of Ilnkrl~rwn
strllcture.
I
geometric expression free of an impowd coordinate system. Then all these methods study a measurement vector, not the outline. h critical review of the scattered literature on this subject may be fourld in Bookstein [la, Chaps. II, III]. There is an alternative to be had from the general theory of plane curves. The curvature of a smooth plant curve is the rate of change of tangent direction with respect to arc-length, the reciprocal of the radius of the “osculating circle,” t’he circle making closest contact with the curve at a point. From curvature as an cbxplicit function of arc-length, measured from any origin on the curve, the whole curve can be reconstructed without error. Then a curve is ~-holly rcpresentcld 1))
FRED
L. BOOKSTEIN
FIG. 6. Two-arc conic spline to digitization of the outline of a pelecypod, Pisidium compressum. Data from Burch [17, p. 141. The open area at the top is a projection of the hinge obscuring the line of closure.
its curvature as a function of arc-length, the intrinsic representation (Stoker [13, Chap. II]). It does not matter at which point of bhe Cartesian plane the integration begins or in which direction the tangent line initially points; vary these and the intrinsic function draws out curves all perfectly congruent. Unfortunately, such a function is never given to us in the course of empirical research. We have, rather, selected points upon such a curve, in the manner of Fig. 4. Before we can proceed at all, these isolated points must be combined into a smooth curve somehow. It is clear that we cannot connect the given points with straight lines, for the resulting function has it’s tangent angle undefined at exactly the points, now vertices of a polygon, at which we need to sample it. In case the data were chain-coded, the situation is even worse, for every azimuth is one of eight values. Nor can we say that the data approximate to a “true” curve, that with greater and greater density of points sampled we would converge to a limiting “true” locus. At the limit, the tangent angle is not a continuous function of arc-length, for arc-length does not exist. The greater the detail in which one examines any real locus, natural or fabricated, the more irregular it appears, the longer the arc-length, and the denser its local extrema of curvature. See Steinhaus [14] and Mandclbrot [15] for a discussion of this remarkably general phenomenon, that arc-length generally increases exponentially with the scale of observation. This paradox is not peculiar to biometrics ; it complicates also the study of seacoasts, galaxies, turbulence, evolutionary histories, and soap. Then to utilize the intrinsic representation theorem in data analysis it is necessary to construct a rectifiable curve fitting the data before we proceed. WC could, for instance, pass a suitably smooth curve through every point of the digitized outline, thus enshrining all measurement error on a par with form itself, Such curves would tend to be a good deal longer than what WC intuitively feel is
i
70
FRED
L.
BOOKSTEIN
the appropriate “length” of the outline. We could contrariwise try to capture the form by a “best-fitting” specimen of some family of simple mathematical curves. These families, however, would have to have a number of parameters which is some multiple of the count of “pieces” of outline (arcs between landmarks) to be satisfactorily fit. They are most simply generated as splines such as WC have been discussing: linked series of curves each of which best fits the data in its vicinity and passes smoothly into its neighbors. Then the value of conic splining is that it supplies a substitute curve, easily fitted closely to the data, with a curvature computable simply from the coefficients of the arcs and with arc-length equal to a familiar elliptic integral offered in the various collections of scientific subroutines. Regular local features of shapeindentat,ions inside the spline or bumps outside-become new shapes which can be analyzed further. The spline captures the “general shape,” and provides an arc-length argument for the local features, in a wholly coordinate-free way. Measures of this general shape can then be derived from inspection of the splined intrinsic function and its variation over sets of shapes, without any ratios, distances, or other “shape variables” having been specified in advance. I forsee t1y-o statistical methods by which these functions might be analyzed. In one, a multivariate measurement vector may be constructed by sampling the intrinsic function directly at finitely many point’s, corresponding to either actual landmark locations or aliquots of arc-length. Then the full information about the curving outline is preserved until the last stages of the analysis, whereas it is lost early if one samples isolated points from the outline. I suggest a specific statistical algorithm in Bookstein [12, Chap. IV]. In the other, the spline coefficient vector V is used itself as the measurement vector. In an Z-arc spline with one “corner,” V exists in projective (21)-space. Means, group contrasts, and the like can be examined there by multivariate analyses upon the cross-product matrix of direction cosines of V. These techniques should apply whenever the derivative of curvature has few changes of sign over the outline of a form. Figures 5, 6, and 7 are examples of fits to suitable biological data. The first of these is a sea-cucumber outline that I digitized from a 3- by j-cm image in an old zoology text. It can be seen how closely a three-arc conic can be made to follow this creature of very odd shape, a shape, in fact, describable only by reference to the smooth spline. The next figure is of a clam shell. The fit of two elliptical arcs to this natural form is close enough that the intrinsic function should differentiate among growth stages of this organism and among congeneric adults, by shell outline alone (should anyone care to digitize more of the images). Figure 7 is of more general interest. Of particular irksomeness in craniometrics is the absence of identifiable landmarks over the vault of the human skull. It is very difficult to describe that form having only quantities (distances, angles, ratios) of rather large reach. Figure 7 shows that an arbitrary specimen of this form is quite closely fit by either two or three conic arcs splincd together-tht: optional third arc makes it possible to represent the little bump above the inion, at the far left. There is then no need for landmarks at all in this part of the skull, but the whole shape itself, as smoothly reconstructed by the spline, may be used
FITTING
CONIC’ SECTIONS
'1'0 S(‘ATTI~;Kl:l,
I)hTA
71
;ts a function-valued or coefficient-valued measure in analysis of gro\vth pat,tcwlF. racial differences, evolution of t’he hominid form, and t.he l&c.
1. K. A. Paton, Conic sections in chr~~mosome analysis, Puttern. ftecuy. 2, 1970, 30. 2. Ii. li. Pat,on, Conic sections in automatic chromosome analysis, Machine Znlelliyrrlce 5, 1970. $11. :;. II. II. Biggerstaff, Three variations iu dental arch form c4nmtcd by :I quadratic eyuatioil, J. Dent. Res. 51, 1972, 1599. 4. A. Albano, Representation of digitized contours in terms of conic arcs and straight-lirrc segments, Computer Graphics and Image Processing 3, 1974, 23. 5. I). I(. Cooper and N. Yalabik, On the computational cost of approximating and recognizing n,lise-perturbed straight lines and quadratic :u-cs in the plane, Il?EE Trans. Comprrters c-25, 1976, 1020. 6. 1:. ( ;nanadesikan, .Ilethotls for Stutisticd DC&I Analysis of J/ ultizuriute Observations, Wiley, New York, 1977. 7. K. Pearson, On lines and planes of closest tit to systems of point,s in space, I’hilos. dlag. Ser. /i 2, 1901, 559. S. C. Ii. Rao, Lineur Statistical Inference and Its Applications, 2nd ed., Wiley, New York, 197:;. 9. I). J. Poiricr, Piecewise regression using cubic splines, J. Amer. Slatisl. Assoc. 68, 1973, 515. 10. S. \\‘old, Bpline functions in data analysis, Technomet. 16, 1974, 1. 11. I). 1‘. Thomas, Pseudospline interpolation for space curves, dfath. Conrp. 30, 1976, 58. 12. F. I,. Bookstein, The Jleasurement of Biological Shape and Shape Change, Springer-Veriag. NPW York, 1978. 18. J. J. Stoker, Di’erential Geometry, Wiley, New York, 19ti9. 14. H. Steinhaus, Length, shape, and area, Colloq. Math. 3, 1951, 1. I ,5. B. AIandclbrot, Fractals: Form, Chance, Dimension, Freeman, San Francisco, 1977. 16. li. Storer, General Zoology, McGralv-Hill, New York, 1951. 17. J. Burch, Freshwafer Sphaeriacean Clnms of ,Yorth Americu, ;\I:dacological Publications, II:imbirrg, Mich., 197,5.