As a part of a dynamics simulator I'm attempting to write, I've created 
a 'real' object center finder.
I know of rune's particle system and the sphere+rod dynamics mechanics, 
but my main goal here is the ability to handle concave objects.

This macro takes the object as the parameter and does a 
trapezoidal-approximation integral of the volume. The tradeoff between 
accuracy and speed is controlled by the 'Divisions_X' and 'Divisions_Y' 
values.
The approach taken was:
1) Find the bounding box of the object
2) Take the x-y face of the bounding box and subdivide, shooting trace 
'rays' in the z direction from the center of each subdivision
4) 'trace' is called again and again from each subsequent hit, ensuring 
that any amount of voids or concavities are found
5) Find the center of each section of the object found, weight by the 
volume of the section, and add to the total volume
6) Divide final center by total volume to find proper center

This is exact for normal rectangular prisms. For curved shapes the 
volume is within 1% accuracy when the subdivision matrix is 50x50, all 
while taking less than a second. The center of mass is even more 
accurate than the volume, in my tests it shows convergence within 1% at 
matrices of size 20x20 or smaller.

If anyone has a better formulation for this macro, please let me know. 
Otherwise, I hope this is possibly useful for someone.

----------------------------
// FindObjectCenter(Object_Identifier) will return a calculated 'real' 
center of mass location
// This assumes a constant density. The calculated volume is also 
available but currently not returned.
// 2009.09.04 Christopher Shake (cshake+pov### [at] gmailcom)
#macro FindObjectCenter(Object_Identifier)
   #local maxcorner = max_extent(Object_Identifier);
   #local mincorner = min_extent(Object_Identifier);
   #local max_x = max(maxcorner.x,mincorner.x);
   #local min_x = min(maxcorner.x,mincorner.x);
   #local max_y = max(maxcorner.y,mincorner.y);
   #local min_y = min(maxcorner.y,mincorner.y);
   #local max_z = max(maxcorner.z,mincorner.z);
   #local min_z = min(maxcorner.z,mincorner.z);
   #local bound_center = (mincorner+maxcorner)/2; // center of bounding 
box, as reference
   // calculate CoM
   #local CoM = <0,0,0>;
   #local dy = abs(max_y-min_y)/Divisions_Y;
   #local dx = abs(max_x-min_x)/Divisions_X;
   #local temp_norm = <0,0,0>;
   #local total_volume = 0;
   #local curr_y = min_y+dy/2;
   #while (curr_y < max_y)
     #local curr_x = min_x+dx/2;
     #while (curr_x < max_x)
       #local more_sections = true; // deals with voids in the object
       #local z_shoot_start = <curr_x,curr_y,min_z-max(dy,dx)>;
       #while(more_sections & z_shoot_start.z <= max_z)
         // shoot ray along this '1-d' 'slice' to find ends
         #local z_hit_min = 
trace(Object_Identifier,z_shoot_start,z,temp_norm);
         #if(vlength(temp_norm)!=0)
           // object exists in this 'slice', z_hit_min is now the first 
'outside' edge
           // shoot another to see the next edge
           #local z_hit_max = 
trace(Object_Identifier,z_hit_min,z,temp_norm);
           #if(vlength(temp_norm)!=0)
             #local point_vol = (z_hit_max.z-z_hit_min.z)*dx*dy;
             #local curr_ctr = (z_hit_max+z_hit_min)/2;
             #local CoM = CoM + curr_ctr*point_vol;
             #local total_volume = total_volume + point_vol;
             // now reset the starting position of the ray shoot in case 
the object has voids
             #local z_shoot_start = z_hit_max;
           #else
             #local z_shoot_start = z_hit_min+z;
           #end
         #else
           #local more_sections = false;
         #end
       #end
       #local curr_x = curr_x + dx;
     #end
     #local curr_y = curr_y + dy;
   #end
   #local CoM = CoM / total_volume;

   // 'CoM' contains the calculated center of mass
   // 'total_volume' contains the calculated volume of the object
   (CoM)
#end

Post a reply to this message