Trying to find a bounding sphere.. doing something wrong, but can't figure out what

Submitted by Malkalypse on Tue, 2009-06-16 14:37

Okay, I've been competely unable to find a Maxscript for finding an object's minumim enclosing sphere (bounding sphere), so I'm trying to write my own.

The basic concept is:

1) Start with the first few points of an object

2) Make spheres based on their positions (i.e. either from between 2 points, or at the circumcenter of a triangle formed by 3 points).

3) Keep the smallest sphere that still includes all of the checked points, delete the rest.

4) Ignore points that do not reach as far out as that smallest sphere (i.e. ones well within the boundaries), then add the next point to be checked

5) Repeat from step 2

Unfortunately, I seem to be having a hard time putting this into practice. My attempts so far have brought me up to the following. I have some code in there to output debugging data to the listener, and don't understand why the code is allowing certain things to happen.

global RO_KnuckleBall

try(destroyDialog RO_KnuckleBall)catch()

global s

-- special thanks to Joel Hewitt (a.k.a. Gravey)
fn circumcenter p1 p2 p3 =
(
   BC = distance p2 p3
   CA = distance p3 p1
   AB = distance p1 p2

   baryCoords = [ (BC^2 * (CA^2 + AB^2 - BC^2)), (CA^2 * (AB^2 + BC^2 - CA^2)), (AB^2 * (BC^2 + CA^2 - AB^2)) ]
   triArea = baryCoords.x + baryCoords.y + baryCoords.z
   baryCoords /= triArea -- normalize the barycentric coordinates

baryCoords.x * p1 + baryCoords.y * p2 + baryCoords.z * p3
)

-- get measurements and create spheres
fn defineSpheres p =
(
   -- midpoints of lines
   midAB = (p[1] + p[2]) / 2
   try(midAC = (p[1] + p[3]) / 2) catch()
   try(midAD = (p[1] + p[4]) / 2) catch()
   try(midBC = (p[2] + p[4]) / 2) catch()
   try(midBD = (p[2] + p[4]) / 2) catch()
   try(midCD = (p[3] + p[4]) / 2) catch()

   -- midpoints of triangles
   try(centerABC = circumcenter p[1] p[2] p[3]) catch()
   try(centerABD = circumcenter p[1] p[2] p[4]) catch()
   try(centerACD = circumcenter p[1] p[3] p[4]) catch()
   try(centerBCD = circumcenter p[2] p[3] p[4]) catch()

   -- spheres defined by 2 points
   try(s2_AB = geoSphere pos:midAB radius:(distance p[1] midAB)) catch()
   try(s2_AC = geoSphere pos:midAC radius:(distance p[1] midAC)) catch()
   try(s2_AD = geoSphere pos:midAD radius:(distance p[1] midAD)) catch()
   try(s2_BC = geoSphere pos:midBC radius:(distance p[2] midBC)) catch()
   try(s2_BD = geoSphere pos:midBD radius:(distance p[2] midBD)) catch()
   try(s2_CD = geoSphere pos:midCD radius:(distance p[3] midCD)) catch()

   -- spheres defined by 3 points
   try(s3_ABC = geoSphere pos:centerABC radius:(distance p[1] centerABC)) catch()
   try(s3_ABD = geoSphere pos:centerABD radius:(distance p[1] centerABD)) catch()
   try(s3_ACD = geoSphere pos:centerACD radius:(distance p[1] centerACD)) catch()
   try(s3_BCD = geoSphere pos:centerBCD radius:(distance p[1] centerBCD)) catch()
)

-- get position and radius of spheres
fn getSphereData =
(
   s = #(#(), #())
   for o in objects where classOf o == geoSphere and o != undefined do
   (
       append s[1] o.pos
       append s[2] o.radius
   )
)

fn knuckleBall_Start tgt =
(
-- define initial points
   p = #()
   for i = 1 to 4 do p[i] = polyOp.getVert tgt i
   defineSpheres p
   getSphereData()

   -- delete spheres with points that fall outside
   for a = 1 to s[1].count do
   (
       format "looking at sphere #%...\n" a

       for b = 1 to p.count do
       (
           if distance s[1][a] p[b] - s[2][a] > .0001 do
           (
               format " point #% (%) falls outside sphere #%\n" b p[b] a

               isDeleted = false
               for o = 1 to objects.count where classOf objects[o] == geoSphere while isDeleted == false do
               (
                   if objects[o].pos == s[1][a] and objects[o].radius == s[2][a] do
                   (
                       format "   deleting object % (matches sphere #%)\n\n" objects[o] a
                       delete objects[o]
                       isDeleted = true
                   )
               ) -- end o loop
           ) -- end if (distance)
       ) -- end b loop
   ) -- end a loop

   -- delete all but the smallest bounding sphere
   getSphereData()
   for o = 1 to objects.count where classOf objects[o] == geoSphere do
   (
       if objects[o].radius != aMin s[2] do delete objects[o]
   )
   getSphereData()

   -- remove redundant points from the array
   for a = 1 to s[1].count do
   (
       for b = 1 to p.count do
       (
           try
           (
               if distance s[1][a] p[b] < s[2][a] do deleteItem p b
           )
           catch(format "NOTE: array item % was NOT deleted\n" b[p])
       )
   )
)

-- continue creating and checking spheres using the remainder of the object's points.
fn knuckleBall_Cont theTarget i =
(
tempPoint = point size:0.1 pos:(polyOp.getVert theTarget i)

   append p (polyOp.getVert theTarget i)
   defineSpheres p
   getSphereData()

   -- delete spheres with points that fall outside
   for a = 1 to s[1].count do -- for all spheres
   (
       format "looking at sphere #%...\n" a

       for b = 1 to p.count do
       (
           if distance s[1][a] p[b] - s[2][a] > .0001
               then
               (
                   format " point #% (%) falls outside sphere\n" b p[b]

                   for o = 1 to objects.count where classOf objects[o] == geoSphere do
                   (
                       if objects[o].pos == s[1][a] do
                       (
                           format " (deleting sphere #%)\n" a
                           delete objects[o]
                       )
                   ) -- end o loop
               ) -- end if/then/else (distance)

       ) -- end b loop
   ) -- end a loop

   -- delete all but the smallest bounding sphere
   getSphereData()
   for o = 1 to objects.count where classOf objects[o] == geoSphere do
   (
       if objects[o].radius != aMin s[2] do delete objects[o]
   )
   getSphereData()

   -- remove redundant points from the array
   for a = 1 to s[1].count do
   (
       for b = 1 to p.count do
       (
           try
           (
               if distance s[1][a] p[b] < s[2][a] do deleteItem p b
           )
           catch (format "item was not deleted\n")
       ) -- end b loop
   ) -- end a loop

) -- end fn knuckleBall_Cont

rollout RO_knuckleBall "Knuckle Ball"
(

pickbutton pck_Tgt "Pick target" width:100
button btn_Cont "Continue" width:100 enabled:false

   global theTarget
   global tempPoint

   local theLoop = 1
   local p

   on pck_Tgt picked obj do
   (
       theTarget = obj
       knuckleBall_Start theTarget
       btn_Cont.enabled = true
   )

   on btn_Cont pressed do
   (
       if tempPoint != undefined and isValidNode tempPoint do
       (
           tempPoint.size = 0.01
       )

       knuckleBall_Cont theTarget theLoop
       theLoop = theLoop + 1
   )
)
createDialog RO_KnuckleBall

General Scripting

Comments

Submitted by AndyOakley on Thu, 2009-07-02 07:48.

This is one of those great

This is one of those great problems that caught my eye and I`ve been thinking about it since I read it.

I`m working on solving it at the moment. Almost there and it's pretty quick too.

Submitted by Anubis on Fri, 2009-06-26 05:29.

Thank you for your comment.

Thank you for your comment. Yes, you are right. Of course, if the problem was easy it would've been solved long ago. But calculation time is also a factor. My version is less accurate but I think it's more practical for dynamic calculations. And... Why not share it?

my recent MAXScripts RSS (archive here)

Submitted by real08121985 on Fri, 2009-06-26 01:36.

On the left side there is

On the left side there is the bsphere of Malkalypse, on the right side the one of Anubis. You see the problem, the one on the right isn't the minimal enclosing sphere. Its a really difficult problem. Malkalypse script is heavy and needs some time to compute, but it returns nice results :)

Submitted by Anubis on Thu, 2009-06-25 16:39.

Hi Malkalypse, I might write

Hi Malkalypse, I might write a little late. You already done by your own way in which I had no time to think, but that it's an interesting problem and now I sat down to think about.
So here my version:

obj = selection[1] -- select EPoly or add EPoly_mod
allVerts = for i in 1 to obj.numVerts collect polyop.getVert obj i
vDist = for i in allVerts collect distance i obj.center
sort vDist
-- Optional step:
-- if obj is not EPoly and you wish to restore it then remove EPloly_mod
bs_radius = vDist[vDist.count] -- BSphere Radius
bs_pos = obj.center -- BSphere Center

my recent MAXScripts RSS (archive here)

Submitted by real08121985 on Wed, 2009-06-17 16:16.

Very nice code, thank you

Very nice code, thank you for this :)

Submitted by Malkalypse on Wed, 2009-06-17 15:38.

(No subject)

Submitted by Malkalypse on Wed, 2009-06-17 15:22.

Okay, I got a script that

Okay, I got a script that works pretty well, certainly good enough for my purposes:

-- special thanks to Joel Hewitt (a.k.a. Gravey) for circumcenter function
fn circumcenter p1 p2 p3 =
(
    BC = distance p2 p3
    CA = distance p3 p1
    AB = distance p1 p2

    baryCoords = [ (BC^2 * (CA^2 + AB^2 - BC^2)), (CA^2 * (AB^2 + BC^2 - CA^2)), (AB^2 * (BC^2 + CA^2 - AB^2)) ]
    triArea = baryCoords.x + baryCoords.y + baryCoords.z
    baryCoords /= triArea -- normalize the barycentric coordinates

baryCoords.x * p1 + baryCoords.y * p2 + baryCoords.z * p3
)

struct pointData (id, pos) -- vertex number and position
struct pointPair (pointA, pointB, dist) -- 2 points and their distance from each other

allPoints = #() -- array for all object vertices
for i = 1 to $.numVerts do allPoints[i] = pointData id:i pos:(polyOp.getVert $ i)

-- find the two points that are furthest from each other
farPoints = pointPair dist:0
for i = 1 to allPoints.count do
(
    a = allPoints[i]

    for j = 1 to $.numVerts do
    (
        b = allPoints[j]

        if distance (a.pos) (b.pos) > farPoints.dist do
        (
            farPoints.pointA = a
            farPoints.pointB = b
            farPoints.dist = distance (a.pos) (b.pos)
        ) -- end if (distance (a.pos) (b.pos))
    ) -- end j loop
) -- end i loop

-- starting position and radius
sPos = ((farPoints.pointA.pos + farPoints.pointB.pos) / 2)
sRadius = ((distance farPoints.pointA.pos farPoints.pointB.pos) / 2)

theCenter = sPos
theRadius = sRadius

-- find vertices that fall outside the initial hypothetical sphere
outPoints = #()
fn getOutPoints =
(
    for i = 1 to allPoints.count do
    (
        if distance allPoints[i].pos theCenter > theRadius do append outPoints allPoints[i]
    )
)
getOutPoints()

--     For every set of 3 points,
--        if the distance between each of the 3 points is larger than the radial hypotenuse of theRadius,
--        find the circumcenter and radius for the 3 points.
--    If the radius of the 3 points is larger than the existing value of theRadius,
--        update the value of theRadius and theCenter
for i = 1 to outPoints.count do
(
    a = outPoints[i].pos

    for j = 1 to allPoints.count do
    (
        b = allPoints[j].pos

        if distance b a > (theRadius * sqrt 2) do
        (
            for k = 1 to allPoints.count do
            (
                c = allPoints[k].pos

                if distance c a > (theRadius * sqrt 2) and distance c b > (theRadius * sqrt 2) do
                (
                    circ = circumcenter a b c
                    r = distance circ a
                    if r > theRadius do
                    (
                        theRadius = r
                        theCenter = circ
                    ) -- end if (r > theRadius)
                ) -- end if (distance c)
            ) -- end k loop
        ) -- end if (distance b)
    ) --end j loop
) -- end i loop

Submitted by real08121985 on Wed, 2009-06-17 12:26.

Hey nice tip :) thank you

Hey nice tip :) thank you very much

Submitted by real08121985 on Wed, 2009-06-17 11:07.

Internettechnology and

Internettechnology and creating websites .. but thats only my studies, not my future. Algorithms and complex problems can be a lot of fun and maxscript is very handy :D

Submitted by Malkalypse on Wed, 2009-06-17 11:36.

If you're ever looking to

If you're ever looking to make some extra cash on the side, there is a shortage of MaxScript programmers at rentacoder.com. Out of a total of 250,000 coders, only 14 have MaxScript listed in their information pages, and I've only been able to contact 3 of them.

(I should mention I don't work for the site and I'm not trying to advertise for them or anything. There don't seem to be any specific forum guidelines against such things, but I generally find advertising in forums to be poor taste.)

About this site