From 98a99cb763949b6a73c8dec76e9a7f3a604f8862 Mon Sep 17 00:00:00 2001 From: Campbell Barton Date: Wed, 22 Jan 2020 01:24:19 +1100 Subject: Fix T73295: Incorrect BMesh volume calculation Use double precision since volume calculation is susceptible to float precision errors. --- source/blender/bmesh/intern/bmesh_query.c | 34 ++++++++++++++++++++++--------- source/blender/bmesh/intern/bmesh_query.h | 2 +- 2 files changed, 25 insertions(+), 11 deletions(-) (limited to 'source/blender/bmesh') diff --git a/source/blender/bmesh/intern/bmesh_query.c b/source/blender/bmesh/intern/bmesh_query.c index a121bd3b9f4..a9b0d4d3918 100644 --- a/source/blender/bmesh/intern/bmesh_query.c +++ b/source/blender/bmesh/intern/bmesh_query.c @@ -2477,39 +2477,53 @@ bool BM_face_is_normal_valid(const BMFace *f) return len_squared_v3v3(no, f->no) < (eps * eps); } -static void bm_mesh_calc_volume_face(const BMFace *f, float *r_vol) +/** + * Use to accumulate volume calculation for faces with consistent winding. + * + * Use double precision since this is prone to float precision error, see T73295. + */ +static double bm_mesh_calc_volume_face(const BMFace *f) { const int tottri = f->len - 2; BMLoop **loops = BLI_array_alloca(loops, f->len); uint(*index)[3] = BLI_array_alloca(index, tottri); - int j; + double vol = 0.0; BM_face_calc_tessellation(f, false, loops, index); - for (j = 0; j < tottri; j++) { + for (int j = 0; j < tottri; j++) { const float *p1 = loops[index[j][0]]->v->co; const float *p2 = loops[index[j][1]]->v->co; const float *p3 = loops[index[j][2]]->v->co; + double p1_db[3]; + double p2_db[3]; + double p3_db[3]; + + copy_v3db_v3fl(p1_db, p1); + copy_v3db_v3fl(p2_db, p2); + copy_v3db_v3fl(p3_db, p3); + /* co1.dot(co2.cross(co3)) / 6.0 */ - float cross[3]; - cross_v3_v3v3(cross, p2, p3); - *r_vol += (1.0f / 6.0f) * dot_v3v3(p1, cross); + double cross[3]; + cross_v3_v3v3_db(cross, p2_db, p3_db); + vol += dot_v3v3_db(p1_db, cross); } + return (1.0 / 6.0) * vol; } -float BM_mesh_calc_volume(BMesh *bm, bool is_signed) +double BM_mesh_calc_volume(BMesh *bm, bool is_signed) { /* warning, calls own tessellation function, may be slow */ - float vol = 0.0f; + double vol = 0.0; BMFace *f; BMIter fiter; BM_ITER_MESH (f, &fiter, bm, BM_FACES_OF_MESH) { - bm_mesh_calc_volume_face(f, &vol); + vol += bm_mesh_calc_volume_face(f); } if (is_signed == false) { - vol = fabsf(vol); + vol = fabs(vol); } return vol; diff --git a/source/blender/bmesh/intern/bmesh_query.h b/source/blender/bmesh/intern/bmesh_query.h index 134e0b99691..fb646b11076 100644 --- a/source/blender/bmesh/intern/bmesh_query.h +++ b/source/blender/bmesh/intern/bmesh_query.h @@ -233,7 +233,7 @@ bool BM_face_is_any_edge_flag_test(const BMFace *f, const char hflag) ATTR_WARN_ bool BM_face_is_normal_valid(const BMFace *f) ATTR_WARN_UNUSED_RESULT ATTR_NONNULL(); -float BM_mesh_calc_volume(BMesh *bm, bool is_signed) ATTR_WARN_UNUSED_RESULT ATTR_NONNULL(); +double BM_mesh_calc_volume(BMesh *bm, bool is_signed) ATTR_WARN_UNUSED_RESULT ATTR_NONNULL(); int BM_mesh_calc_face_groups(BMesh *bm, int *r_groups_array, -- cgit v1.2.3