We consider fitting data points in space by a set of two concentric
spheres. This problem ought to occur within computational metrology. A heuristic algorithm is developed and its efficiency is demonstrated by some numerical example

H. Späth

Hrčak - Portal of scientific journals of Croatia

Mathematical Communications 11(2006), 63-68 63Data fitting with a set of two concentric spheresHelmuth Spa¨th∗Abstract. We consider fitting data points in space by a set oftwo concentric spheres. This problem ought to occur within computa-tional metrology. A heuristic algorithm is developed and its efficiency isdemonstrated by some numerical example.Key words: data fitting, two concentric spheresAMS subject classifications: 65D10Received June 8, 2006 Accepted June 13, 20061. Statement of the problemOne standard problem in computational metrology [1, Problem C6] is to find formeasured points in space two concentric spheres of minimum difference in radii thatcontain all the points between them. But assuming all points to have errors it oughtto be more adequate to determine two spheres such that a suitable sum of squarederrors attains a minimum where simultaneously each point has to be associatedwith one of the spheres.A set of two concentric spheres is given by one common center (a, b, c) and tworadii r and R with R > r > 0. Some parametric representation for the smallersphere is given byx(u, v) = a+ r cosu sin v ,y(u, v) = b+ r sinu sin v , (1)z(u, v) = c+ r cos v, 0 ≤ u < 2π, 0 < v ≤ πand the larger one is obtained when replacing r by R.2. Smallest distance of some given point to some sphereTo find the minimal squared distance p2j = p2j(u, v) of some given point (xj , yj , zj),j = 1, . . . ,m, onto the smaller sphere requires to minimizep2j(u, v) = (a+ r cosu sin v − xj)2 + (b + r sinu sin v − yj)2+(c+ r cos v − zj)2 (2)∗Department of Mathematics, University of Oldenburg, 26 111 Oldenburg, Germany, e-mail:spaeth@mathematik.uni-oldenburg.de64 H. Spa¨thwith respect to u and v. The corresponding function q2j = q2j (u, v) for the largersphere is received by replacing r by R in (2). The necessary conditions for a mini-mum of (2) (and similarly for q2j ) are12r∂p2j∂u= sin v(cos u(yj − b)− sinu(xj − a)) = 0 , (3)12r∂p2j∂v= cos v(cos u(xj − a) + sinu(yj − b))− sin v(zj − c) = 0 . (4)For (3) it results eithersin v = 0 (implying v = π and p2j constant) (5)ortg u =yj − bxj − a . (6)From (4) it resultstg v =cosu(xj − a) + sinu(yj − b)zj − c . (7)Both expressions (6) and (7) have two solutionsu = atan(yj − b)(xj − a) and u := u+ π , (8)v = atan(cosu(xj − a) + sinu(yj − b)zj − c)and v := v + π, (9)stemming from the fact that e.g.tg u =sinucosu=sin(u + π)cos(u+ π)=− sinu− cosu ,i.e. we have to choose the right signs of sinu and cosu to get the absolute minimum.As expressions (6) and (7) do not depend on r, they are the same for q2j .The decision process is as follows [2]. For j = 1, . . . ,m(i) Use (6) to determine u.(ii) Calculate suj = sinu, cuj = cosu.(iii) If (yj − b) < 0, set suj = −suj.(iv) If (xj − a) < 0, set cuj = −cuj.(v) Use (7) to determine v.(10)(vi) Calculate svj = sin v, cvj = cos v.(vii) If (cuj(xj − a) + suj(yj − b)) < 0, set svj = −svj.(viii) If (zj − c) < 0, set cvj = −cvj.(ix) Store the resulting right values into arrays, i.e. set su(j) = suj, cu(j) = cuj,sv(j) = svj, cv(j) = cvj, and use these values corresponding to the absoluteminima of p2j and q2j , respectively, when evaluating these expressions.Data fitting with a set of two concentric spheres 653. The objective functionIn order to set up a reasonable least squares principle we try to find two subsets Vand W of {1, . . . ,m} with V ∩W = ∅, V ∪W = {1, . . . ,m}, |V | ≤ 4, |W | ≤ 4 suchthat the sum of squared minimal distancesS(V,W, a, b, c, r, R) =∑j∈Vminu,vp2j(u, v) +∑j∈Wminu,vq2j (v, v) (11)is minimized. Here V is the set of given points associated with the smaller sphereandW is that one for the larger sphere. Thus we have a mixed combinatorial (V,W )and continuous (a, b, c, r, R) minimization problem without an analytic solution.4. Some heuristic algorithmWe propose the following iterative method where (beginning with suitable startingvalues for a, b, c, r, R) alternatively (V,W ) and (a, b, c, r, R) are improved such thatS is descending. This will not guarantee convergence to some local or even to aglobal minimum. But numerical experience shows that the method will normallywork in the desired sense, i.e. will find a global minimum.Step 1 [One special way of estimating starting values].Set t = 0 (iteration counter) and S(0) =∞. For an initial center(a(0), b(0), c(0)) we propose the means of all given points, i.e.a(0) = x =1mm∑j=1xj , b(0) = y =1mm∑j=1yj , c(0) = z =1mm∑j=1zj . (12)Then we calculateh2 = minj((xj − a(0))2 + (yj − b(0))2 + (zj − c(0))2),(13)H2 = maxj((xj − a(0))2 + (yj − b(0))2 + (zj − c(0))2),i.e. the minimal and the maximal squared distance of all the given points (j =1, . . . ,m) to the defined starting center. If we now take r(0) = h and R(0) = H asestimates, then all given points would lie between the inner and the outer sphere.To avoid this (because not all points with errors are estimated to be such) we reduceR(0) and increase r(0) by settingR(0) = f ∗H, r(0) = 1f∗ h with f < 1 (14)and additionally f such that R(0) > r(0).Step 2 [Determination of V (t) and W (t)].Set V (t) = W (t) = ∅. For j = 1, . . . ,m we use su(j), cu(j), sv(j), cv(j) from(10) depending on (a(t), b(t), c(t)) to calculate p2j and q2j via (2). If p2j < q2j , thenj is adjoined to V (t) otherwise to W (t). Also |V | = |V (t)|, |W | = |W (t)|, andS(t) = S(V (t),W (t), a(t), b(t), c(t), r(t), R(t)) are calculated.66 H. Spa¨thStep 3 [Updating of a, b, c, r, R].For fixed V = V (t) and W =W (t) and corresponding values for su(j), cu(j), sv(j),cv(j), j = 1, . . . ,m from Step 2 we will minimizeS(V,W, a, b, c, r, R) =∑j∈Vmin p2j +∑j∈Wmin q2j (15)w.r.t. a, b, c, r, R. The necessary conditions accordingly are∂S∂a=∂S∂b=∂S∂c=∂S∂r=∂S∂R= 0and result in the following linear systems of five equations and five unknowns a, b, c, r, R:m 0 0 F (V ) F (W )0 m 0 G(V ) G(W )0 0 m H(V ) H(W )F (V ) G(V ) H(V ) |V | 0F (W ) G(W ) H(W ) 0 |W |abcrR=mxmymzE(V )E(W ), (16)whereF (V ) =∑j∈Vcu(j) sv(j) ,G(V ) =∑j∈Vsu(j) sv(j) ,(17)H(V ) =∑j∈Vcv(j) ,E(V ) =∑j∈Vxj cu(j) sv(j) + yj su(j) sv(j) + zj cv(j) .The coefficient matrix in (16) is symmetric. If it is nonsingular (empirically true)there will be a unique solution of (16) and we can proceed.Step 4. We set t := t + 1, a(t) = a, b(t) = b, c(t) = c, r(t) = r, R(t) = R andgo back to Step 2 if there was some descent for S or if t is below some given upperbound. Otherwise stop.5. Some numerical exampleTo get test data we started with (a, b, c) = (1, 2, 3), r = 4, R = 5. Then wegenerated 8 pointsxj = a+ r cosu sin vyj = b+ r sinu sin v (j = 1, . . . , 8)zj = c+ r cos vData fitting with a set of two concentric spheres 67where for each j a new pair of pseudo random numbers (u, v) with u equally distrib-uted in [0, 2π) and v equally distributed in (0, π] was used. The resulting data pointswere disturbed by rounding them to the next integers. Similarly for j = 9, . . . , 16the same procedure was done with R instead of r. Putting together these datapoints we gotx 1 -1 0 -1 2 2 4 5 1 -1 0 -1 3 2 5 6y 2 1 3 4 2 3 1 1 2 1 3 5 2 3 1 1z -1 0 -1 6 -1 -1 6 4 -2 -1 -2 7 -2 -2 6 5Thus it can be expected that V = {1, . . . , 8} andW = {9, . . . , 16} should resultand also (a, b, c) near (1, 2, 3) and (r,R) near (4, 5). We used for the factor f inStep 1 four values f = .7, .8, .9, 1. As the results were identical and also the course ofthe iterations was very similar, we will give details only for f = .8. To characterizeV and W we introduce an integer vector p with p(j) = 1 if j ∈ V and p(j) = 2 ifj ∈W (j = 1, . . . ,m = 16).t p0–30 1 1 1 2 1 1 2 2 1 1 1 2 1 1 2 231–32 1 1 1 2 1 1 2 1 1 1 1 2 1 1 2 233 1 1 1 1 1 1 2 1 1 1 1 2 2 1 2 234 1 1 1 1 1 1 2 1 1 1 2 2 2 1 2 235–138 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2Table 1.t s a b c r R0 6.5376 1.6875 2.1875 1.3125 2.9263 5.512533 3.46524 1.4118 2.1721 2.3017 4.0188 5.240334 2.53933 1.4017 2.7472 2.5183 4.1544 5.179335 1.56106 1.4035 2.7647 2.6623 4.1878 5.090550 .18400 1.4569 2.7644 3.0095 4.1580 5.136575 .18354 1.4498 2.7474 3.0113 4.1562 5.1360100 .18347 1.4470 2.7404 3.0120 4.1555 5.1358138 .18345 1.4456 2.7369 3.0124 4.1551 5.1357Table 2.Table 1 shows the change of the vector p during the iterations. Thus we startedwith V = {1, 2, 3, 5, 6, 9, 10, 11, 13, 14} and W = {4, 7, 8, 12, 15, 16} and ended up(as expected) with V = {1, . . . , 8} and W = {9, . . . , 16}. Table 2 shows for t = 0the initial values after Step 1 and the slow changes after t = 50. Other examplesshowed a similar behaviour. Sometimes the initial V and W were so good thatthose sets did not change through the iterations on (a, b, c, r, R).68 H. Spa¨thReferences[1] V. Srinivasan, How tall is the pyramid of Cheops . . . And other problems incomputational metrology, SIAM News 29(1996), 8–9.[2] H. Spa¨th, Least-squares fitting with spheres, J. Opt. Theory and Appl.96(1998), 191–199.

Data fitting with a set of two concentric spheres

Abstract

Similar works

Full text

Available Versions

Hrčak - Portal of scientific journals of Croatia