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.

Caseation

Probably not NP-complete because it's not a decision problem. Sounds vaguely Knapsack-y though, so I'd wager NP-hard. – Indiana

@Indiana A lot of NP-hard optimization problems (including KP) become NP-complete once coupled with a threshold; in this case, it would be "is there a solution involving n spheres and an error metric no greater than e"? If we take the union-of-convexes approach and force the approximation to be a superset of the original solid, in fact, there's a pretty trivial reduction to 3SAT. – Caseation

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.

Cragsman

Welp, at first glance that paper does seem to address my need. :) I was aware of the connection to the medial axis transform, but I couldn't convince myself that, in the case of a non-infinite number of spheres,it would necessarily be an optimal idea to put them all on medial axis points. Consider a capsule formed by the convex hull of two touching balls; the best way to approximate that with, say, five spheres would be one on each end of the medial axis line segment, and three arranged in some clever manner off-axis between them. I'll read the paper more closely, though. – Caseation

@Sneftel: I think in your capsule example, rather than three off-axis, it would be better to center three on-axis. In any case, the crust algorithm is for arbitrary shapes. Geometrically simple examples like a capsule will often have better ad hoc solutions exploiting their symmetry. – Holpen

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.)

Pandiculation

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.

Phalange

Answer 4

A good start would be to develop an algorithm to:

1) Find the largest radius of an inscribed sphere.

2) Then consider the subtract volume

3) Approximate the subtract volume by a polyhedral.

4) Subdivide that polyhedral into smaller (finer) polyhedrals.

5) Redo step 1.

There might be some variations, but it seems to answer your question. As you can see, the error function decreases as the number of spheres increases, so you can't minimize both: that is a trade-off. But you can fix one and try to optimize another, e.g. fix the error function to be an acceptable tolerance, and minimize the number of spheres to do it, or vice versa.

Marte

That algorithm is reasonable, except that it doesn't allow spheres to overlap each other; that seems like it would severely limit the quality of the solutions produced. – Caseation

There is a closed form expression to calculate the sphere-sphere intersection volume. I guess you can have some tolerance to let your spheres intersecting each other, but that would increase the number of spheres you want to use. So the question is: to what limit you can accept this sphere-sphere intersection? Do you have any specific value of fraction of volume that you would like the spheres to cover the polyhedral? 95%? 99%? – Marte

Erm, perhaps I wasn't clear enough in the original question. There's no problem with spheres intersecting each other; the approximation volume is their union. Obviously one sphere being entirely contained in another would be a waste of a sphere, but in general allowing sphere intersection is crucial for good-quality results, and your solution doesn't seem to support that (because the recursion is only on the result of the subtraction). – Caseation

I guess so. It doesn't give you a detail answer but you can always improve the idea to the stage of sufficiently enough to publish a paper. The problem you are facing arises in some fields, for example, CFD where only sphere primitive is allowed. You would have to dive into the topic to get some results, but you are absolutely right. – Marte

I'd modify this proposed solution to include the constraint that the center of any sphere created after the first step be in one of the segments connecting the center of any previous sphere to the farthest points in the hull. This would give a set of "candidate" spheres for any region of remaining volume to fill; from this set the largest spheres can be selected, keeping the number of spheres down. May or may not lead to something useful.. – Tinware