Sphere Fitting with Applications to Machine Tracking

: We suggest a provable and practical approximation algorithm for ﬁtting a set P of n points in R d to a sphere. Here, a sphere is represented by its center x ∈ R d and radius r > 0. The goal is to minimize the sum ∑ p ∈ P | (cid:107) p − x (cid:107) − r | of distances to the points up to a multiplicative factor of 1 ± ε , for a given constant ε > 0, over every such r and x . Our main technical result is a data summarization of the input set, called coreset, that approximates the above sum of distances on the original (big) set P for every sphere. Then, an accurate sphere can be extracted quickly via an inefﬁcient exhaustive search from the small coreset. Most articles focus mainly on sphere identiﬁcation (e.g., circles in 2 D image) rather than ﬁnding the exact match (in the sense of extent measures), and do not provide approximation guarantees. We implement our algorithm and provide extensive experimental results on both synthetic and real-world data. We then combine our algorithm in a mechanical pressure control system whose main bottleneck is tracking a falling ball. Full open source is also provided.


Introduction
Approximating or fitting a set of points in a d-dimensional space by a simple shape is a fundamental problem in geometric optimization with many applications in robotics [1][2][3], computational geometry [4][5][6], machine learning [7,8], computer vision [9,10], image processing [11], computational metrology, etc. A tangible example is finding the exact location of the flaming ring as shown in Figure 1.
A possible model, which we use in this paper, is to minimize the sum of distances (called the 1 -norm of the distance vectors). The sum of distances is known to be more robust to noise than the sum of squared distances. For example, the number that minimizes the sum of distances to a set of n numbers is its median, and increasing one of the input numbers to infinity would not change this solution. However, this is not the case for the mean, which minimizes the sum of squared distances.
For every x ∈ R d and r > 0 we define to be the measure function for median-sphere. Most existing techniques are heuristics with no provable bounds on the approximation error, or have some understanding of the input points.
In this paper, we consider the special case of approximating points by a sphere, and hope to generalize the technique to other shapes in the future. Let x * , r * ∈ arg min x∈R d ,r>0 ∑ p∈P | p − x − r | represent an optimal median-sphere. Our goal is to compute x and r such that ν(P, x, r) ≤ (1 + ε)ν(P, x * , r * ).
(1) Figure 1. An RGB image that captured a flame of fire [12]. The goal is to approximate the center of the flaming ring.

Thinnest Sphere Annulus
For fitting a sphere to a set of points, the acceptable compromise is the approximation of the thinnest sphere annulus (called the ∞ -norm of the distance vectors), that minimizes maximum of distances to the input set, max p∈P | p − x − r | over every r ∈ R and x ∈ R d .
The main drawback of this optimization function is sensitivity to a single point (outlier) in the input data set that may completely change the optimal solution even compared to the sum of squared distances that was discussed in the previous section.

Median-Sphere
A better fitting solution is the median-sphere (called the 1 -norm of the distance vectors), that minimizes maximum of distances to the input set, max p∈P | p − x − r | over every r ∈ R and x ∈ R d . Although there is a lot of work for the corresponding ∞ problem, less is known for the 1 case.

Least Mean Square
Another measure is the least mean square, that minimizes sum of square of distances to the input set, ∑ p∈P p − x 2 over every x ∈ R d .

Heuristics for Circles Detection in 2D Images
Most significant work has been done for two Dimensional data in the context of image processing. Extensive work tends to compute circles for the purpose of object detection in images. However, most of these heuristics do no provide approximation guarantees and suitable only for d = 2.
A common heuristic to fit points to a circle is the well-known Circle Hough transform (CHT) [13]. The approach of this technique is to find imperfect instances of spheres by a voting procedure. This voting rates candidates of possible spheres (i.e., center and radius). The main drawback of this technique is its inefficient running time of complexity O(n d ).
In [14], a circle detection was applied by Learning Automata (LA), which is a heuristic method to solve complex multi-modal optimization problems. The detection process is considered a multi-modal optimization problem, allowing the detection of multiple circular shapes through only one optimization procedure. The drawbacks of this algorithm is that it is designated to 2D images.
A variant of a Hough transform algorithm, called EDCircles [15], aims to approximate points by arcs which are based on their Edge Drawing Parameter Free (EDPF). The general idea of the algorithm is to extract line segments in an image, convert them into circular arcs and then combine these arcs to detect circles, which is performed by an arc join heuristic algorithm with no approximation guarantees. Moreover, its running time is still inefficient when dealing with big data and large dimensional space (i.e., O(n d )) in the worst case.
In [16], a circle detector based on region-growing of gradient and histogram distribution of Euclidean distances is presented. Region-growing of the gradient is applied to generate arc support regions from a single point. The corresponding square fitting areas are defined to accelerate detection and decrease storage. A histogram is then used to count frequency of the distances that participates in the accumulator and the parameters of each circle are acquired.
In [17], the proposed detection method is based on geometric property and polynomial fitting in polar coordinates instead of Cartesian coordinates. It is tailored for two-dimensional LIDAR (laser radar sensors) data. They use Support Vector Machines (SVM) to detect the target circular object from natural lidar coordinates features representing segments from the image. For specific examples, they measure up to 99.79% detections with execution time of 16.9 milliseconds per image.
A circle detection algorithm that is based on a random sampling of isosceles triangles (ITs) is presented in [18]. This algorithm provides a distinctive probability distribution for circular shapes using a small number of iterations. Although they introduced an accurate and robust detection at sharp images, it does not handle corrupted and moving images and does not have provable guarantees.
Another variant of the Hough transform algorithm [19], called Vector Quantization (VQHT), tries to detect circles by decomposing their edges in the image into many sub-images by using Vector Quantization (VQ) algorithm. The edge points that reside in each sub-image are considered as one circle candidate group. Then the VQHT algorithm is applied for fast circle detection. In [20,21], a randomized iterative work-flow is suggested, which exploits geometrical properties of isophotes in the image to select the most meaningful edge pixels and to classify them in subsets of equal isophote curvature.
In [22], the authors propose a method to recognize a circular form using a geometric symmetry through a rotational scanning process.
The random sample consensus (RANSAC) algorithm introduced by Fishler and Bolles in [23] is a widely used robust estimator that has become a standard in the field of computer vision. The RANSAC algorithm proceeds as follows: Repeatedly, subsets of the input data are randomly selected, and model parameters (i.e., optimal median-sphere) fitting these subsets are computed. In the second step, the quality of the sphere is evaluated on the input data. The output of the algorithm is the sphere that has the highest match of all the spheres.
In [24,25] the authors presents variants of probabilistic and randomized sampling to determine sphere in 3D points over clouds.

Coreset
One of the classical techniques in developing approximation algorithms is the extraction of "small" amount of "most relevant" information from the given data (called a coreset), and then applying traditional optimal algorithms on this extracted data which ensure provable approximation results. The exact definition of coreset changes from paper to paper.
Moreover, generic techniques enable coreset maintenance of streaming, distributed and dynamic data. A survey on coresets is given in e.g., [26,27].
In [28] the author suggests a (1 + ε) approximation for the median-sphere that takes O(n + log f (d) (n) + (1/ε) f (d) ) time for constant integer d ≥ 1, where f (d) is a function of d.
In [29] the authors suggest (1 + ε)-approximation for the thinnest-sphere annulus fitting in O(nd + 1/ε 3d ) time using linearization technique and ε-kernel of a point set P, for every d ≥ 1.

Coreset Definition
In this paper, a coreset for an input set of points is a weighted subset, such that the sum of distances to every given sphere from the coreset and the original input is the same, up to multiplicative (1 + ε) factor. Here, ε ∈ (0, 1) is an error parameter, and the size of the coreset usually depends on ε, but not always [30].

Our Main Contribution
We suggest an algorithm that computes an ε-coreset for the family of spheres in R d for the 1 -norm of distances. The cardinality of the coreset is O(d 2 log 2 n/ε 2 ). See Theorem 2 for details.
Using existing algorithms, we can then extract a (1 + ε)-approximation to the optimal sphere of the original input in additional (d log n) O(1) time. That is, the algorithm outputs (1 + ε)-approximation for (1), i.e., a pair x ∈ R d and r ≥ 0 such that We highlight the strengths of the technique by implementing and testing it for different values of number n of points in the input set and for distorted and noisy input points. Our solution holds for any dimension, and provides several examples for the case of points in two and three dimensional space.
In theory, the price for using our approach is that we get only an approximated solution to the desired shape in the worst case by running the optimal solution on the coreset. In practice, we do not compute the optimal solution also on the original data, but run heuristics that search for a local minimum instead. In these cases, we see that the quality of the results on running the heuristics on the coreset is better than running them on the full large set. Intuitively, the coreset "cleans" the data, removes noise and reduces the number of bad local minima. Hence, in practice the resulting fitting error using the coreset is smaller. For example, by running the common "RANSAC" heuristic on our coreset, our experimental results show that the results are more accurate and faster; See Figure 4 in Section 3.4.
Moreover, while most of the existing solutions aim to detect a circle in two-dimensional space, our algorithm returns provable approximation for every d ≥ 2.
As an example of a robotics application, we demonstrate our algorithms on a mechanical pressure control system which requires an efficient and extremely accurate tracking of the dynamics of a falling ball using a laser profiler, see Section 4.

Overview
In Section 2.1 we formalize the sphere fitting problem. In Section 2.3 we present and describe our shape fitting algorithms. In Section 2.4 we review the theory and the analysis of the algorithms. In Section 3 we show experimental results using our algorithm, and measures its accuracy compared to other algorithms. In Section 4 we present an operating mechanical pressure control system that uses our algorithms. In Section 5 we conclude our work and suggest open problems.

Median-Sphere
The median-sphere problem is to compute a sphere that best fits a given set of input points in the Euclidean d-dimensional space, i.e., a sphere that minimizes the sum of distances to these points, over every unit sphere in R d . For the case of two-dimensional space (i.e., points on the plane), the sphere is a circle of points on the plane.

Notation
For a center (point) x ∈ R d and radius r ≥ 0, we define sphere of radius r around x by We denote by S d−1 the union over every sphere in R d . Let P = {p 1 , p 2 . . . p n } ⊆ R d be a set of n points. For every S ∈ S d−1 , and p ∈ P, we define distance between p and S by dist(p, S) = min q∈S q − p .
For every S ∈ S d−1 and a weight function u : P → [0, ∞), we denote the weighted sum of distances to the sphere by Hence, smaller cost implies a better fit of the sphere to the points. For an unweighed set P, we use a default function u : P → {1} that maps a weight of 1 to every point in P.

Problem Statement
Given a set P of n points in R d , the median-sphere is a sphere that minimizes the sum of distances to these points over every sphere in R d . We denote it by We now define an α-approximation to the median-sphere problem.

Thinnest-Sphere
The thinnest-sphere problem is to compute a sphere that minimizes the maximum of distances to a given set of input points in the Euclidean d-dimensional space, over every unit sphere in R d . Formally, for a given set P of n points in R d , we denote the thinnest-sphere by

Algorithms
In this section we present our main algorithms. In Section 2.3.1 we present a polynomial time algorithm for computing the optimal median-sphere in time n O(d) . In Section 2.3.2 we present our main algorithms for constructing an ε-coreset of P for sphere median. Finally, we apply the exact (inefficient) algorithm on the coreset.

Optimal Solution for the Median-Sphere
In order to solve the median-sphere problem we restate the problem as real polynomial system. Let p ∈ R d , and S ∈ S d−1 be a sphere of radius r > 0 which is centered at x ∈ R d . Hence, dist(p, S) =| p − x − r | is the distance between p and S. That is, dist(p, S) = r − p − x if the point p is inside the sphere, and dist(p, S) = p − x − r if the point p is outside the sphere.
In what follows, sgn(p, S) denotes the sign of dist(p, S). More precisely, sgn(p, S) = 1 if p − x ≥ r, and sgn(p, S) = −1 otherwise. Therefore, Recall that the number that minimizes the sum of distances to a set of numbers is a median of this set. Hence, opt 1 (P) is a sphere that minimizes the sum of distances to the input points, and so must separate P such that the difference between the number of points inside the sphere and the number of points outside the sphere must be less than the number of points that are placed on the sphere. That is, the number of possible permutations for constructing a set of equations (i.e., to determine a possible sign for each point) is limited by the number of spheres that can be defined using a tuple of d + 1 points, i.e., n (d+1) spheres. Hence, we define a solution for the median-sphere problem by finding the sphere that minimizes the sum of distances to each of the n (d+1) polynomial systems, each consisting of n equations. The optimum among these polynomial systems is the S * ∈ S d−1 (that is defined by its center x * and radius r * ) that minimizes the sum of distances ∑ n i=1 dist(p i , S * ). In Algorithm 1, we present such algorithm for the median sphere.
Let P be a set of n points in R d . We denote by P d+1 a tuple of d + 1 different points from P. We denote by SPHERE(P d+1 ) the d-dimensional sphere that defined by P d+1 . We denote by POLYNOMIAL(P, u, S) a system of n equations (polynomial system) that optimizes the median-sphere problem. We denote by SOLVER( f ) an optimization for polynomial system f , e.g., using Maple [31] or Wolfram Alpha [32], etc... Set S * = S 8 Return S * Clearly, such a solution is not applicable to a large set of points. Instead, we offer the following solution in Algorithm 2 (MEDIANSPHERECORESET), which reduces the amount of data and then runs the existing solution on the resulting small set, called a coreset.

Coreset Construction
In order to reduce the number of input points, we suggest Algorithm 2 (called MEDIANSPHERECORESET) to construct a coreset for the median sphere by removing layers of the thinnest-sphere surface using Algorithm 3 (THINNESTSPHERECORESET). We first give a formal definition for such a coresets. Definition 2 (Coreset for median-sphere). Let P ⊆ R d be a finite set of points, and u : P → [0, ∞) be a weight function. For ε ∈ [0, 1), C ⊆ P, and a weight function w : C → (0, ∞), the pair (C, w) is an ε-coreset for the median-sphere of P if for every sphere S ∈ S d−1 | cost(P, u, S) − cost(C, w, S) |≤ ε · cost(P, u, S).
Definition 3 (coreset for thinnest-sphere). Let P ⊆ R d be a finite set of points. A subset C ⊆ P, is an ε-coreset for the thinnest-sphere of P if for every sphere S ∈ S d−1 , For every point p ∈ P, we define the sensitivity [33,34] of p with respect to every sphere S ∈ S d−1 as its relative contribution to the overall cost(P, u, S), Here we assume that the sup is over every S ∈ S d−1 such that ∑ q∈P dist(q, S) > 0.
Definition 4 (Query space [35]). Let P be a finite set, u : P → [0, ∞) be a weight function, and P = (P, u) denote a weighted set. Let Q be a (possibly infinite) set called query set, f : P × Q → [0, ∞) be called a cost function, and loss be a function that assigns a non-negative real number for every real vector. The tuple (P , Q, f , loss) is called a query space. For every q ∈ Q we define the overall fitting error of P to q by f loss (P , q) := loss((u(p) f (p, q)) p∈P ) = loss(u(p 1 ) f (p 1 , q), . . . , u(p n ) f (p n , q)).
The dimension of a query space (P, u, S d−1 , cost) is the VC-dimension of the range space that it induced, as defined below. The classic VC-dimension was defined for sets and subset and here we generalize it to query spaces, following [34].
Definition 5 (Dimension for a query space [34][35][36]). For a set P and set ranges of subsets of P, the VC-dimension of (P, ranges) is the size |C| of the largest subset C ⊆ P such that | {C ∩ range | range ∈ ranges} |= 2 |C| .
It was proven in [35] that an ε-coreset can be obtained for a given problem (e.g., median-sphere), with high probability, by sampling points i.i.d. from the input points with respect to their sensitivity. Then we assign each sample point with a weight that is inverse proportional to its sensitivity. The number of sampled points should be proportional to the VC-dimension of the query space and the sum of sensitivities ∑ p∈P s(p). Formally, Theorem 1. Let ((P, u), Q, f , · 1 ) be a query space of dimension d , and let n be the size of P. Let s : P → (0, 1) be a sensitivity function for every p ∈ P, and let t = ∑ p∈P s(p) be the total sensitivity of the set P. Let ε, δ ∈ (0, 1). Let c > 0 be a universal constant that can be determined from the proof [35]. Let There is an algorithm CORESET(P, u, s, m) in [35] that outputs weighted set (C, w) such that | C |= m, and, with probability at least 1 − δ, (C, w) is an ε-coreset for ((P, u), Q, f , · 1 ).
In Algorithm 2, we present such a coreset construction for the median sphere.

Overview of Algorithm 2
The input for Algorithm 2 is a set P of n points in R d , a weight function u : P → [0, ∞), and an error parameter ε ∈ (0, 1]. The output of the algorithm is an ε-coreset for median-sphere. In Line 3, we compute (iteratively) a coreset Q for thinnest-sphere of P, see Algorithm 3. In Line 5, we assign a sensitivity s(p) for every point p in Q. The sensitivity of each point depends on the iteration number in which it was selected. In Line 7, we update the set P by removing the coreset Q from P. In Lines 8-11, we sample points to obtain coreset. Set P := P \ Q // Remove the surface of the thinnest-sphere from P. 8 while | C |< m do 9 Randomly sample a point p from P into the coreset (for a total of m =

Overview of Algorithm 3
The input for Algorithm 3 is a set P of n points in R d . The output of the algorithm is a coreset of size d + 2 for thinnest-sphere, as in Definition 3. The thinnest-sphere problem depends on the external sphere surface (convex hull) and the internal sphere surface (interior convex hull) that bound the set P. Since the external sphere surface and the internal sphere surface are concentric spheres, it can be proved that the thinnest-sphere defined by exactly d + 2 points where at least one of the points placed on the internal sphere surface and one of the points placed on the external sphere surface. The thinnest-sphere can be computed in O(dn + 1/ε 3d ) time using Theorem 6.8 in [29].

Algorithm 3: THINNESTSPHERECORESET(P)
Input: A set P = {p 1 , p 2 , . . . , p n } ⊆ R d of points. Output: A subset C ⊆ P which is a coreset of P for the thinnest-sphere, see Definition 3.

Median-Sphere Approximation
We propose an approximation algorithm to fit a sphere (i.e., median-sphere) to a set P of n points in R d that is performed in the following two steps: (1) Coreset construction using Algorithm 2 to reduce the original data set to few points for each measure and (2) run the optimal Algorithm 1 on the resolved small coreset, which ensures (1 + ε)-approximation of the optimal solution.

Analysis
In this section, we analyze the algorithms of the previous section in terms of time complexity, coreset's size and the quality of the approximation. The approximation that is states in Theorem 2 that follows is with respect to the optimal solution opt 1 (P, u) such that the calculated cost of our (sphere) solution is compared directly with the calculated cost of the optimal (sphere) solution.

Accuracy:
Our proof uses the technique of [37] to bounds the maximum distance of each point in P to any sphere.
Let c = n/(d + 2) be the total number of iterations of the "while" loop in Lines 2-7. For every integer i ∈ [1, c]. let Q i := Q (Line 3) be the output of THINNESTSPHERECORESET(P) at the i-th iteration. That is, Q i defines concentric external and internal spheres that bounds the remaining set P. Therefore, by the triangle inequality, we can bound the maximum distance from P \Q i to Q i by Therefore, we can bound the sensitivity of every point q ∈ Q i by where the first inequality holds since the right denominator is a subset of the left denominator, and the second inequality holds since the numerator is uphold 1 − ε of the denominator by definition. Hence, the total sensitivity is where the first equality holds by Line 5 of Algorithm 2. The middle inequality refers to the size of Algorithm 3. The right equality holds by a trivial geometric sequence. Substituting in Theorem 1 the input set P, the VC-dimension d = (d + 2), which raised from a trivial classification (cardinality of the largest set of points that the algorithm can shatter), and the query set Q = S d−1 , yield that with probability at least 1 − δ, the pair (C, w) is an ε-coreset of P for median sphere, and therefore, the size of our coreset is O(d 2 log 2 (n) log(1/δ)/ε 2 ).

Hardware
We implemented Algorithms 2 and 3 from the previous section. The implementation was done in Python language [38] using the packages Numpy [39], Scipy [40] and OpenCV [41]. We ran all the tests on a standard Intel 4810MQ, 2.8 GHz CPU laptop. Open code with scripts for reproducing our experiments is provided in [42].

Data Sets
We generated synthetic data sets of dimensions d = 2, 3, 4. The data sets are of different sizes, and with different types of noise characteristics. We also gather some RGB images that captured rings with a flame of fire [12], that simulates noises and abnormalities.

Experiments
We run the above algorithms and compared their performance with and without our coresets for both real and synthetic datasets.
The results show that our algorithm significantly improves both the running time and accuracy of existing heuristics by applying them to the proposed coresets ("Improved RANSAC"). Of course, for the case of computing the optimal solution (via exhaustive search), the result that is extracted from the coreset cannot be better than the one from the original input. However, as expected by our analysis, the approximation (1 + ε) error is small while the running time is shorter by orders of magnitude.

Coreset'S Size
The coreset size returned by Algorithm 2 is determined by the required approximation error ε, and log 2 dependency on the size of the origin set as required in Line 8. Examples are given in Figure 2.  Figure 3 shows that the runtime of our algorithm is far better than other algorithms, such as the "RANSAC" algorithm that considered fast. Furthermore, using our coreset algorithm to select important points significantly improves these algorithms, i.e., "Improved RANSAC". Running time. Experiment on synthetic data sets for several expected value of error factor (ε). The measured running time performed for the optimal solution on coreset vs. "RANSAC" algorithm on the full data set and "RANSAC" algorithm on the coreset ("Improved RANSAC").

Accuracy
Coresets enable complex algorithms to run in a very short time while ensuring quality results. In addition, it can be seen that coresets significantly improve existing algorithms such as "Improved RANSAC" . Figures 4-6 show a comparison for fitting accuracy in 2, 3 and 4 dimensions.  The 3D experiments result. The following four graphs show comparison results between coreset sampling by our algorithm with uniform sampling and the same amount of sample data, compared to the RANSAC algorithm (on the full set and on the coreset) that runs at the same time. Each graph displays results for a differently-sized data set.

Fitting a Sphere over Images
Other visible examples allow us to show how our algorithm is able to identify and validate the "exact" circle within a set of points representing diverse images. Figure 7 shows samples of such fitting over 2D images.

Example System-Mechanical Pressure Control
In this section, we present an experimental mechanical pressure control system that uses our coreset construction to track a falling ball; See Figure 8. The goal of the system is to learn and compare the kinetic effect on different types of materials, namely rubber membranes that consist of an aluminum layer with a thin layer of thick silicone oil. This is by using pressure sensors located below these materials, and by tracking the movement of a ball (marble) that is falling from a filter to a ground (board) that is made by one of these materials.
Comparing the path on different boards and from different heights implies the properties of each of the materials which are useful for many applications in material engineering [43][44][45]. The tracking is done via λ = 200 laser beams that scan the depth in each direction till hitting an object; see Figure 9. In our case, the object is either the ball (in a depth 1/2 that depends on its diameter) or the white sticker behind the falling path of the ball.
The beam is directed toward this path, so in time t 0 , before the ball begins to fall, the depth vector is roughly v 0 = {1, · · · , 1} ∈ R λ . Here, the i-th entry of v 0 denotes the depth that is sensed by the i-th laser beam. This result stays the same while the ball does not enter the scanned frame of the beams, time t 1 , so v 1 = v 0 . When the ball enters the frame of beams, say on t 2 , the first beams hit the ball so v 1 ∼ {0.335, 0.34, 0.35, 0.37, 1, · · · , 1} ∈ R λ . In the middle of the fall, time t 3 , all the ball enters the frame of the scanner. Then, in time t 4 , the ball hits the board, and finally, after hitting the board, time t 5 , the ball begins to jump and oscillate. These sets of depth vectors are based on the motion of the ball which is strongly depends on the material of the board.

Fitting Method
We apply our circle fitting algorithm on the image that corresponds to its vector v, where the x-axis represents the index i which is integer in λ, and the y-axis represents the value v i of the i-th index of v ∈ {1, · · · , λ}. The results for some of the vectors over time are shown in Figure 10, where on the right is an illustration of the system state and on the left a graph corresponding to the vector v.
Each of these thousands of measurements required very high precision of circle fitting. We then run our Algorithm 2 to construct a coreset, thereby significantly reducing the number of points in each data set (from hundreds of points to only twenty points), which allows efficient solution using a polynomial system for those thousands of tests. Time t 2 . The marble partially enters the beamprojecting area. Some beams hit a part of the entered marble and the rest of the beams hit the white sticker.
Time t 3 .
The marble covered by the beam-projecting area. Some beams hit the entered marble and the rest of the beams hit the white sticker.
Time t 4 . The marble hits the board, so the marble is partially covered by the beams.
Time t 5 . The marble oscillates on the board upwards.

Fitting Method Results
In Figure 11, we show some examples of our fitting results. Figure 11. Circle fitting: The following images show examples for fitting circles (red) to a set of points (blue) with high accuracy as required by the system. Although the samples contain noise and distortions, using our system, we were able to fit an accurate circle to all the samples.

Conclusions
We suggested a provable and practical algorithm that fits a sphere for a given set of points in the Euclidean d-dimensional space. The algorithm is based on first computing a coreset that is tailored for this problem, and then running the existing inefficient solution on the small coreset. The result is provable for (1 + ε) approximation in a provably fast and practical running time O(n) for constant d and ε > 0.
When running the optimal (inefficient) optimal solution on our coreset we indeed obtain running time that is smaller by orders of magnitudes (minutes instead of days) in the price of ε ∼ 0.001 approximation error. Finally, we run experiments on a real-world physical system that was the main motivation for writing this paper.
Future research directions include generalization to other shapes, proving lower bounds on the size and computation time of the corresponding coreset, and applying the algorithms on more systems.