I have a mesh, I want to find the smallest possible bounding sphere that
encompasses all the points. I know this is a colossal combinatorial problem
and to test every possible bounding sphere were talking about an algorithm
that runs in O(k*(k-1)^2*(k-2)^3*(k-3)^4) time, which is completely
unreasonable even for small meshes. I did some reading into the subject, but
all I could find was a bit of really poorly documented code and a couple of
brief mentions of Welzl's algorithm for computing bounding volumes. Anyone
have any thoughts on a method of finding this? so far I have two possible
solutions, but neither of them are super accurate. I've been pondering
heuristic methods, I know how they work but I really and truly do not have
any idea as to how I might go about implementing them.

Solution 1:

#declare Centre = <0, 0, 0>;
#declare i = 0;
#while(i < Pts)
    #declare Centre = Centre + Point[i];
    #declare i = i + 1;
#end

#declare Centre = Centre / Pts;

#declare Radius = 0;
#declare i = 0;
#while(i < Pts)
    #if(vlength(Points[i] - Centre) > r)
        #declare Radius = vlength(Points[i] - Centre);
    #end
    #declare i = i + 1;
#end

Solution 2:

#declare Centre = (min_extent + max_extent) / 2;

#declare Radius = 0;
#declare i = 0;
#while(i < Pts)
    #if(vlength(Points[i] - Centre) > r)
        #declare Radius = vlength(Points[i] - Centre);
    #end
    #declare i = i + 1;
#end

--
Kevin
http://www.geocities.com/qsquared_1999/
#macro _(r)#if(r<12)#local i=asc(substr("oqshilacefg",r,1))-97;
disc{<mod(i,7)-3,div(i,7)-1,6>,z,.4pigment{rgb 10}}_(r+1)
#end#end _(1)//KL

Post a reply to this message

> I have a mesh, I want to find the smallest possible bounding sphere that
> encompasses all the points.

This is a really interesting problem.

After thinking about it a bit, the only solution I can come up with is that
if you find the two points which are farthest apart, and then center the
sphere right between them, with a radius of (their distance
apart)/2*sqrt(3), all of the points will be in that sphere. (The sqrt(3) is
necessary because there could be a third point forming an equilateral
triangle with the first two.)

However, I'm tempted to think that there's a possible solution that will
find it perfectly. Maybe something like starting with a tiny (r=0) sphere in
the midst of the points, and growing its radius and moving it around to make
it encompass one point after another... I dunno.

 - Slime
[ http://www.slimeland.com/ ]

Post a reply to this message

Wasn't it Kevin Loney who wrote:
>I have a mesh, I want to find the smallest possible bounding sphere that
>encompasses all the points. I know this is a colossal combinatorial problem
>and to test every possible bounding sphere were talking about an algorithm
>that runs in O(k*(k-1)^2*(k-2)^3*(k-3)^4) time, which is completely
>unreasonable even for small meshes. I did some reading into the subject, but
>all I could find was a bit of really poorly documented code and a couple of
>brief mentions of Welzl's algorithm for computing bounding volumes. Anyone
>have any thoughts on a method of finding this? so far I have two possible
>solutions, but neither of them are super accurate. I've been pondering
>heuristic methods, I know how they work but I really and truly do not have
>any idea as to how I might go about implementing them.
>
>Solution 1:
>
>#declare Centre = <0, 0, 0>;
>#declare i = 0;
>#while(i < Pts)
>    #declare Centre = Centre + Point[i];
>    #declare i = i + 1;
>#end
>
>#declare Centre = Centre / Pts;
>
>#declare Radius = 0;
>#declare i = 0;
>#while(i < Pts)
>    #if(vlength(Points[i] - Centre) > r)
>        #declare Radius = vlength(Points[i] - Centre);
>    #end
>    #declare i = i + 1;
>#end
>
>Solution 2:
>
>#declare Centre = (min_extent + max_extent) / 2;
>
>#declare Radius = 0;
>#declare i = 0;
>#while(i < Pts)
>    #if(vlength(Points[i] - Centre) > r)
>        #declare Radius = vlength(Points[i] - Centre);
>    #end
>    #declare i = i + 1;
>#end


If you don't need incredible accuracy, then a stochastic approach might
be reasonable.

I started with your Solution 2, having observed that your Solution 1
would be likely to give poor results for meshes that have the points
concentrated towards one end - like the heads of Poser figures.

I turned the Radius calculation into a macro

#macro Rad(C)
 #local i=0;
 #local Radius=0;
 #while(i < Pts)
  #if(vlength(Points[i] - C) > Radius)
        #local Radius = vlength(Points[i] - C);
    #end
    #local i = i + 1;
 #end
 Radius
#end

I displaced the centre in a random direction and recalculated the
Radius. If the new radius was smaller than the previous one, I used the
new centre as my new starting point and displaced again. I repeated the
process, gradually reducing the size of the random displacement.

#declare K=seed(123);
#declare D=vlength(min_extent(Fred)-Centre)/2;
#while (D>0.001)
  #declare Centre2=Centre + <(rand(K)-0.5)*D, (rand(K)-0.5)*D, 
        (rand(K)-0.5)*D>;
  #declare Radius2 = Rad(Centre2);  
  #if (Radius2 < Radius)
    #declare Radius=Radius2;
    #declare Centre=Centre2;
  #end
  #declare D=D*0.95;
#end


Some experiments showed that it was not uncommon to have a long unlucky
streak of random choices and end up with poor results, so I tried also
stepping in the opposite direction, and this seemed to give a
considerable improvement.

#declare K=seed(123);
#declare D=vlength(min_extent(Fred)-Centre)/2;
#while (D>0.001)
  #declare Delta = <(rand(K)-0.5)*D, (rand(K)-0.5)*D, (rand(K)-0.5)*D>;
  #declare Centre2=Centre + Delta;
  #declare Radius2 = Rad(Centre2);  
  #if (Radius2 < Radius)
    #declare Radius=Radius2;
    #declare Centre=Centre2;
    #debug concat("D ",str(D,1,3)," Radius ",str(Radius,1,10),"\n")
  #end
  #declare Centre2=Centre - Delta;
  #declare Radius2 = Rad(Centre2);  
  #if (Radius2 < Radius)
    #declare Radius=Radius2;
    #declare Centre=Centre2;
    #debug concat("D ",str(D,1,3)," Radius ",str(Radius,1,10),"\n")
  #end
  #declare D=D*0.95;
#end

When Pts is large, you may want to change that "#declare D=D*0.95;" to
something like "#declare D=D*0.90;".

When Pts is large, you might consider removing all points that are
closer to the original Centre than Radius/Sqrt(3), since these interior
points can't possibly affect the bounding sphere.

-- 
Mike Williams
Gentleman of Leisure

Post a reply to this message

> Some experiments showed that it was not uncommon to have a long unlucky
> streak of random choices and end up with poor results, so I tried also
> stepping in the opposite direction, and this seemed to give a
> considerable improvement.

Perhaps, if stepping in a certain direction is successful (it reduces the
radius), you should make the next step in the same direction.

 - Slime
[ http://www.slimeland.com/ ]

Post a reply to this message

On Sat, 28 Dec 2002 21:17:11 -0700, "Kevin Loney" <klo### [at] pt2mcom>
wrote:

>I have a mesh, I want to find the smallest possible bounding sphere that
>encompasses all the points. 

Check the Graphic Gems collection, there are various algorithms
implemented in Pascal, Assembler and C (and maybe other languages) for
finding the smallest sphere containing n points.


Peter Popov ICQ : 15002700
Personal e-mail : pet### [at] vipbg
TAG      e-mail : pet### [at] tagpovrayorg

Post a reply to this message

Wasn't it Slime who wrote:
>> Some experiments showed that it was not uncommon to have a long unlucky
>> streak of random choices and end up with poor results, so I tried also
>> stepping in the opposite direction, and this seemed to give a
>> considerable improvement.
>
>Perhaps, if stepping in a certain direction is successful (it reduces the
>radius), you should make the next step in the same direction.

Oddly, experiments show that tends to make it worse.

-- 
Mike Williams
Gentleman of Leisure

Post a reply to this message

Kevin Loney wrote:

> subject, but all I could find was a bit of really poorly documented code
> and a couple of brief mentions of Welzl's algorithm for computing bounding
> volumes. Anyone have any thoughts on a method of finding this? so far I

Hey! That Welzl is a professor at my university! I even had a lecture with 
him... So I did a bit research about this algorithm you mentioned. 
Following the link on Welzl's page (http://www.inf.ethz.ch/personal/emo/) I 
found this: "Smallest enclosing disks (balls and ellipsoids)":
 http://www.inf.ethz.ch/personal/emo/ps-files/SmallEnclDisks-LNCS555.ps

Haven't read through it but maybe it's useful for you

- Micha

-- 
objects.povworld.org - The POV-Ray Objects Collection
book.povworld.org    - The POV-Ray Book Project

Post a reply to this message

I don't know if this is a good approach, but I think it can be very fast:

Find the bounding box, which is "cheap", and now find the smallest sphere
that includes this box.

I think it really is up to your code to decide if it is worth to use a more
exhaustive approach (and more expensive), if that would speed up other parts
of your code.

I hope this helps,

Fernando.

Post a reply to this message

I'm not sure if that will lead to the smallest sphere possible, because the
rotation of the bounding box will matter a lot. If the rotation is wrong,
the bounding box will be too big and so will be the sphere.

Post a reply to this message

In article <3e11c6bb$1@news.povray.org>,
 "Apache" <apa### [at] yahoocom> wrote:

> I'm not sure if that will lead to the smallest sphere possible, because the
> rotation of the bounding box will matter a lot. If the rotation is wrong,
> the bounding box will be too big and so will be the sphere.

It isn't even close, the sphere will be very "loose" in almost every 
case. It will only be minimal if the object touches the bounding box at 
two opposite corners and the box is a perfect cube. It might work well 
as a first approximation, though. I haven't done any research on it, but 
here's what I'd try:

Calculate the sphere enclosing the box.

Compute the minimum radius required to contain all points.

There will now be at least one point on the sphere. Go to the step with 
the number of contacting points:

1: Move the sphere toward the contact point, reducing radius to keep it 
on the sphere surface, until one other point is on the surface.

2: The same as step 1, but move toward the midpoint of the two existing 
points.

3: The 3 points define a triangle. Do the same thing as the first two 
steps, but toward a point in the triangle that is equidistant from all 3 
points.

4: Stop. The sphere is pretty tight, maybe the best possible.

There are probably mathematical solutions for all of these steps, but it 
might be easier to do it iteratively, moving the sphere in small steps 
or doing a binary search.

-- 
Christopher James Huff <cja### [at] earthlinknet>
http://home.earthlink.net/~cjameshuff/
POV-Ray TAG: chr### [at] tagpovrayorg
http://tag.povray.org/

Post a reply to this message