Fix T41019: Calculate Mass does not calculate actual volume.
authorLukas Tönne <lukas.toenne@gmail.com>
Fri, 11 Jul 2014 10:06:13 +0000 (12:06 +0200)
committerLukas Tönne <lukas.toenne@gmail.com>
Fri, 11 Jul 2014 10:16:32 +0000 (12:16 +0200)
This was a ToDo item, for mesh-based rigid body shapes (trimesh, convex)
the operator was simply using the bounding box volume, which can grossly
overestimate the volume and mass.

Calculating the actual volume of a mesh is not so difficult after all,
see e.g.
http://research.microsoft.com/en-us/um/people/chazhang/publications/icip01_ChaZhang.pdf

This patch also allows calculating the center-of-mass in the same way.
This is currently unused, because the rigid body system assumes the CoM
to be the same as the geometric object center. This is fine most of the
time, adding such user settings for "center-of-mass offset" would also
add quite a bit of complexity in user space, but it could be necessary
at some point. A number of other physical properties could be calculated
using the same principle, e.g. the moment of inertia.

source/blender/blenkernel/BKE_mesh.h
source/blender/blenkernel/BKE_rigidbody.h
source/blender/blenkernel/intern/mesh_evaluate.c
source/blender/blenkernel/intern/rigidbody.c
source/blender/blenlib/BLI_math_geom.h
source/blender/blenlib/intern/math_geom.c
source/blender/editors/physics/rigidbody_object.c

index 587dea5095ca853357d0422af8ce981b768e161c..975a8fd63df744f4312754441d47a974478c7f94 100644 (file)
@@ -212,6 +212,10 @@ bool BKE_mesh_center_median(struct Mesh *me, float cent[3]);
 bool BKE_mesh_center_bounds(struct Mesh *me, float cent[3]);
 bool BKE_mesh_center_centroid(struct Mesh *me, float cent[3]);
 
+void BKE_mesh_calc_volume(struct MVert *mverts, int numVerts,
+                          struct MFace *mfaces, int numFaces,
+                          float *r_vol, float *r_com);
+
 /* tessface */
 void BKE_mesh_loops_to_mface_corners(
         struct CustomData *fdata, struct CustomData *ldata,
index 86be3bfe77015955e9174d799aeba1b4d3450a9c..c946f3ac9e881ffae6ee35ea05ec81d531381a38 100644 (file)
@@ -69,6 +69,9 @@ void BKE_rigidbody_world_groups_relink(struct RigidBodyWorld *rbw);
 /* 'validate' (i.e. make new or replace old) Physics-Engine objects */
 void BKE_rigidbody_validate_sim_world(struct Scene *scene, struct RigidBodyWorld *rbw, bool rebuild);
 
+void BKE_rigidbody_calc_volume(struct Object *ob, float *r_vol);
+void BKE_rigidbody_calc_center_of_mass(struct Object *ob, float r_com[3]);
+
 /* -------------- */
 /* Utilities */
 
index cb0386b120378ffcb9a6c049fe64c9a267dc4f7a..7e547ec4c2553e9cc363d0be3ac901c91fdcdc53 100644 (file)
@@ -1094,6 +1094,125 @@ bool BKE_mesh_center_centroid(Mesh *me, float cent[3])
 /** \} */
 
 
+/* -------------------------------------------------------------------- */
+
+/** \name Mesh Volume Calculation
+ * \{ */
+
+static bool mesh_calc_center_centroid_ex(MVert *mverts, int UNUSED(numVerts),
+                                         MFace *mfaces, int numFaces,
+                                         float center[3])
+{
+       float totweight;
+       int f;
+       
+       zero_v3(center);
+       
+       if (numFaces == 0)
+               return false;
+       
+       totweight = 0.0f;
+       for (f = 0; f < numFaces; ++f) {
+               MFace *face = &mfaces[f];
+               MVert *v1 = &mverts[face->v1];
+               MVert *v2 = &mverts[face->v2];
+               MVert *v3 = &mverts[face->v3];
+               MVert *v4 = &mverts[face->v4];
+               float area;
+               
+               area = area_tri_v3(v1->co, v2->co, v3->co);
+               madd_v3_v3fl(center, v1->co, area);
+               madd_v3_v3fl(center, v2->co, area);
+               madd_v3_v3fl(center, v3->co, area);
+               totweight += area;
+               
+               if (face->v4) {
+                       area = area_tri_v3(v3->co, v4->co, v1->co);
+                       madd_v3_v3fl(center, v3->co, area);
+                       madd_v3_v3fl(center, v4->co, area);
+                       madd_v3_v3fl(center, v1->co, area);
+                       totweight += area;
+               }
+       }
+       if (totweight == 0.0f)
+               return false;
+       
+       mul_v3_fl(center, 1.0f / (3.0f * totweight));
+       
+       return true;
+}
+
+void BKE_mesh_calc_volume(MVert *mverts, int numVerts,
+                          MFace *mfaces, int numFaces,
+                          float *r_vol, float *r_com)
+{
+       float center[3];
+       float totvol;
+       int f;
+       
+       if (r_vol) *r_vol = 0.0f;
+       if (r_com) zero_v3(r_com);
+       
+       if (numFaces == 0)
+               return;
+       
+       if (!mesh_calc_center_centroid_ex(mverts, numVerts, mfaces, numFaces, center))
+               return;
+       
+       totvol = 0.0f;
+       for (f = 0; f < numFaces; ++f) {
+               MFace *face = &mfaces[f];
+               MVert *v1 = &mverts[face->v1];
+               MVert *v2 = &mverts[face->v2];
+               MVert *v3 = &mverts[face->v3];
+               MVert *v4 = &mverts[face->v4];
+               float vol;
+               
+               vol = volume_tetrahedron_signed_v3(center, v1->co, v2->co, v3->co);
+               if (r_vol) {
+                       totvol += vol;
+               }
+               if (r_com) {
+                       /* averaging factor 1/4 is applied in the end */
+                       madd_v3_v3fl(r_com, center, vol); // XXX could extract this
+                       madd_v3_v3fl(r_com, v1->co, vol);
+                       madd_v3_v3fl(r_com, v2->co, vol);
+                       madd_v3_v3fl(r_com, v3->co, vol);
+               }
+               
+               if (face->v4) {
+                       vol = volume_tetrahedron_signed_v3(center, v3->co, v4->co, v1->co);
+                       
+                       if (r_vol) {
+                               totvol += vol;
+                       }
+                       if (r_com) {
+                               /* averaging factor 1/4 is applied in the end */
+                               madd_v3_v3fl(r_com, center, vol); // XXX could extract this
+                               madd_v3_v3fl(r_com, v3->co, vol);
+                               madd_v3_v3fl(r_com, v4->co, vol);
+                               madd_v3_v3fl(r_com, v1->co, vol);
+                       }
+               }
+       }
+       
+       /* Note: Depending on arbitrary centroid position,
+        * totvol can become negative even for a valid mesh.
+        * The true value is always the positive value.
+        */
+       if (r_vol) {
+               *r_vol = fabsf(totvol);
+       }
+       if (r_com) {
+               /* Note: Factor 1/4 is applied once for all vertices here.
+                * This also automatically negates the vector if totvol is negative.
+                */
+               if (totvol != 0.0f)
+                       mul_v3_fl(r_com, 0.25f / totvol);
+       }
+}
+
+
 /* -------------------------------------------------------------------- */
 
 /** \name NGon Tessellation (NGon/Tessface Conversion)
index 09293503ad5a4d636a818cd1a7f7e1154cfcf10e..2aaf8adea57d8ecf1b9cdbfe0c04a24c94f6d244 100644 (file)
@@ -57,6 +57,7 @@
 #include "BKE_effect.h"
 #include "BKE_global.h"
 #include "BKE_library.h"
+#include "BKE_mesh.h"
 #include "BKE_object.h"
 #include "BKE_pointcache.h"
 #include "BKE_rigidbody.h"
@@ -455,6 +456,182 @@ static void rigidbody_validate_sim_shape(Object *ob, bool rebuild)
 
 /* --------------------- */
 
+/* helper function to calculate volume of rigidbody object */
+// TODO: allow a parameter to specify method used to calculate this?
+void BKE_rigidbody_calc_volume(Object *ob, float *r_vol)
+{
+       RigidBodyOb *rbo = ob->rigidbody_object;
+
+       float size[3]  = {1.0f, 1.0f, 1.0f};
+       float radius = 1.0f;
+       float height = 1.0f;
+
+       float volume = 0.0f;
+
+       /* if automatically determining dimensions, use the Object's boundbox
+        *      - assume that all quadrics are standing upright on local z-axis
+        *      - assume even distribution of mass around the Object's pivot
+        *        (i.e. Object pivot is centralised in boundbox)
+        *      - boundbox gives full width
+        */
+       // XXX: all dimensions are auto-determined now... later can add stored settings for this
+       BKE_object_dimensions_get(ob, size);
+
+       if (ELEM3(rbo->shape, RB_SHAPE_CAPSULE, RB_SHAPE_CYLINDER, RB_SHAPE_CONE)) {
+               /* take radius as largest x/y dimension, and height as z-dimension */
+               radius = MAX2(size[0], size[1]) * 0.5f;
+               height = size[2];
+       }
+       else if (rbo->shape == RB_SHAPE_SPHERE) {
+               /* take radius to the the largest dimension to try and encompass everything */
+               radius = max_fff(size[0], size[1], size[2]) * 0.5f;
+       }
+
+       /* calculate volume as appropriate  */
+       switch (rbo->shape) {
+               case RB_SHAPE_BOX:
+                       volume = size[0] * size[1] * size[2];
+                       break;
+
+               case RB_SHAPE_SPHERE:
+                       volume = 4.0f / 3.0f * (float)M_PI * radius * radius * radius;
+                       break;
+
+               /* for now, assume that capsule is close enough to a cylinder... */
+               case RB_SHAPE_CAPSULE:
+               case RB_SHAPE_CYLINDER:
+                       volume = (float)M_PI * radius * radius * height;
+                       break;
+
+               case RB_SHAPE_CONE:
+                       volume = (float)M_PI / 3.0f * radius * radius * height;
+                       break;
+
+               case RB_SHAPE_CONVEXH:
+               case RB_SHAPE_TRIMESH:
+               {
+                       if (ob->type == OB_MESH) {
+                               DerivedMesh *dm = rigidbody_get_mesh(ob);
+                               MVert *mvert;
+                               MFace *mface;
+                               int totvert, totface;
+                               
+                               /* ensure mesh validity, then grab data */
+                               if (dm == NULL)
+                                       return;
+                       
+                               DM_ensure_tessface(dm);
+                       
+                               mvert   = (dm) ? dm->getVertArray(dm) : NULL;
+                               totvert = (dm) ? dm->getNumVerts(dm) : 0;
+                               mface   = (dm) ? dm->getTessFaceArray(dm) : NULL;
+                               totface = (dm) ? dm->getNumTessFaces(dm) : 0;
+                               
+                               if (totvert > 0 && totface > 0) {
+                                       BKE_mesh_calc_volume(mvert, totvert, mface, totface, &volume, NULL);
+                               }
+                               
+                               /* cleanup temp data */
+                               if (dm && ob->rigidbody_object->mesh_source == RBO_MESH_BASE) {
+                                       dm->release(dm);
+                               }
+                       }
+                       else {
+                               /* rough estimate from boundbox as fallback */
+                               /* XXX could implement other types of geometry here (curves, etc.) */
+                               volume = size[0] * size[1] * size[2];
+                       }
+                       break;
+               }
+
+#if 0 // XXX: not defined yet
+               case RB_SHAPE_COMPOUND:
+                       volume = 0.0f;
+                       break;
+#endif
+       }
+
+       /* return the volume calculated */
+       if (r_vol) *r_vol = volume;
+}
+
+void BKE_rigidbody_calc_center_of_mass(Object *ob, float r_com[3])
+{
+       RigidBodyOb *rbo = ob->rigidbody_object;
+
+       float size[3]  = {1.0f, 1.0f, 1.0f};
+       float height = 1.0f;
+
+       zero_v3(r_com);
+
+       /* if automatically determining dimensions, use the Object's boundbox
+        *      - assume that all quadrics are standing upright on local z-axis
+        *      - assume even distribution of mass around the Object's pivot
+        *        (i.e. Object pivot is centralised in boundbox)
+        *      - boundbox gives full width
+        */
+       // XXX: all dimensions are auto-determined now... later can add stored settings for this
+       BKE_object_dimensions_get(ob, size);
+
+       /* calculate volume as appropriate  */
+       switch (rbo->shape) {
+               case RB_SHAPE_BOX:
+               case RB_SHAPE_SPHERE:
+               case RB_SHAPE_CAPSULE:
+               case RB_SHAPE_CYLINDER:
+                       break;
+
+               case RB_SHAPE_CONE:
+                       /* take radius as largest x/y dimension, and height as z-dimension */
+                       height = size[2];
+                       /* cone is geometrically centered on the median,
+                        * center of mass is 1/4 up from the base
+                        */
+                       r_com[2] = -0.25f * height;
+                       break;
+
+               case RB_SHAPE_CONVEXH:
+               case RB_SHAPE_TRIMESH:
+               {
+                       if (ob->type == OB_MESH) {
+                               DerivedMesh *dm = rigidbody_get_mesh(ob);
+                               MVert *mvert;
+                               MFace *mface;
+                               int totvert, totface;
+                               
+                               /* ensure mesh validity, then grab data */
+                               if (dm == NULL)
+                                       return;
+                       
+                               DM_ensure_tessface(dm);
+                       
+                               mvert   = (dm) ? dm->getVertArray(dm) : NULL;
+                               totvert = (dm) ? dm->getNumVerts(dm) : 0;
+                               mface   = (dm) ? dm->getTessFaceArray(dm) : NULL;
+                               totface = (dm) ? dm->getNumTessFaces(dm) : 0;
+                               
+                               if (totvert > 0 && totface > 0) {
+                                       BKE_mesh_calc_volume(mvert, totvert, mface, totface, NULL, r_com);
+                               }
+                               
+                               /* cleanup temp data */
+                               if (dm && ob->rigidbody_object->mesh_source == RBO_MESH_BASE) {
+                                       dm->release(dm);
+                               }
+                       }
+                       break;
+               }
+
+#if 0 // XXX: not defined yet
+               case RB_SHAPE_COMPOUND:
+                       volume = 0.0f;
+                       break;
+#endif
+       }
+}
+
+/* --------------------- */
+
 /**
  * Create physics sim representation of object given RigidBody settings
  *
@@ -1414,6 +1591,8 @@ struct RigidBodyOb *BKE_rigidbody_copy_object(Object *ob) { return NULL; }
 struct RigidBodyCon *BKE_rigidbody_copy_constraint(Object *ob) { return NULL; }
 void BKE_rigidbody_relink_constraint(RigidBodyCon *rbc) {}
 void BKE_rigidbody_validate_sim_world(Scene *scene, RigidBodyWorld *rbw, bool rebuild) {}
+void BKE_rigidbody_calc_volume(Object *ob, float *r_vol) { if (r_vol) *r_vol = 0.0f; }
+void BKE_rigidbody_calc_center_of_mass(Object *ob, float r_com[3]) { zero_v3(r_com); }
 struct RigidBodyWorld *BKE_rigidbody_create_world(Scene *scene) { return NULL; }
 struct RigidBodyWorld *BKE_rigidbody_world_copy(RigidBodyWorld *rbw) { return NULL; }
 void BKE_rigidbody_world_groups_relink(struct RigidBodyWorld *rbw) {}
index f4bcc81084615abf8cb682fd57910ba78f52c16c..ed0777aceea654cf009ad131ad34e4dcdb938e65 100644 (file)
@@ -72,6 +72,7 @@ MINLINE float plane_point_side_v3(const float plane[4], const float co[3]);
 /********************************* Volume **********************************/
 
 float volume_tetrahedron_v3(const float v1[3], const float v2[3], const float v3[3], const float v4[3]);
+float volume_tetrahedron_signed_v3(const float v1[3], const float v2[3], const float v3[3], const float v4[3]);
 
 bool is_quad_convex_v3(const float v1[3], const float v2[3], const float v3[3], const float v4[3]);
 bool is_quad_convex_v2(const float v1[2], const float v2[2], const float v3[2], const float v4[2]);
index 5f3ab5eb73e823cc70b29a2057aaa6ce89475e12..efb881cde1b108ef7624eb41daecf06d9d3ebd53 100644 (file)
@@ -238,6 +238,18 @@ float volume_tetrahedron_v3(const float v1[3], const float v2[3], const float v3
        return fabsf(determinant_m3_array(m)) / 6.0f;
 }
 
+/**
+ * The volume from a tetrahedron, normal pointing inside gives negative volume
+ */
+float volume_tetrahedron_signed_v3(const float v1[3], const float v2[3], const float v3[3], const float v4[3])
+{
+       float m[3][3];
+       sub_v3_v3v3(m[0], v1, v2);
+       sub_v3_v3v3(m[1], v2, v3);
+       sub_v3_v3v3(m[2], v3, v4);
+       return determinant_m3_array(m) / 6.0f;
+}
+
 
 /********************************* Distance **********************************/
 
index b46d79fc82ea2cfedb312544ba1f5342bf0f2083..13a3d7a523f6ee7f117d389cfb25e32a12cd3317 100644 (file)
@@ -485,78 +485,6 @@ static EnumPropertyItem *rigidbody_materials_itemf(bContext *UNUSED(C), PointerR
 
 /* ------------------------------------------ */
 
-/* helper function to calculate volume of rigidbody object */
-// TODO: allow a parameter to specify method used to calculate this?
-static float rigidbody_object_calc_volume(Object *ob)
-{
-       RigidBodyOb *rbo = ob->rigidbody_object;
-
-       float size[3]  = {1.0f, 1.0f, 1.0f};
-       float radius = 1.0f;
-       float height = 1.0f;
-
-       float volume = 0.0f;
-
-       /* if automatically determining dimensions, use the Object's boundbox
-        *      - assume that all quadrics are standing upright on local z-axis
-        *      - assume even distribution of mass around the Object's pivot
-        *        (i.e. Object pivot is centralised in boundbox)
-        *      - boundbox gives full width
-        */
-       // XXX: all dimensions are auto-determined now... later can add stored settings for this
-       BKE_object_dimensions_get(ob, size);
-
-       if (ELEM3(rbo->shape, RB_SHAPE_CAPSULE, RB_SHAPE_CYLINDER, RB_SHAPE_CONE)) {
-               /* take radius as largest x/y dimension, and height as z-dimension */
-               radius = MAX2(size[0], size[1]) * 0.5f;
-               height = size[2];
-       }
-       else if (rbo->shape == RB_SHAPE_SPHERE) {
-               /* take radius to the the largest dimension to try and encompass everything */
-               radius = max_fff(size[0], size[1], size[2]) * 0.5f;
-       }
-
-       /* calculate volume as appropriate  */
-       switch (rbo->shape) {
-               case RB_SHAPE_BOX:
-                       volume = size[0] * size[1] * size[2];
-                       break;
-
-               case RB_SHAPE_SPHERE:
-                       volume = 4.0f / 3.0f * (float)M_PI * radius * radius * radius;
-                       break;
-
-               /* for now, assume that capsule is close enough to a cylinder... */
-               case RB_SHAPE_CAPSULE:
-               case RB_SHAPE_CYLINDER:
-                       volume = (float)M_PI * radius * radius * height;
-                       break;
-
-               case RB_SHAPE_CONE:
-                       volume = (float)M_PI / 3.0f * radius * radius * height;
-                       break;
-
-               /* for now, all mesh shapes are just treated as boxes...
-                * NOTE: this may overestimate the volume, but other methods are overkill
-                */
-               case RB_SHAPE_CONVEXH:
-               case RB_SHAPE_TRIMESH:
-                       volume = size[0] * size[1] * size[2];
-                       break;
-
-#if 0 // XXX: not defined yet
-               case RB_SHAPE_COMPOUND:
-                       volume = 0.0f;
-                       break;
-#endif
-       }
-
-       /* return the volume calculated */
-       return volume;
-}
-
-/* ------------------------------------------ */
-
 static int rigidbody_objects_calc_mass_exec(bContext *C, wmOperator *op)
 {
        int material = RNA_enum_get(op->ptr, "material");
@@ -589,7 +517,7 @@ static int rigidbody_objects_calc_mass_exec(bContext *C, wmOperator *op)
                        /* mass is calculated from the approximate volume of the object,
                         * and the density of the material we're simulating
                         */
-                       volume = rigidbody_object_calc_volume(ob);
+                       BKE_rigidbody_calc_volume(ob, &volume);
                        mass = volume * density;
 
                        /* use RNA-system to change the property and perform all necessary changes */