From: amara@quake.stanford.edu (Amara Graps) Newsgroups: sci.math.num-analysis Subject: Re: Fun and Exciting Ellipse Fitting Problem... Date: Sun, 24 Mar 1996 16:57:29 -0800 In article <315174CF.5E0BF2D2@lincoln.ac.nz>, Dave Lane wrote: > Hello, > > I've got a fun and challenging problem here that's been haunting me for > the past couple of weeks. [...] > What I want to do is use these points, as well as the angular > constraints to compute the center (c,d), angle of orientation (angle "m" > which I define as the angle between the major axis and the x-axis), and > the lengths of the major and minor axes (a & b) of the *best-fit* > ellipse. [...] Here, I post a preview of part of Steve Sullivan's April 1996 issue of his Numerical Analysis FAQ. I wrote a new section: "Fitting a circle to a set of points" several weeks ago, that can maybe give you some good pointers for your ellipse-fitting problem. The circle/ellipse fitting problem has come up here on this list 3 or 4 times in the last year. In my write up, I used answers from myself and the following 4 people: Peter Somlo somlo@zeta.org.au Pawel Gora gora@if.uj.edu.pl Zdislav V. Kovarik kovarik@mcmail.cis.McMaster.CA Daniel Pfenniger pfennige@scsun.unige.ch So here you go. Amara (Amara Graps agraps@netcom.com) --------------------------------- q230.8. Fitting a circle to a set of points Problem: Given N Points: (x1,y1); (x2,y2).. (xN,yN), we want to find for best-fit circle: (X0,Y0), R. (Note: for fitting an *ellipse*, substitute the equation for an ellipse for the equation for a circle in the "Brute Force Approach"). ------------- Brute Force Approach (leads to a non-linear system): [Amara Graps] ------------- Idea: Minimize by least squares the root-mean-squared error of the radius in the equation for a circle. In this method, one minimizes the sum of the squares of the *length* differences, which gives you a non-linear problem with all the associated pitfalls. (x-X0)^2 + (y-Y0)^2 = R^2 Equation for a circle R = SQRT[ (x-X0)^2 + (y-Y0)^2) ] Radius of the circle where: (X0,Y0) = center of circle (x,y) = point coordinates R = radius 1) Get first estimate for (X0,Y0,R). (Details: Find some points to make first estimates- either solve the circle exactly (3 equations, 3 unknowns) to get a first estimate of the center and radius, or just do a kind of centroid calculation on all points- to get a rough estimate of the center and radius.) 2) Calculate r (r1, r2,.. rN) for each of your N points from the equation for a radius of a circle. 3) Calculate the root-mean-squared error For example, for 5 given points on the circle: RMS error = SQRT[ [ (r1-R)^2 + (r2-R)^2 + (r3-R)^2 + (r4-R)^2 + (r5-R)^2] / 3 ] (dividing by "3" because we have 3 free parameters.) 4) Change (X0,Y0,R) slightly in your minimization algorithm to try a new (X0,Y0,R). (Details: Because minimization algorithms can get very computationally intensive, if one's problem is a simple one, I would look for a "canned" minimization routine. Some commercial computer programs for plotting and spreadsheets do this sort of thing. For example, the Excel spreadsheet has a built-in "solver" that will perform minimzation. Other possibilties for software: Matlab with the optimization toolbox, MACSYMA, ODRPACK from Netlib, the recent L-BFGS-B package (from ftp://eecs.nwu.edu/pub/lbfgs) which allows you to specify bounds on the variables.) 5) Calculate r (r1, r2 etc.) again from new (X0,Y0,R) above. 6) Calculate RMS error again. 7) Compare current RMS error with previous RMS error. If it doesn't vary by more some small amount, say 10^{-3} then you're done, otherwise continue steps 4 -- 7. Other References: Berman & Somlo: "Efficient procedures for fitting circles..." IEEE Instrum & Meast., IM-35, no.1, pp. 31-35, 1986. [Peter Somlo] ----------------- Other (more elegant) approaches that reduce to a linear system. ----------------- If you choose to minimize the squares of the *area* differences, you get a linear problem, which is a much safer way. [Pawel Gora, Zdislav V. Kovarik, Daniel Pfenniger, Condensed by Amara Graps] 1. Function to minimize is (sum of the area differences): Q = Sum { [ (xi - X0)^2 + (yi -Y0)^2 - R^2 ]^2 , i=1..N } 2. A necessary condition is the system of 3 equations with 3 unknowns X0, Y0, R. Calculate the partial derivatives of Q, with respect to X0, Y0, R. (all d's are partial) dQ / dX0 = 0 dQ / dY0 = 0 dQ / dR = 0 3. Developing we get the linear least-squares problem: | x1 y1 1 | | a | | -x1^2-y1^2 | | x2 y2 1 | | b | =~ | -x2^2-y2^2 | | x3 y3 1 | | c | | -x3^2-y3^2 | | x4 y4 1 | | -x4^2-y4^2 | | x5 y5 1 | | -x5^2-y5^2 | ..... ......... (for example, for 5 points) where a = 2 X0; b = 2 Y0 ; c = X0^2 + Y0^2 - R^2. Take any good least-squares algorithm to solve it, yielding a,b,c. So the final circle solution will be given with X0 = a/2; Y0 = b/2; R^2 = X0^2+Y0^2 - c. By the way with 5 points you can also find the unique quadratic form (ellipse, hyperbola) which passes exactly through 5 points. With more than 5 points one can do a linear least-squares as well. The problem is then to minimize: | x1^2-y1^2 x1y1 x1 y1 1 | | a | | -x1^2-y1^2 | | x2^2-y2^2 x2y2 x2 y2 1 | | b | =~ | -x2^2-y2^2 | | x3^2-y3^2 x3y3 x3 y3 1 | | c | | -x3^2-y3^2 | | x4^2-y4^2 x4y4 x4 y4 1 | | e | | -x4^2-y4^2 | | x5^2-y5^2 x5y5 x5 y5 1 | | f | | -x5^2-y5^2 | ..... ......... At obsftp.unige.ch in /fit there are Fortran programs fcircle.f, fellipse.f and the Lawson & Hanson least-squares routines ls.f showing how to implement these least-squares problems. ----------------- The robust or L1 or least-first-power approximation [Zdislav V. Kovarik]. ----------------- If you try to minimize W(a,b,r) = SUM(j=1:N) ABS (((x_j-a)^2 + (y_j-b)^2)^(1/2) - r) all you have to do is set up the 10 (i.e. 5 choose 3) circles passing through every choice of 3 of 5 points, calculate W(a,b,r) for each of them and pick the minimizing circle. The extra feature is that this procedure separates and discards points which are conspicuously far from the majority trend. (Of course, it becomes increasingly time-consuming when the number of given points increases.) [Amara Graps] This method of determining the minimum bounding circle from a set of _circles_ is solved, and with code available at: http://www.intergalact.com/circles.html -- *************************************************************** Amara Graps amara@quake.stanford.edu Solar Oscillation Investigations Stanford University http://quake.stanford.edu/~amara/amara.html *************************************************************** "In theory, there is no difference between theory and practice. In practice, however, there is." -- Franz Kafka ============================================================================== From: sande@haven.ios.com (Gordon Sande) Newsgroups: sci.math.num-analysis Subject: Re: Fun and Exciting Ellipse Fitting Problem... Date: Mon, 25 Mar 96 10:54:59 GMT The problem of fitting a circle or ellipse from a number of points along some portion of the boundary arrives here or in comp.graphics.algorithms every so often. The ODRPack (Orthogonal Distance Regression Package) in NetLib gives the fitting of a fragment of the ellipse as an example. There are some subtleties to the nonlinear problem which ODRPack has worked out. The graphics problem is associated with very low error in fitting to font outlines and wants to have low computational cost. There are many alternate "criteria" for judging the "fit" which produce different answers. They all match for no error and some make little sense for larger error. ODRPack is usable when the errors are larger but requires more computation. It may even be usable for the font fitting problem with current faster machines. The ODRPack example is not immediately and exactly the problem here, but ODRPack should be able to solve the nonlinear implicit equations required. Gordon Sande In Article <315174CF.5E0BF2D2@lincoln.ac.nz>, Dave Lane wrote: >Hello, > >I've got a fun and challenging problem here that's been haunting me for >the past couple of weeks. I'm writing this hoping against hope that >someone with a mind greater than my own will either take it on as an >intellectual conquest or will have dealt with this sort of thing before >and point me in the right direction. Here it is-> > >Suppose I have the following: > >From the origin of an x-y coordinate system, I have 3 rays extending. >For visualisation purposes, I'll tell you that the angle between the >outermost two is less than 45 degrees. I know the angles (and therefore >the eq'ns of each ray, y = x*sin(angle)) and I'll call them a1, a2, a3. > >In the 2D space bounded by the two outer rays there is an ellipse of >unknown size and orientation. I also don't know where its center, (c,d) >is, though I *do* know that it lies along the center ray (d = >c*sin(a2)), and that the two outer rays are tangent to it. > >To allow for a more or less definite solution, I also have a number of >(x,y) points which fall on or very near (within measurement error) the >side of the ellipse which is closest to the origin. > >What I want to do is use these points, as well as the angular >constraints to compute the center (c,d), angle of orientation (angle "m" >which I define as the angle between the major axis and the x-axis), and >the lengths of the major and minor axes (a & b) of the *best-fit* >ellipse. (I'd like to try to find a closed form least squares solution, >but until them, I'm trying to set up a working iterative calculation >inefficient though it might be.) > >For convenience, here are some other relevant expressions. I'll define >the equation of the ellipse as follows: > >p^2/a^2 + q^2/b^2 = 1 > >Note, that for the arbirarily oriented ellipse centered on a point >(c,d), p = x*cos(m)-y*sin(m)+c and q = x*sin(m)+y*cos(m)+d > >Also, if it's of any use to anyone, a parametric representation of the >same expression is > >x = a*cos(m)*cos(g) - b*sin(m)*sin(g) + c >y = a*sin(m)*cos(g) + b*cos(m)*sin(g) + d > >where g = {0 to 2*pi}. > >I just hope this problem with trigger an intellectual frenzy in one of >you addicted puzzlers out there. Please email your responses to me at >the address directly below. Thanks in advance for all of your >dedication and sleepless nights ;). >Cheers, > >Dave > >-- >+======================================================================+ >| Dave Lane dlane@lincoln.ac.nz | >| NZ Forest Research Institute | >| DNRE, Lincoln University Wk: +64 3 325 2811 x8793 | >| P.O. Box 84 Hm: 343 3192 | >| Lincoln, NZ Fx: 325 3845 | >+======================================================================+ Gordon Sande Sande & Associates, Inc. (201) 864-4564 600 Sanderling Court sande@haven.ios.com Secaucus, NJ 07094 ============================================================================== From: dwells@nrao.edu (Don Wells) Newsgroups: sci.math.num-analysis Subject: Re: Fun and Exciting Ellipse Fitting Problem... Date: 27 Mar 1996 19:30:43 GMT "GS" == Gordon Sande writes: GS> The problem of fitting a circle or ellipse from a number of points along GS> some portion of the boundary arrives here or in comp.graphics.algorithms GS> every so often. The ODRPack (Orthogonal Distance Regression Package) in GS> NetLib gives the fitting of a fragment of the ellipse as an example.. GS> The ODRPack example is not immediately and exactly the problem here, but GS> ODRPack should be able to solve the nonlinear implicit equations required. Yes, ODRPACK is able to solve this problem, except for the constaint equations. ODRPACK is particularly useful in repetitive production situations, because its Fortran subroutines can be called from application programs, and because its internal variables are visible, which permits customization. I append bibliographic references for ODRPACK. Another good approach to this class of problems is Gaussfit. Gaussfit is a standalone program which is especially convenient to use for one-time problems and exploratory work. Gaussfit solves ODR problems for nonlinear implicit equations, just as ODRPACK does. In addition, Gaussfit accepts constraint equations. The latter feature would be useful for the problem described in the query which started this thread of discussion. Furthermore, Gaussfit can operate in various "robust" estimation modes to reject outliers. See the appended reference, which contains the anonFTP URL. Gaussfit was constructed by the astrometry instrument team of the Hubble Space Telescope project. I use both ODRPACK and Gaussfit to fit three-dimensional paraboloidal surfaces to XYZ data points using orthogonal distance regression. The two programs give essentially the same answers, as you would expect. The Gaussfit model which I use is appended below. The model fits the common vertex equation of the conic sections, with numerical eccentricity 'eps' held constant (it is set to unity in an input file). If the declaration statement 'constant eps' is changed to 'parameter eps', this model should be able to fit three-dimensional ellipsoidal surfaces to XYZ data points. (For an ellipsoid fit, the 'fl' parameter produced by this model is $b^2 / {2 a}$, where b is the semiminor axis and a is the semimajor axis.) My particular problem involves a plane of symmetry, and so there are no vx, ry or rz parameters in this model, but the appropriate generalization is obvious (see ODRPACK code below). In summary, it appears to me that this 3-D model would be a solution for the 2-D problem posed in the query which started this thread if the special constraint equations were added (see the Gaussfit manual). The ODRPACK model which I use is also appended. It also is capable of fitting arbitrary 3-D conic sections (including ellipsoids), and it includes the three parameters which are omitted in the Gaussfit model. It also includes analytic derivatives, which were derived by Mathematica. ODRPACK does not require the analytical derivatives, of course, and it can even be faster to command it to use the usual numerical derivative approximation when the expressions are as complicated as in this case. (Note: Gaussfit contains its own symbolic algebra facilities, which it uses to do analytic differentiation of its model and constraint equations.) -Don Wells @Manual{jfmm88, title = "User's Manual -- GaussFit: {A} system for least squares and robust estimation", author = "William H. Jefferys and Michael J. Fitzpatrick and Barbara E. McArthur and James E. McCartney", organization = "The {U}niversity of {T}exas at {A}ustin", edition = "1.0-12/2/88", year = 1988, month = "December", note = "See URL=\verb|ftp://clyde.as.utexas.edu/pub/gaussfit/|" } @Manual{bbrs92f, title = "User's Reference Guide for {ODRPACK} Version 2.01---Software for Weighted Orthogonal Distance Regression", author = "Paul T. Boggs and Richard H. Byrd and Janet E. Rogers and Robert B. Schnabel", organization = "National Institute of Standards and Technology", address = "Gaithersburg, MD 20899", year = 1992, month = "June", note = "ODRPACK is a software package for {\em weighted orthogonal distance regression}, i.e., for finding the set of parameters that minimize the sum of the weighted orthogonal distances from a set of observations to the curve or surface determined by the parameters. It can also be used to solve the nonlinear ordinary least squares problem. See URL=\verb|http://netlib.att.com/magic/netlib_find?db=0&pat=odrpack| (or URL=\verb|ftp://netlib.att.com/netlib/odrpack/|). Also see \cite{bbs87} and \cite{bdbs89}. This is report NISTIR 92-4834." } @Article{bbs87, author = "Paul T. Boggs and Richard H. Byrd and Robert B. Schnabel", title = "A Stable and efficient algorithm for nonlinear orthogonal distance regression", journal = "SIAM Journal on Scientific and Statistical Computing", year = 1987, volume = 8, number = 6, pages = "1052-1078", month = "November" } @Article{bdbs89, author = "P. T. Boggs and J. R. Donaldson and R. H. Byrd and R. B. Snabel", title = "Algorithm 676, {ODRPACK} -- Software for Weighted orthogonal distance regression", journal = "ACM Transactions On Mathematical Software", year = 1989, volume = 15, pages = "348-364" } -=-=-=-=- cut here -=-=-=-=- /* GaussFit program to fit paraboloids to the distorted GBT surface */ /* Don Wells, NRAO-CV, 08-December-1994,3/16/95 */ /* -=-=-=-=-=-=-=-=- Begin Copyright Notice -=-=-=-=-=-=-=-=-=-=-=-=-=- Copyright (C) 1995 Associated Universities, Inc. Washington DC, USA. This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. Correspondence concerning GBT software should be addressed as follows: GBT Operations, National Radio Astronomy Observatory, P. O. Box 2, Green Bank, WV 24944-0002 USA -=-=-=-=-=-=-=-=-=- End Copyright Notice -=-=-=-=-=-=-=-=-=-=-=-=-=- */ observation x; /* + = away from plane of symmetry */ observation y; /* + = away from feed arm */ observation z; /* + = up, height of surface */ parameter vy; /* vertex position */ parameter vz; parameter rx; /* paraboloid tilt about x-axis (radians) */ parameter fl; /* focal length */ constant eps; /* eccentricity */ main() { variable x0, y0, z0, srx, crx, xp, yp, zp, xpp, ypp, zpp; x0 = 0.0; y0 = vy; z0 = vz; srx = sin(rx); crx = cos(rx); /* COMPUTE PREDICTED VALUES */ while (import()) { xp = x - x0; yp = y - y0; zp = z - z0; /* positive rx is right-hand-rule about +X: */ xpp = xp ; ypp = +crx * yp +srx * zp; zpp = -srx * yp +crx * zp; /* Common vertex equation of the conics (VNR, p.315): */ export( (xpp^2 + ypp^2) - (((4.0 * fl) - (1.0 - eps^2) * zpp) * zpp) ); } } -=-=-=-=- cut here -=-=-=-=- SUBROUTINE FCN(N,M,NP,NQ, + LDN,LDM,LDNP, + BETA,XPLUSD, + IFIXB,IFIXX,LDIFX, + IDEVAL,F,FJACB,FJACD, + ISTOP) c implicit undefined (a-z) c C SUBROUTINE ARGUMENTS C ==> N NUMBER OF OBSERVATIONS C ==> M NUMBER OF COLUMNS IN EXPLANATORY VARIABLE C ==> NP NUMBER OF PARAMETERS C ==> NQ NUMBER OF RESPONSES PER OBSERVATION C ==> LDN LEADING DIMENSION DECLARATOR EQUAL OR EXCEEDING N C ==> LDM LEADING DIMENSION DECLARATOR EQUAL OR EXCEEDING M C ==> LDNP LEADING DIMENSION DECLARATOR EQUAL OR EXCEEDING NP C ==> BETA CURRENT VALUES OF PARAMETERS C ==> XPLUSD CURRENT VALUE OF EXPLANATORY VARIABLE, I.E., X + DELTA C ==> IFIXB INDICATORS FOR "FIXING" PARAMETERS (BETA) C ==> IFIXX INDICATORS FOR "FIXING" EXPLANATORY VARIABLE (X) C ==> LDIFX LEADING DIMENSION OF ARRAY IFIXX C ==> IDEVAL INDICATOR FOR SELECTING COMPUTATION TO BE PERFORMED C <== F PREDICTED FUNCTION VALUES C <== FJACB JACOBIAN WITH RESPECT TO BETA C <== FJACD JACOBIAN WITH RESPECT TO ERRORS DELTA C <== ISTOP STOPPING CONDITION, WHERE C 0 MEANS CURRENT BETA AND X+DELTA WERE C ACCEPTABLE AND VALUES WERE COMPUTED SUCCESSFULLY C 1 MEANS CURRENT BETA AND X+DELTA ARE C NOT ACCEPTABLE; ODRPACK SHOULD SELECT VALUES C CLOSER TO MOST RECENTLY USED VALUES IF POSSIBLE C -1 MEANS CURRENT BETA AND X+DELTA ARE C NOT ACCEPTABLE; ODRPACK SHOULD STOP C INPUT ARGUMENTS, NOT TO BE CHANGED BY THIS ROUTINE: INTEGER I,IDEVAL,ISTOP,L,LDIFX,LDM,LDN,LDNP,M,N,NP,NQ DOUBLE PRECISION BETA(NP),XPLUSD(LDN,M) INTEGER IFIXB(NP),IFIXX(LDIFX,M) C OUTPUT ARGUMENTS: DOUBLE PRECISION F(LDN,NQ),FJACB(LDN,LDNP,NQ),FJACD(LDN,LDM,NQ) c c -=-=-=-=-=-=-=- Begin Copyright Notice -=-=-=-=-=-=-=-=-=-=-=-=-=-=- c c Copyright (C) 1995 Associated Universities, Inc. Washington DC, USA. c c This program is free software; you can redistribute it and/or modify c it under the terms of the GNU General Public License as published by c the Free Software Foundation; either version 2 of the License, or c (at your option) any later version. c c This program is distributed in the hope that it will be useful, but c WITHOUT ANY WARRANTY; without even the implied warranty of c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU c General Public License for more details. c c You should have received a copy of the GNU General Public License c along with this program; if not, write to the Free Software c Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. c c Correspondence concerning GBT software should be addressed as follows: c GBT Operations, National Radio Astronomy Observatory, P. O. Box 2, c Green Bank, WV 24944-0002 USA c c -=-=-=-=-=-=-=- End Copyright Notice -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- c real*8 x0, y0, z0, phi, theta, psi, fl, eps, & sphi, cphi, stheta, ctheta, spsi, cpsi, & x, y, z, xp, yp, zp, xpp, ypp, zpp c x0 = beta(1) y0 = beta(2) z0 = beta(3) phi = beta(4) theta = beta(5) psi = beta(6) fl = beta(7) eps = beta(8) c sphi = sin(phi) cphi = cos(phi) stheta= sin(theta) ctheta= cos(theta) spsi = sin(psi) cpsi = cos(psi) c C COMPUTE PREDICTED VALUES IF (MOD(IDEVAL,10).GE.1) THEN DO 110 L = 1,NQ DO 100 I = 1,N c x = xplusd(i,1) y = xplusd(i,2) z = xplusd(i,3) c xp = x - x0 yp = y - y0 zp = z - z0 c xpp = xp +spsi * yp -stheta * zp ypp = -spsi * xp + yp +sphi * zp zpp = stheta * xp -sphi * yp + zp c c Common vertex equation of the conics (VNR, p.315): f(i,l) = (xpp**2 + ypp**2) & - (((4d0 * fl) - (1d0 - eps**2) * zpp) * zpp) c 100 CONTINUE 110 CONTINUE END IF c C COMPUTE DERIVATIVES WITH RESPECT TO BETA IF (MOD(IDEVAL/10,10).GE.1) THEN DO 210 L = 1,NQ DO 200 I = 1,N c x = xplusd(i,1) y = xplusd(i,2) z = xplusd(i,3) c xp = x - x0 yp = y - y0 zp = z - z0 c fjacb(i,1,l) = & 2d0*spsi*(yp + zp*sphi - xp*spsi) - & (1d0 - eps**2)*stheta* & (zp - yp*sphi + xp*stheta) - & 2d0*(xp + yp*spsi - zp*stheta) + & stheta*(4d0*fl - (1d0 - eps**2)* & (zp - yp*sphi + xp*stheta)) fjacb(i,2,l) = & -2d0*(yp + zp*sphi - xp*spsi) + & (1d0 - eps**2)*sphi*(zp - yp*sphi + & xp*stheta) - & 2d0*spsi*(xp + yp*spsi - zp*stheta) - & sphi*(4d0*fl - (1d0 - eps**2)* & (zp - yp*sphi + xp*stheta)) fjacb(i,3,l) = & 4d0*fl - 2d0*sphi*(yp + zp*sphi - xp*spsi) - & 2d0*(1d0 - eps**2)*(zp - yp*sphi + xp*stheta) + & 2d0*stheta*(xp + yp*spsi - zp*stheta) fjacb(i,4,l) = & 2d0*zp*cphi*(yp + zp*sphi - xp*spsi) - & (1d0 - eps**2)*yp*cphi* & (zp - yp*sphi + xp*stheta) + & yp*cphi*(4d0*fl - & (1d0 - eps**2)*(zp - yp*sphi + xp*stheta)) fjacb(i,5,l) = & (1d0 - eps**2)*xp*ctheta* & (zp - yp*sphi + xp*stheta) - & 2d0*zp*ctheta*(xp + yp*spsi - & zp*stheta) - & xp*ctheta*(4d0*fl - & (1d0 - eps**2)*(zp - yp*sphi + xp*stheta)) fjacb(i,6,l) = & -2d0*xp*cpsi*(yp + zp*sphi - xp*spsi) + & 2d0*yp*cpsi*(xp + yp*spsi - zp*stheta) fjacb(i,7,l) = & -4d0*(zp - yp*sphi + xp*stheta) fjacb(i,8,l) = & -2d0*eps*(zp - yp*sphi + xp*stheta)**2 c 200 CONTINUE 210 CONTINUE END IF c C COMPUTE DERIVATIVES WITH RESPECT TO DELTA IF (MOD(IDEVAL/100,10).GE.1) THEN DO 310 L = 1,NQ DO 300 I = 1,N c x = xplusd(i,1) - x0 y = xplusd(i,2) - y0 z = xplusd(i,3) - z0 c xp = x - x0 yp = y - y0 zp = z - z0 c fjacd(i,1,l) = & -2d0*spsi*(yp + zp*sphi - xp*spsi) + & (1d0 - eps**2)*stheta* & (zp - yp*sphi + xp*stheta) + & 2d0*(xp + yp*spsi - zp*stheta) - & stheta*(4d0*fl - (1d0 - eps**2)* & (zp - yp*sphi + xp*stheta)) fjacd(i,2,l) = & 2d0*(yp + zp*sphi - xp*spsi) - & (1d0 - eps**2)*sphi*(zp - yp*sphi + & xp*stheta) + & 2d0*spsi*(xp + yp*spsi - zp*stheta) + & sphi*(4d0*fl - (1d0 - eps**2)* & (zp - yp*sphi + xp*stheta)) fjacd(i,3,l) = & -4d0*fl + 2d0*sphi*(yp + zp*sphi - xp*spsi) + & (1d0 - eps**2)*(zp - yp*sphi + xp*stheta) - & (-1d0 + eps**2)*(zp - yp*sphi + xp*stheta) - & 2d0*stheta*(xp + yp*spsi - zp*stheta) c 300 CONTINUE 310 CONTINUE END IF c RETURN END -- Donald C. Wells Associate Scientist dwells@nrao.edu http://fits.cv.nrao.edu/~dwells National Radio Astronomy Observatory +1-804-296-0277 520 Edgemont Road, Charlottesville, Virginia 22903-2475 USA ============================================================================== From: ph+@cs.cmu.edu (Paul Heckbert) Newsgroups: sci.math Subject: Re: best fitting a circle Date: 7 Jun 1996 18:39:31 GMT In article <4p9ihh$dpn@earth.alpha.net>, Richard Lukas wrote: >I have a file containing approximately 100 (x,y) coordinates. I know >that these points are sampled from a semicircular arc. I need a >procedure that will allow me to find the radius and centerpoint >of the circle these points fall on. I recommend the following paper @ARTICLE{Pratt87, author = "Vaughan Pratt", title = "Direct Least-Squares Fitting of Algebraic Surfaces", pages = "145-152", journal = "Computer Graphics (SIGGRAPH '87 Proceedings)", volume = "21", number = "4", year = "1987", month = "July", } ============================================================================== From: Amara Graps Newsgroups: sci.math,sci.physics Subject: Re: Fitting a circle to a set of points Date: Fri, 16 Aug 1996 18:11:05 -0700 Paul Knowles wrote: > Hi Folks, > > I'm working with some tracking algorithms for particles in a > magnetic field, and one of the things i have to optimize is > the fit of several (x,y) pairs to the unique circle through > the points. [...] Sure. There are a many ways. From the Numerical Analysis FAQ, available at: ftp://rtfm.mit.edu/pub/usenet/news.answers/num-analysis/faq/part1 Hope this helps. Amara ["q230.8" quoted from above -- djr] ***************************************************************** Amara Graps amara@quake.stanford.edu Solar Oscillation Investigations Stanford University http://quake.stanford.edu/~amara/amara.html ***************************************************************** ============================================================================== From: gerry@me.pdx.edu (Gerald Recktenwald) Newsgroups: sci.stat.math,sci.math.num-analysis,aus.mathematics Subject: Re: Help! Need Best-Fit Circle Equation. Date: Thu, 02 Apr 1998 21:27:16 -0800 In article <3523A24A.F9C10F3B@conversionsinc.com>, Kurt Feldbush wrote: > I am looking for an equation or algorithm to find the best fit >circle or ellipse to a set of data points. I have searched everywhere I >could think of with no luck. Can anyone help? Have a look at "Least-squares fitting of circles and ellipses", W.Gander, G.H. Golub, and R. Strebel, BIT, vol 34 (1994), pp. 448-578 Gander et al compare different methods. MATLAB source codes are available at ftp.inf.etch.ch. Gerry -- +-+-+-+-+-+-+-+-... Gerald Recktenwald (503) 725-4296 gerry@me.pdx.edu http://www.me.pdx.edu/~gerry/ PSU Mechanical Engineering, P.O. Box 751, Portland, OR 97207 ============================================================================== From: islam Newsgroups: sci.math.num-analysis Subject: Re: least-square-fit for an ellipse Date: Tue, 28 Jul 1998 11:56:32 +0100 Hi - This has already been done a number of times. for a good reveiw of the subject+code see: http://vision.dai.ed.ac.uk/maurizp/ElliFitDemo/Paper/icip.html http://www.dai.ed.ac.uk/students/maurizp/ElliFitDemo/demo.html also available on the web somewhere: Book Author/Editor :: Gander W; Golub GH; Strebel R Book Title :: Fitting Circles and Ellipses; Least Squares Solutions Publisher :: Departement Informatik. Instititut fur Wissenschaftliches Rechnen, ETH Zurich ____________________________________________ Suhail A Islam Biomolecular Modelling Laboratory Imperial Cancer Research Fund, P.O. Box 123 44 Lincoln's Inn Fields, London WC2A 3PX Tel: (0171) 269 3380 Fax: (0171) 269 3258 email: islam@icrf.icnet.uk http://www.icnet.uk/bmm/ ____________________________________________ > > I need to know an easy(!) way to calculate the parameters of an ellipse (the > two focal points and the two axis) by using least-square-fit. Assume there > are enough given points on the circumference of the ellipse. > > Here are the two equations for an ellipse I tried to get a solution with: > A*x^2 + B*y^2 + C*x + D*y + E = 0 > and > (x-x0)^2 / B^2 + (y-y0)^2 / B^2 = 1 > > As I've to implement a solution in c, any algorithms in c are welcome. > > Thanks for your efforts, > > Dennis