Approximation of a solid by a union of spheres

Question

I've got a 3D solid, represented as the union of a set of polyhedral convex hulls. (Or a single convex, if that makes things easier.) I'd like to approximate that solid as the union of a set of spheres, in a way which minimizes both the number of spheres in the set and the error in the approximation. (The latter objective is deliberately vague: any reasonable error metric will do. Likewise, the way in which the objectives are combined is up in the air; either the number of spheres or the error metric could be constrained, or some function of the two could be minimized. I don't want to specify myself into a corner.)

The approximation does not need to entirely contain or be entirely contained by the original set. Each sphere may have an arbitrary radius.

This feels like the sort of problem that's NP-complete, and in any case unlikely to be practical using exact methods, so I'm assuming the solution lies in the realm of stochastic optimization. It feels like some variant of k-means might fit (assigning uncovered locations to their closest spheres, and refining the spheres to cover some of them), but I'm not sure how to handle multiply-covered locations, or how to find the local, not-necessarily-covering-everything optimum even for a single sphere. Also, for iterative methods efficiency is important, and doing 3D boolean operations is not going to be efficient.

Answer 1

The problem is not simple, but has been studied previously. The central concept is the medial axis, which can be viewed as a representation of an object by an infinite union of balls. A key paper addressing your question is:

"The power crust, unions of balls, and the medial axis transform." Nina Amenta, Sunghee Choi, Ravi Krishna Kolluri. Computational Geometry. Volume 19, Issues 2–3, July 2001, Pages 127–153. (Journal link.)

---

 
(Images source: From Point Clouds to Power Crusts.)

---

A second paper is

Cazals, Frédéric, et al. "Greedy Geometric Algorithms for Collection of Balls, with Applications to Geometric Approximation and Molecular Coarse‐Graining." Computer Graphics Forum. Vol. 33. No. 6. 2014. (PDF download.)

whose 1st sentence is "Choosing balls which best approximate a 3D object is a non-trivial problem."! Their primary application is to molecular models, which might be far from your interests.

Answer 2

Hm, my best idea so far involves support vector machines. Turn your object into a whole bunch of (probably evenly spaced) points within and on the surface of the object. Train an SVDD model using a linear kernel (see libsvm for an SVDD implementation). The decision function of the model then represents an implicit surface defined by the support vectors of the model (and rho). Turning down the cost will get you more support vectors, turning it up gets you fewer.

Unfortunately, the nature of SVMs is such that the area covered by nearby support vectors will, uh, 'blob' together, sort of like this:

(sorry, my intuition for SVMs is entirely geometric/visual.)

Now, you don't have nice crisp spheres, but (massive hand waving!) hopefully the algorithm chose a useful distribution of centers for spheres.

Finally, you can concoct a function that computes error as a function of radii for spheres centers on all those points. Then just feed that into a nonlinear optimizer and tell it to minimize. Bam.

If you're willing to throw more CPU power at it, you could run another layer of error minimization over top of that one, which reruns the entire above process for different support vector costs, attempting to minimize some combination of error and cost. (Perhaps error/cost.)

Answer 3

This is what I came up with. This approach is more of an iterative 3D boolean operation so it might not be what you're looking for. The surface seems more difficult so I concentrated on that.

Overview

Basically add spheres inside the shape in positions that maximize the coverage of the surface. We convert the sphere into a 3D array of signed byte values. These values are points and will be gobbled up with spheres. We add one spheres at a time inside the object and then grow/shrink it in different directions to "eat" as many points as possible. The goal is to rack up as many points as possible per sphere. Points are earned by summing the points in the area of the sphere. With the addition of each sphere we count then count that area as used and set the Array values to 0.

  (A)       (B) ZZZZZZ   (C) ZZZZZZ   (D) ZZZZZZ   (E) ZZZZZZ   (F) ZZZZZZ   
     /\         ZX33XZ       ZX33XZ       ZX33XZ       ZX33XZ       ZX33XZ   
    /  \       ZX3223XZ     ZX3223XZ     ZX##23XZ     ZX  ##XZ     ZX    XZ  
   /    \     ZX321123XZ   ZX321123XZ   ZX####23XZ   ZX  ####XZ   ZX      XZ 
  |      |   ZX32111123XZ ZX32111123XZ ZX######23XZ ZX  ######XZ ZX        XZ
  |      |   ZX32111133XZ ZX32111133XZ ZX######23XZ ZX  ######XZ ZX        XZ
  |      |   ZX32222223XZ ZX##222223XZ ZX3####223XZ ZX3  ####3XZ ZX3     ##XZ
  |------|   ZX33333333XZ ZX##333333XZ ZX33##3333XZ ZX33  ##33XZ ZX33    ##XZ
   X= -1     ZXXXXXXXXXXZ ZXXXXXXXXXXZ ZXXXXXXXXXXZ ZXXXXXXXXXXZ ZXXXXXXXXXXZ
   Y= -2     ZZZZZZZZZZZZ ZZZZZZZZZZZZ ZZZZZZZZZZZZ ZZZZZZZZZZZZ ZZZZZZZZZZZZ

(A) The shape that we want to fill. (2D used here but 3D would be similar)

(B) Convert the shape in a 3D grid of points. Surface gets largest number and as you work to the center the numbers settle on low positive numbers(like 1); Outside the shape gets negative numbers; deep interior gets 1

(C) Add a small sphere. (we will grow it)

(D) Expand the sphere until we gobble up the maximum number of points

(E) Add the Next sphere and grow it.

(F) Add another sphere. This one is small.

Process

5) first break down the shape into a 3D block resolution.

10) Then give the most "points" to the blocks around the surface. High points with the block that actually touches the surface and lower points as you move inward or outward. As you go outward the points should quickly become negative as this will prevent protruding spheres. As you move inward from the surface the points should settle at 1 so that the central area would be all ones. These points can be setup in different ways to produce different results.

15) Find a location on the inside edge of the shape that is outside a sphere. While being near the edge is not required it does minimize the search area. If a location cannot be found goto 80. If a location cannot be found that is near

20) Draw a sphere with a zero radius inside the shape that is not covered

25) Move the sphere up/down until the points are maximized (simulated annealing)

26) Move the sphere in/out until the points are maximized

27) Move the sphere left/right until the points are maximized

28) Move the Top of the sphere up/down until the points are maximized (simulated annealing/hill climbing)

29) Move the Bottom the sphere up/down until the points are maximized

30) Move the Left side of the sphere in/out until the points are maximized

31) Move the Right side of the sphere in/out until the points are maximized

32) Move the Front of the sphere in/out until the points are maximized

50) If any changes in 25-32 then repeat (goto 25)

70) Subtract out the points from the last added sphere. Set all values to zero for internal values(positive numbers) and do not adjust external values(negative numbers). Goto 15.

80) (optional) Fill in an internal gaps. Basically visit each element to make sure it is less then or equal to 0. If positive then goto 20 with that element. Else, if none found then goto 85. Note:all outside values should be negative and all internal values that are in a sphere should be 0.

85) Finished

Notes

• Since there would a grid of values, a 1000x1000x1000 workspace would consume up 1GB.
• Not an exact solution
• Could take lots of compute for higher resolutions.
• For efficiency, smaller spheres can have their pixel ranges pre-recorded so that the sphere does not need to be calculated for each iteration.
• A lower resolution(i.e. 500x500x500) version could be completed first and then the location and size of the spheres could be applied to a larger 1000x1000x1000. Also for very large shapes sub-sections could be solved.