 |
|
|
|
|
|
| |
| |
|
|
|
|
| |
| |
|
|
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] gmail com)
#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
|
|
| |
| |
|
|
|
|
| |
| |
|
|
CShake wrote:
> 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.
Could you add it to the Object Collection please? It might be useful for
people in the future.
-Mike
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
SharkD wrote:
> Could you add it to the Object Collection please? It might be useful for
> people in the future.
>
> -Mike
I plan on putting the completed collection, but wanted some feedback to
see if it worked well enough before adding it. I think everything except
for the function name itself complies with the collection requirements.
I figure that most people don't have a use for center of mass right now,
except possibly trying to get an object to rotate around the center, but
by changing the output variable on this macro I can see the 'total
volume' being useful to some.
My next macro will be finding the moment of inertia of an object, so I
will probably package these together since they're similar and upload it
as the first group.
If anyone wants to look at the actual math here for any obvious things
I've missed I'd appreciate it too.
Chris
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
CShake wrote:
>...
> I figure that most people don't have a use for center of mass right now,
> except possibly trying to get an object to rotate around the center, but
> by changing the output variable on this macro I can see the 'total
> volume' being useful to some.
>
> My next macro will be finding the moment of inertia of an object, so I
> will probably package these together since they're similar and upload it
> as the first group.
>
> If anyone wants to look at the actual math here for any obvious things
> I've missed I'd appreciate it too.
>
Very nice, Chris! To the contrary, I think there are at least a few people who
can use this, I being one of them.
The only issue I see immediately is with very complicated objects such as those
based on fractals, noise, or some odd mathematical formula. Some fractal
objects, for example, may have an infinite number of intersections along a test
ray and should, at the very least, increase the calculation time a great deal.
I hope you don't mind if I make a companion macro based on your code to
calculate the CoM and volume using a more naive approach specifically for such
pathological objects. If I can get it to work, the approach I'm thinking of
would be less accurate for a given number of subdivisions, but would always
finish in n^3 iterations.
Would you mind posting your test scenes, or include some when you post this to
the object database submission? It would be helpful for comparison.
Thanks for posting this.
~David
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
waggy wrote:
> The only issue I see immediately is with very complicated objects such as those
> based on fractals, noise, or some odd mathematical formula. Some fractal
> objects, for example, may have an infinite number of intersections along a test
> ray and should, at the very least, increase the calculation time a great deal.
True. I picked the 'orientation' of what direction to shoot the rays as
somewhat arbitrary. For most objects it should work fine, but I am aware
specifically that an image-based heightfield (since it normally rises up
from the x-z plane) will create a large number of intersections, and
would be more efficient to shoot along y. Considering the likelihood of
this sort of comparison, I may actually switch it to shooting along y,
since for 'normal' object it wouldn't make much of a difference.
> I hope you don't mind if I make a companion macro based on your code to
> calculate the CoM and volume using a more naive approach specifically for such
> pathological objects. If I can get it to work, the approach I'm thinking of
> would be less accurate for a given number of subdivisions, but would always
> finish in n^3 iterations.
I did have in mind a few different ways of shooting rays, specifically
doing a sort of spherical shoot where theta and phi are the control
variables, but evenly distributing the rays around the surface (and
finding the volume of each 'chunk') would be a little more complicated.
At least this method works in k*n^2, where k is something like the
number of voids in the object, which is slightly more optimized than n^3
in the case of most simpler objects (CSG, etc). I thought about trying
to use an adaptive sampling method, but the extra code would be
significant (and I don't have a good handle on how to do it yet, so it
would involve more research too).
It could be helpful to do a sort of basic test if the object is a
primitive, and then use the appropriate mathematical formula, but the
amount of time you'd run into that 'in the wild' would be nearly 0.
> Would you mind posting your test scenes, or include some when you post this to
> the object database submission? It would be helpful for comparison.
#declare BoxObj = box{-1,1}
#declare SphereObj = sphere{0,1}
#declare CSGObj = difference{object{BoxObj} object{SphereObj}}
#declare TorusObj = torus{2,0.2}
#declare Torus2 = union{object{TorusObj} object{TorusObj rotate 90*x
translate 1*x}}
Run the macro for each, send the calculated values to debug, and compare
to hand calculations for the exact volume (for the primitives). That's
basically what I did.
My intended use is with a lot of toruses in a chain, hopefully making it
possible to apply gravity to the chain and then drape it over arbitrary
fixed objects. That's why I used this formulation, as I don't anticipate
using fractal-based surfaces. If you can come up with something that
would handle *any* object and still be somewhat efficient, that would be
great!
Thanks for the input!
Chris
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
CShake wrote:
> True. I picked the 'orientation' of what direction to shoot the rays as
> somewhat arbitrary. For most objects it should work fine, but I am aware
> specifically that an image-based heightfield (since it normally rises up
> from the x-z plane) will create a large number of intersections, and
> would be more efficient to shoot along y. Considering the likelihood of
> this sort of comparison, I may actually switch it to shooting along y,
> since for 'normal' object it wouldn't make much of a difference.
It might be a good idea to do that. It may be worthwhile putting into the usage
notes that the following construction works fine to rotate the object just for
use with this macro. (I checked it.)
FindObjectCenter(object{BoxObj rotate<1,2,3>})
> I did have in mind a few different ways of shooting rays, specifically
> doing a sort of spherical shoot where theta and phi are the control
> variables, but evenly distributing the rays around the surface (and
> finding the volume of each 'chunk') would be a little more complicated.
> At least this method works in k*n^2, where k is something like the
> number of voids in the object, which is slightly more optimized than n^3
> in the case of most simpler objects (CSG, etc). I thought about trying
> to use an adaptive sampling method, but the extra code would be
> significant (and I don't have a good handle on how to do it yet, so it
> would involve more research too).
> It could be helpful to do a sort of basic test if the object is a
> primitive, and then use the appropriate mathematical formula, but the
> amount of time you'd run into that 'in the wild' would be nearly 0.
Initially I thought it might be worthwhile to keep track of how many times the
#while(more_sections... loop is iterated for each ray, then give up and use
another method if the number of iterations exceeds 10 or 100 times the number
of divisions. But I agree, it would probably be better to just post a
different method and let the scene designer choose which one to use. Yours is
likely to be superior for 99 out of a hundred uses, but it didn't look too
difficult to provide an alternative for that poor slob with the oddball case.
The main reason I considered writing an alternative is that I like to use noise
functions to make rock-shaped isosurface objects. Although the ones I usually
need are simple bumpy rocks (so I'll end up using your macro), I seem to
remember getting the idea after seeing rather complex asteroid-like or volcanic
ejecta-shaped objects and thought someone may need to be able to determine their
masses for some kind of animation...
I got this alternative version working, and as expected it is much slower and
less accurate than yours for simple shapes (with a given number of x and y
divisions, and a similar number of z divisions for my version). I'll test it
with more gnarly objects, see if I can find some cases where the alternative
gives better performance, and post the results and the finished code here.
Thanks for providing your test objects. Since I make so many mistakes, I always
prefer to verify known working cases. :)
Also, thank you for letting me start with your nice, clean code. It really
helped to get it working quickly. (Proper attribution is given in the code
comments, of course.)
~David
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
> True. I picked the 'orientation' of what direction to shoot the rays as
> somewhat arbitrary. For most objects it should work fine, but I am aware
> specifically that an image-based heightfield (since it normally rises up
> from the x-z plane) will create a large number of intersections, and
> would be more efficient to shoot along y. Considering the likelihood of
> this sort of comparison, I may actually switch it to shooting along y,
> since for 'normal' object it wouldn't make much of a difference.
>
The orientation of the tests can be passed as a parameter,somewhat like
the direction vector used by trace, or a global variable.
In the later case, you can have a default orientation that is used if
the user don't set it.
That way, if the object been tested have a special shape, the direction
can be chosen to be more effecient. After all, not all hightfields are
used horizontaly.
Alain
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
CShake schrieb:
> I did have in mind a few different ways of shooting rays, specifically
> doing a sort of spherical shoot where theta and phi are the control
> variables, but evenly distributing the rays around the surface (and
> finding the volume of each 'chunk') would be a little more complicated.
You could start from a platonic polyhedron as an approximation - a cube
would be sufficiently "spherical" for this purpose I guess :-)
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
waggy wrote:
> I got this alternative version working, and as expected it is much slower and
> less accurate than yours for simple shapes (with a given number of x and y
> divisions, and a similar number of z divisions for my version). I'll test it
> with more gnarly objects, see if I can find some cases where the alternative
> gives better performance, and post the results and the finished code here.
>
> Thanks for providing your test objects. Since I make so many mistakes, I always
> prefer to verify known working cases. :)
I've now written my moment of inertia macro, and couldn't find an
optimized way to do it in n^2, so here it is in n^3. My issue here is
that I'm getting the wrong values for the moments for everything but a
box, and I can't figure out why. Spheres and Toruses are giving almost
exactly half the official value. This macro would have even less of an
audience than my last, but again any math help would be appreciated.
I need to read up more on how to use moment of inertial tensors for
calculating the moment about an arbitrary axis, and I left my dynamics
textbooks at home this semester...
Also, this macro requires the center of mass to already be computed.
Chris
(I'm interchanging volume and mass here in the terminology, sorry. It
will be better once density is involved.)
// Macro, description, and test values:
#declare X_Resolution=20;
#declare Y_Resolution=20;
#declare Z_Resolution=20;
// Finds the moment of inertia for an arbitrary object.
// Return value is <Ix,Iy,Iz> where each is the moment around that axis
where it passes
// through the center of mass of the object, allowing for "easy"
calculation of any moment
// with the parallel axis theorem
#macro FindMomentsOfInertia(Object_Identifier, Center_Of_Mass)
// Set up reference vars
#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 dx = abs(max_x-min_x)/X_Resolution;
#local dy = abs(max_y-min_y)/Y_Resolution;
#local dz = abs(max_z-min_z)/Z_Resolution;
#local dv = dx*dy*dz;
#local Moments = <0,0,0>;
#local Volume = 0;
// Go through entire 3D matrix
#local curr_x = min_x+dx/2;
#while (curr_x < max_x)
#local curr_y = min_y+dy/2;
#while (curr_y < max_y)
#local curr_z = min_z+dz/2;
#while (curr_z < max_z)
#if(inside(Object_Identifier,<curr_x,curr_y,curr_z>))
// I = sum(m*r^2) where m is each point mass, r is distance
to axis
#local Volume = Volume+dv;
#local Moments = Moments+dv*<pow(abs(Center_Of_Mass.x-curr_x),2),
pow(abs(Center_Of_Mass.y-curr_y),2),
pow(abs(Center_Of_Mass.z-curr_z),2)>;
#end
#local curr_z = curr_z+dz;
#end
#local curr_y = curr_y+dy;
#end
#local curr_x = curr_x+dx;
#end
//(Moments)
DebugVector3("Calculated moments",Moments)
#debug concat("Calculated volume: ",str(Volume,0,2),"\n")
#end
#declare BoxObj=box{-1,1}
#declare SphereObj=sphere{0,1.5}
#declare TorusObj=torus{1,0.5}
#debug "Box: \n"
#debug concat("Real moments: ",str(1/12*4*(pow(2,2)+pow(2,2)),0,2),"\n")
#debug concat("Real volume: ",str(pow(2,3),0,2),"\n")
FindMomentsOfInertia(BoxObj,<0,0,0>)
#debug "Sphere: \n"
#local sphereMass = 4/3*pi*pow(1.5,3);
#debug concat("Real moments: ",str(2/5*sphereMass*pow(1.5,2),0,2),"\n")
#debug concat("Real volume: ",str(sphereMass,0,2),"\n")
FindMomentsOfInertia(SphereObj,<0,0,0>)
#debug "Torus: \n"
#local torusMass = 2*pow(pi,2)*1*pow(0.5,2);
#debug concat("Real moments:
",str(torusMass*(pow(1,2)+3/4*pow(0.5,2)),0,2),0,2),"\n")
#debug concat("Real volume: ",str(torusMass,0,2),"\n")
FindMomentsOfInertia(TorusObj,<0,0,0>)
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
| |
|
|
Alain wrote:
> The orientation of the tests can be passed as a parameter,somewhat like
> the direction vector used by trace, or a global variable.
> In the later case, you can have a default orientation that is used if
> the user don't set it.
>
> That way, if the object been tested have a special shape, the direction
> can be chosen to be more effecient. After all, not all hightfields are
> used horizontaly.
>
> Alain
Good idea, I'll look at implementing that. It will be a lot more work to
make it an arbitrary axis (with how I'm using the grid from one side of
the bounding box), but I figure that a choice of x,y,z should be enough
for most uses. Anything more and it's just as easy to rotate the object,
call the macro, then rotate it back and transform the returned location.
Chrsi
Post a reply to this message
|
|
| |
| |
|
|
|
|
| |
|
|
