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