In [1] the problem of estimating the parameters of a rotated parabola fitted to measured points in the plane was examined. The
corresponding method, also used in [2,3], is extended
here to the case of a rotated paraboloid. Fitting by such a surface occurs in computational metrology e.g. when some parabolic reflector will be checked to be a good one

H. Späth

Hrčak - Portal of scientific journals of Croatia

Mathematical Communications 6(2001), 173-179 173Least squares fitting with rotated paraboloids∗Helmuth Spa¨th§Abstract. In [1] the problem of estimating the parameters of a ro-tated parabola fitted to measured points in the plane was examined. Thecorresponding method, also used in [2, 3], is extended here to the case ofa rotated paraboloid. Fitting by such a surface occurs in computationalmetrology e.g. when some parabolic reflector will be checked to be a goodone.Key words: least squares, paraboloid, rotated paraboloid, descentmethodAMS subject classifications: 65D10, 65C60Received October 1, 2001 Accepted December 28, 20011. The modelFitting the given data(xi, yi, zi), i = 1, . . . ,m (1)with some rotated paraboloid e.g. appears in computational metrology when aparabolic reflector is measured. A paraboloid with the z-axis as a rotation axis andthe origin as the vertex is given byz = d(x2 + y2), |d| > 0 . (2)To be able to consider rotations and also for our numerical method it is moreconvenient to use the parametric formx = v cosuy = v sinuz = dv2(3)that fulfills (2). Considering a translation of the origin to (a, b, c) and rotationsA(β) in the x− z plane and B(γ) in the y − z plane we finally have xyz = 1 0 00 cos γ − sin γ0 sin γ cos γ cosβ 0 − sinβ0 1 0sinβ 0 cosβ a+ v cosub+ v sinuc+ dv2 (4)∗This paper is dedicated to Vreni and Dr. Hubert Huschke, CH-Brissago, Lago Maggiore, whosekind invitation in August 2001 enabled me to complete it.§Department of Mathematics, University of Oldenburg, Postfach 2503, D-26 111 Oldenburg,Germany, e-mail: spaeth@mathematik.uni-oldenburg.de174 H.Spa¨thwhere0 ≤ v <∞, 0 ≤ u < 2π (5)and the unknowns are a, b, c, d, β, γ.Instead of rotating the translated model (2) we prefer to rotate the given data(1). This can be done in two steps by xiyizi = 1 0 00 cos γ sinγ0 − sinγ cos γ xiyizi (i = 1, . . . ,m) (6)and  x˜iy˜iz˜i = cosβ 0 sinβ0 1 0− sinβ 0 cosβ xiyizi (i = 1, . . . ,m) . (7)For later purposes we note thatx˜′i = − sinβ xi + cosβ ziy˜′i = 0z˜′i = − cosβ xi − sinβ zi(8)(Here ′ means the derivative w.r.t. β.) andx˜′i = cosβ xi − sinβ(cos γ yi + sin γ zi)y˜′i = − sinγ yi + cos γ ziz˜′i = − sinβ xi − cosβ(cos γ yi + sin γ zi) .(9)(Here ′ means the derivative w.r.t. γ.)Now let some point on the paraboloid, i.e. (a + vi cosui, b+ vi sinui, c + dv2i )with unknown values (ui, vi) (i = 1, . . . ,m) correspond to each given and rotated(so far with unknown angles β and γ) data point (x˜i, y˜i, z˜i). Then the minimizationofS(a, b, c, d, β, γ, u1, . . . , um, v1, . . . , vm)=12m∑i=1(x˜i − a− vi cosui)2 + (y˜i − b− vi sinui)2 + (z˜i − c− dv2i )2(10)means to minimize the (half) sum of squared orthogonal distances from the ro-tated data to the unrotated paraboloid. The equivalent would be to minimize thesum of squared orthogonal distances from the original data to points on the ro-tated model (4). Anyway, we have introduced a lot of further unknowns, i.e. 2m,namely u1, . . . , um, v1, . . . , vm. But this will simplify our numerical method to bedeveloped.2. The general algorithmAt first we will discuss the algorithm for a more general case. Then we will specifyit for our problem. Let the function to be minimized be more generallyT = T (w1, . . . , wM ) ≥ C > −∞Least squares fitting with rotated paraboloids 175where the M unknowns w1, . . . , wM are numbered in some way such that thereexist N groups of variablesw1 = (w1, . . . , w1),w2 = (w1+1, . . . , w2),...wN = (wN−1+1, . . . , wN ),of sizes L − L−1 (L = 1, . . . , N, 0 = 0, N =M) with the following property: ForL = 1, . . . , N and givenw(t+1)1 ,w(t+1)2 , . . . ,w(t+1)L−1 ,w(t)L+1, . . . ,w(t)Nin the (t+1)-th iteration it should be possible to find a global minimum w∗L of thefunctionT (wL) = T (w(t+1)1 , . . . ,w(t+1)L−1 ,wL,w(t)L+1, . . . ,w(t)N ) .Then we set w(t+1)L = w∗L and proceed. (Necessary conditions for a minimum are∂T∂wL = 0, but we suppose also some means to identify a global minimum.) Theabove mentioned property would implyT (w(t)1 ,w(t)2 ,w(t)3 , . . . ,w(t)N )≥ T (w(t+1)1 ,w(t)2 ,w(t)3 , . . . ,w(t)N )≥ T (w(t+1)1 ,w(t+1)2 ,w(t)3 , . . . ,w(t)N )≥ . . .≥ T (w(t+1)1 ,w(t+1)2 ,w(t+1)3 , . . . ,w(t+1)N ) .Thus we would have a descent when moving from t to t+1. For t = 0 starting valueshave to be given. It will depend on these values to which minimum the algorithmwill converge.Now this algorithm will be used for the special objective function S of (10). Itwill turn out that due to its properties the group sizes are always one and that itwill be very easy to find global minima as desired. Just as in the general case it ispossible to choose a suitable sequence of w1, . . . ,wN in order to eventually improveconvergence.3. The algorithm for the rotated paraboloidStep 0: Let starting values a, b, c, d, β, γ be given ((ui, vi), i = 1, . . . ,m will not beneeded.)Step 1: Using β and γ as given we calculate (xi, yi, zi) and (x˜i, y˜i, z˜i) (i = 1, . . . ,m)using (6) and (7) in turn.Step 2: For each i = 1, . . . ,m the necessary condition ∂S∂ui = 0 results for vi = 0(vi = 0 makes no sense) insinui(x˜i − a)− cosui(y˜i − b) = 0 . (11)176 H.Spa¨thIf∂2S∂u2i= cosui(x˜i − a) + sinui(y˜i − b) > 0,then the minimum isui = atan(y˜i − bx˜i − a), (12)otherwise ui has to be replaced by ui + π (i = 1, . . . ,m).Step 3: For each i = 1, . . . ,m the necessary condition ∂S∂vi = 0 results (using thoseui from Step 2) in2d2v3i + (1− 2d(z˜i − c))vi − (cos ui(x˜i − a) + sinui(y˜i − b)) = 0 . (13)As d = 0, this is a third degree polynomial equation in vi that has eitherone real root or three real roots (see also [1]).In the first case the root must correspond to the unique global minimumbecause limvi→±∞ S(vi) = ∞. In the second case one has to select thatvalue out of three that minimizes the i-th term of S. (Note that S isseparable w.r.t. either ui or vi for each i = 1, . . . ,m).Step 4: The necessary condition ∂S∂β = 0 delivers (using (9))H sinβ −G sinβ = 0, (14)whereH =n∑i=1xi(a+ vi cosui) + zi(c+ dv2i ) ,G =m∑i=1zi(a+ vi cosui)− xi(c+ dv2i ) .If∂2S∂β2= H cosβ +G sinβ > 0,thenβ = atan(GH), (15)else β has to be replaced by β + π.Step 5: The necessary condition ∂S∂γ = 0 delivers (using (8))U cos γ + V sin γ = 0, (16)whereU =m∑i=1sinβ yi(a+ vi cosui)− zi(b+ v sinui) + cosβ yi(c+ dv2i ) ,V =m∑i=1sinβ zi(a+ vi cosui)− yi(b+ vi sinui) + cosβ zi(c+ dv2i ) .Least squares fitting with rotated paraboloids 177If∂2S∂γ2= −U sin γ + V cos γ > 0,thenγ = atan(−UV), (17)else γ has to be replaced by γ + π.Step 6: Using the new values for β and γ we now calculate new values for (xi, yi, zi)and (x˜i, y˜i, z˜i) (i = 1, . . . ,m) applying (6) and (7).Step 7: The necessary conditions ∂S∂a =∂S∂b =∂S∂c = 0 give in turna =1mm∑i=1(x˜i − vi cosui) , (18)b =1mm∑i=1(y˜i − vi sinui) , (19)c =1mm∑i=1(z˜i − dv2i ) . (20)These values (18), (19), and (20) correspond to global minima.Step 8: Finally ∂S∂d = 0 gives (using (20)) the global minimumd =m∑i=1v2i (z˜i − c)m∑i=1v4i(21)w.r.t. d.Step 9: Calculate the current value of S to compare it with the one in the nextiteration and compare also the values of the unknowns in two successiveiterations (e.g. relative error less than given ε). If accuracy is not sufficient,then go back to Step 2. Otherwise calculate the residuals of the fit, i.e. xiyizi−B(γ)A(β) a+ vi cosuib+ vi sinuic+ dv2i (i = 1, . . . ,m) , (22)and also the translation for the original data by abc := B(γ)A(β) abc . (23)178 H.Spa¨th4. Numerical examplesAt first we produced data points (xi, yi, zi), i = 1, . . . ,m = 20 on an unrotatedparaboloid by definingxi = 5r1, yi = 5r2, zi = .5(x2i + y2i ),where r1 and r2 are different and equally distributed pseudorandom numbers in[−1, 1] for each new i. These data were rotated by A(.5) and B(−.5) and afterwardstranslated by (a, b, c) = (1, 2, 3) to give the first data set. Four further data setswere derived similarly by adding g∗r to the new xi, yi, zi where again r ∈ [−1, 1] waspseudo-randomly varying with each component and with each i = 1, . . . ,m = 20.The number g was 0 for the first data set and g = .1, .25, .6, 1 for the four otherones, respectively. The data of all five data sets were rounded to three digits afterthe decimal point before using them. For a global minimum of S we thus wouldexpect S ≈ 0 for g = 0 and S increasing with g.To test our algorithm we used ten different starting values for (a, b, c, d, β, γ),namely (r1, r2, r3, r4, r5, r6), where rk (k = 1, . . . , 6) were different pseudorandomnumbers equally distributed in [0, 1], different for each of the five data sets and alsodifferent for each of the ten sets. The results are found in Table 1.g a b c d β γ S it ri0 1.0002 2.0002 3.0001 .5000 .5000 -.5000 .00000124 1285 450− 2400.1 .9774 1.9922 2.9793 .5044 .4978 -.5005 .07953690 870 700− 1200.25 .8283 2.0519 3.0584 .4861 .4967 -.5023 .21590963 1570 975− 2300.6 2.7480 1.9731 1.3587 .6586 .5733 -.4020 .53725064 1810 1575− 21001. 1.3638 1.4461 2.2132 .5974 .4831 -.5242 3.12157106 910 600− 1175Table 1.Astonishingly for each value of g = 0, .1, .25, .6, 1 we received for each of the tenstarting values the same value for S (thus most probably the global minimum) andalso for the unknowns a, b, c, d, β, γ, u1, . . . , um, v1, . . . , vm. The range of the num-ber of iterations ri to get four exact digits after the decimal point for all unknownsand also the corresponding average number it of iterations seem rather high at firstglance. But on the other hand, the overall computing time for all five data sets andeach time ten starting values was about one minute on a PC and thus remarkablylow. Considering the value of S this one normally was very fast decreasing duringthe first few iterations and then it took a very large number of iterations to receivethe attended four digits accuracy.5. ConclusionsThe described algorithm to fit the measured data with a rotated paraboloid canbe implemented easily and it seems to behave well with arbitrary starting values,Least squares fitting with rotated paraboloids 179though a global minimum cannot be guaranteed. The same situation is with theGAUSS-NEWTON (see [4]) or the NEWTON method where you need the Jacobianand/or the Hessian matrix, too.) Our algorithm can also be realized in a similarway e.g. for spheres [2], for ellipsoids [3], cylinders, and half cones.References[1] H.Spa¨th, An algorithm for orthogonal squared distance fitting with parabolas,in: Proceedings of the IMACS-GAMM International Symposium on NumericalMethods and Error-Bounds held at the University of Oldenburg, July 9–12,1995,(G.Alefeld and J.Herzberger, eds.), Akademie-Verlag, Berlin, 1996, 261–269.[2] H.Spa¨th, Orthogonal least squares fitting by conic sections, in: Recent Ad-vances in Total Least Squares and Errors-in-Variables Techniques, (S. vanHuffel,ed.), SIAM, Philadelphia, 1997, 259–264.[3] H.Spa¨th, Least squares fitting with spheres, J. of Optimization Theory andApplications 96(1998), 191–199.[4] D.A.Turner, I. J. Anderson, J. C.Mason, M.G.Cox, An algorithmfor fitting an ellipsoid to data, unpublished manuscript (available fromhttp://citeseer.nj.nec.com/322908.html ).

Least squares fitting with rotated paraboloids

Abstract

Similar works

Full text

Available Versions

Hrčak - Portal of scientific journals of Croatia