From 5ea27bec1fc1e65ce85cccd3079883a41970cd6c Mon Sep 17 00:00:00 2001 From: Campbell Barton Date: Tue, 14 Jun 2016 16:47:05 +1000 Subject: BMesh Decimate: use doubles to calculate optimized position This allows the error threshold for calculating the optimized location to be much lower. Resolves visible artifacts w/ 1m-tri happy-buddha example. --- source/blender/blenlib/BLI_quadric.h | 7 +- source/blender/blenlib/intern/quadric.c | 86 ++++++++++++++++------ .../blender/bmesh/tools/bmesh_decimate_collapse.c | 44 +++++++---- 3 files changed, 94 insertions(+), 43 deletions(-) (limited to 'source') diff --git a/source/blender/blenlib/BLI_quadric.h b/source/blender/blenlib/BLI_quadric.h index eaf9c7a0738..d46a221d108 100644 --- a/source/blender/blenlib/BLI_quadric.h +++ b/source/blender/blenlib/BLI_quadric.h @@ -39,8 +39,7 @@ typedef struct Quadric { /* conversion */ void BLI_quadric_from_plane(Quadric *q, const double v[4]); -void BLI_quadric_to_tensor_m3(const Quadric *q, float m[3][3]); -void BLI_quadric_to_vector_v3(const Quadric *q, float v[3]); +void BLI_quadric_to_vector_v3(const Quadric *q, double v[3]); void BLI_quadric_clear(Quadric *q); @@ -50,7 +49,7 @@ void BLI_quadric_add_qu_ququ(Quadric *r, const Quadric *a, const Quadric *b); void BLI_quadric_mul(Quadric *a, const double scalar); /* solve */ -double BLI_quadric_evaluate(const Quadric *q, const float v_fl[3]); -bool BLI_quadric_optimize(const Quadric *q, float v[3], const float epsilon); +double BLI_quadric_evaluate(const Quadric *q, const double v[3]); +bool BLI_quadric_optimize(const Quadric *q, double v[3], const double epsilon); #endif /* __BLI_QUADRIC_H__ */ diff --git a/source/blender/blenlib/intern/quadric.c b/source/blender/blenlib/intern/quadric.c index 588cd9c2cb5..c9d27494455 100644 --- a/source/blender/blenlib/intern/quadric.c +++ b/source/blender/blenlib/intern/quadric.c @@ -63,26 +63,68 @@ void BLI_quadric_from_plane(Quadric *q, const double v[4]) q->d2 = v[3] * v[3]; } -void BLI_quadric_to_tensor_m3(const Quadric *q, float m[3][3]) +#if 0 /* UNUSED */ + +static void quadric_to_tensor_m3(const Quadric *q, double m[3][3]) { - m[0][0] = (float)q->a2; - m[0][1] = (float)q->ab; - m[0][2] = (float)q->ac; + m[0][0] = q->a2; + m[0][1] = q->ab; + m[0][2] = q->ac; - m[1][0] = (float)q->ab; - m[1][1] = (float)q->b2; - m[1][2] = (float)q->bc; + m[1][0] = q->ab; + m[1][1] = q->b2; + m[1][2] = q->bc; - m[2][0] = (float)q->ac; - m[2][1] = (float)q->bc; - m[2][2] = (float)q->c2; + m[2][0] = q->ac; + m[2][1] = q->bc; + m[2][2] = q->c2; } -void BLI_quadric_to_vector_v3(const Quadric *q, float v[3]) +#endif + +/** + * Inline inverse matrix creation. + * Equivalent of: + * + * \code{.c} + * quadric_to_tensor_m3(q, m); + * invert_m3_db(m, eps); + * \endcode + */ +static bool quadric_to_tensor_m3_inverse(const Quadric *q, double m[3][3], double epsilon) { - v[0] = (float)q->ad; - v[1] = (float)q->bd; - v[2] = (float)q->cd; + const double det = + (q->a2 * (q->b2 * q->c2 - q->bc * q->bc) - + q->ab * (q->ab * q->c2 - q->ac * q->bc) + + q->ac * (q->ab * q->bc - q->ac * q->b2)); + + if (fabs(det) > epsilon) { + const double invdet = 1.0 / det; + + m[0][0] = (q->b2 * q->c2 - q->bc * q->bc) * invdet; + m[1][0] = (q->bc * q->ac - q->ab * q->c2) * invdet; + m[2][0] = (q->ab * q->bc - q->b2 * q->ac) * invdet; + + m[0][1] = (q->ac * q->bc - q->ab * q->c2) * invdet; + m[1][1] = (q->a2 * q->c2 - q->ac * q->ac) * invdet; + m[2][1] = (q->ab * q->ac - q->a2 * q->bc) * invdet; + + m[0][2] = (q->ab * q->bc - q->ac * q->b2) * invdet; + m[1][2] = (q->ac * q->ab - q->a2 * q->bc) * invdet; + m[2][2] = (q->a2 * q->b2 - q->ab * q->ab) * invdet; + + return true; + } + else { + return false; + } +} + +void BLI_quadric_to_vector_v3(const Quadric *q, double v[3]) +{ + v[0] = q->ad; + v[1] = q->bd; + v[2] = q->cd; } void BLI_quadric_clear(Quadric *q) @@ -105,26 +147,22 @@ void BLI_quadric_mul(Quadric *a, const double scalar) mul_vn_db((double *)a, QUADRIC_FLT_TOT, scalar); } -double BLI_quadric_evaluate(const Quadric *q, const float v_fl[3]) +double BLI_quadric_evaluate(const Quadric *q, const double v[3]) { - const double v[3] = {UNPACK3(v_fl)}; return ((q->a2 * v[0] * v[0]) + (q->ab * 2 * v[0] * v[1]) + (q->ac * 2 * v[0] * v[2]) + (q->ad * 2 * v[0]) + (q->b2 * v[1] * v[1]) + (q->bc * 2 * v[1] * v[2]) + (q->bd * 2 * v[1]) + (q->c2 * v[2] * v[2]) + (q->cd * 2 * v[2]) + (q->d2)); } -bool BLI_quadric_optimize(const Quadric *q, float v[3], const float epsilon) +bool BLI_quadric_optimize(const Quadric *q, double v[3], const double epsilon) { - float m[3][3]; + double m[3][3]; - BLI_quadric_to_tensor_m3(q, m); - - if (invert_m3_ex(m, epsilon)) { + if (quadric_to_tensor_m3_inverse(q, m, epsilon)) { BLI_quadric_to_vector_v3(q, v); - mul_m3_v3(m, v); - negate_v3(v); - + mul_m3_v3_db(m, v); + negate_v3_db(v); return true; } else { diff --git a/source/blender/bmesh/tools/bmesh_decimate_collapse.c b/source/blender/bmesh/tools/bmesh_decimate_collapse.c index fe8b132a2a5..3ea7dae536b 100644 --- a/source/blender/bmesh/tools/bmesh_decimate_collapse.c +++ b/source/blender/bmesh/tools/bmesh_decimate_collapse.c @@ -68,7 +68,9 @@ #endif #define BOUNDARY_PRESERVE_WEIGHT 100.0f -#define OPTIMIZE_EPS 0.01f /* FLT_EPSILON is too small, see [#33106] */ +/* Uses double precision, impacts behavior on near-flat surfaces, + * cane give issues with very small faces. 1e-2 is too big, see: T48154. */ +#define OPTIMIZE_EPS 1e-8 #define COST_INVALID FLT_MAX typedef enum CD_UseFlag { @@ -140,8 +142,8 @@ static void bm_decim_build_quadrics(BMesh *bm, Quadric *vquadrics) } -static void bm_decim_calc_target_co( - BMEdge *e, float optimize_co[3], +static void bm_decim_calc_target_co_db( + BMEdge *e, double optimize_co[3], const Quadric *vquadrics) { /* compute an edge contraction target for edge 'e' @@ -158,10 +160,22 @@ static void bm_decim_calc_target_co( return; /* all is good */ } else { - mid_v3_v3v3(optimize_co, e->v1->co, e->v2->co); + optimize_co[0] = 0.5 * ((double)e->v1->co[0] + (double)e->v2->co[0]); + optimize_co[1] = 0.5 * ((double)e->v1->co[1] + (double)e->v2->co[1]); + optimize_co[2] = 0.5 * ((double)e->v1->co[2] + (double)e->v2->co[2]); } } +static void bm_decim_calc_target_co_fl( + BMEdge *e, float optimize_co[3], + const Quadric *vquadrics) +{ + double optimize_co_db[3]; + bm_decim_calc_target_co_db(e, optimize_co_db, vquadrics); + copy_v3fl_v3db(optimize_co, optimize_co_db); +} + + static bool bm_edge_collapse_is_degenerate_flip(BMEdge *e, const float optimize_co[3]) { BMIter liter; @@ -240,8 +254,6 @@ static void bm_decim_build_edge_cost_single( const float *vweights, const float vweight_factor, Heap *eheap, HeapNode **eheap_table) { - const Quadric *q1, *q2; - float optimize_co[3]; float cost; if (eheap_table[BM_elem_index_get(e)]) { @@ -279,15 +291,17 @@ static void bm_decim_build_edge_cost_single( } /* end sanity check */ + { + double optimize_co[3]; + bm_decim_calc_target_co_db(e, optimize_co, vquadrics); - bm_decim_calc_target_co(e, optimize_co, vquadrics); - - q1 = &vquadrics[BM_elem_index_get(e->v1)]; - q2 = &vquadrics[BM_elem_index_get(e->v2)]; - - cost = (BLI_quadric_evaluate(q1, optimize_co) + - BLI_quadric_evaluate(q2, optimize_co)); + const Quadric *q1, *q2; + q1 = &vquadrics[BM_elem_index_get(e->v1)]; + q2 = &vquadrics[BM_elem_index_get(e->v2)]; + cost = (BLI_quadric_evaluate(q1, optimize_co) + + BLI_quadric_evaluate(q2, optimize_co)); + } /* note, 'cost' shouldn't be negative but happens sometimes with small values. * this can cause faces that make up a flat surface to over-collapse, see [#37121] */ @@ -1155,7 +1169,7 @@ static bool bm_decim_edge_collapse( return false; } - bm_decim_calc_target_co(e, optimize_co, vquadrics); + bm_decim_calc_target_co_fl(e, optimize_co, vquadrics); /* check if this would result in an overlapping face */ if (UNLIKELY(bm_edge_collapse_is_degenerate_flip(e, optimize_co))) { @@ -1426,7 +1440,7 @@ void BM_mesh_decimate_collapse( goto invalidate; } - bm_decim_calc_target_co(e, optimize_co, vquadrics); + bm_decim_calc_target_co_fl(e, optimize_co, vquadrics); if (e_index_mirr == e_index) { optimize_co[symmetry_axis] = 0.0f; -- cgit v1.2.3