diff options
Diffstat (limited to 'source/blender/blenkernel/intern/mesh_remap.c')
| -rw-r--r-- | source/blender/blenkernel/intern/mesh_remap.c | 295 |
1 files changed, 244 insertions, 51 deletions
diff --git a/source/blender/blenkernel/intern/mesh_remap.c b/source/blender/blenkernel/intern/mesh_remap.c index 08823c1c4a1..c818c6b9f19 100644 --- a/source/blender/blenkernel/intern/mesh_remap.c +++ b/source/blender/blenkernel/intern/mesh_remap.c @@ -51,6 +51,250 @@ /* -------------------------------------------------------------------- */ +/** \name Some generic helpers. + * \{ */ + +static bool mesh_remap_bvhtree_query_nearest( + BVHTreeFromMesh *treedata, BVHTreeNearest *nearest, + const float co[3], const float max_dist_sq, float *r_hit_dist) +{ + /* Use local proximity heuristics (to reduce the nearest search). */ + if (nearest->index != -1) { + nearest->dist_sq = min_ff(len_squared_v3v3(co, nearest->co), max_dist_sq); + } + else { + nearest->dist_sq = max_dist_sq; + } + /* Compute and store result. If invalid (-1 index), keep FLT_MAX dist. */ + BLI_bvhtree_find_nearest(treedata->tree, co, nearest, treedata->nearest_callback, treedata); + + if ((nearest->index != -1) && (nearest->dist_sq <= max_dist_sq)) { + *r_hit_dist = sqrtf(nearest->dist_sq); + return true; + } + else { + return false; + } +} + +static bool mesh_remap_bvhtree_query_raycast( + BVHTreeFromMesh *treedata, BVHTreeRayHit *rayhit, + const float co[3], const float no[3], const float radius, const float max_dist, float *r_hit_dist) +{ + BVHTreeRayHit rayhit_tmp; + float inv_no[3]; + + rayhit->index = -1; + rayhit->dist = max_dist; + BLI_bvhtree_ray_cast(treedata->tree, co, no, radius, rayhit, treedata->raycast_callback, treedata); + + /* Also cast in the other direction! */ + rayhit_tmp = *rayhit; + negate_v3_v3(inv_no, no); + BLI_bvhtree_ray_cast(treedata->tree, co, inv_no, radius, &rayhit_tmp, treedata->raycast_callback, treedata); + if (rayhit_tmp.dist < rayhit->dist) { + *rayhit = rayhit_tmp; + } + + if ((rayhit->index != -1) && (rayhit->dist <= max_dist)) { + *r_hit_dist = rayhit->dist; + return true; + } + else { + return false; + } +} + +/** \} */ + +/** + * \name Auto-match. + * + * Find transform of a mesh to get best match with another. + * \{ */ + +/** + * Compute a value of the difference between both given meshes. + * The smaller the result, the better the match. + * + * We return the inverse of the average of the inversed shortest distance from each dst vertex to src ones. + * In other words, beyond a certain (relatively small) distance, all differences have more or less the same weight + * in final result, which allows to reduce influence of a few high differences, in favor of a global good matching. + */ +float BKE_mesh_remap_calc_difference_from_dm( + const SpaceTransform *space_transform, const MVert *verts_dst, const int numverts_dst, DerivedMesh *dm_src) +{ + BVHTreeFromMesh treedata = {NULL}; + BVHTreeNearest nearest = {0}; + float hit_dist; + + float result = 0.0f; + int i; + + bvhtree_from_mesh_verts(&treedata, dm_src, 0.0f, 2, 6); + nearest.index = -1; + + for (i = 0; i < numverts_dst; i++) { + float tmp_co[3]; + + copy_v3_v3(tmp_co, verts_dst[i].co); + + /* Convert the vertex to tree coordinates, if needed. */ + if (space_transform) { + BLI_space_transform_apply(space_transform, tmp_co); + } + + if (mesh_remap_bvhtree_query_nearest(&treedata, &nearest, tmp_co, FLT_MAX, &hit_dist)) { + result += 1.0f / (hit_dist + 1.0f); + } + else { + /* No source for this dest vertex! */ + result += 1e-18f; + } + } + + result = ((float)numverts_dst / result) - 1.0f; + +// printf("%s: Computed difference between meshes (the lower the better): %f\n", __func__, result); + + return result; +} + +/* This helper computes the eigen values & vectors for covariance matrix of all given vertices coordinates. + * + * Those vectors define the 'average ellipsoid' of the mesh (i.e. the 'best fitting' ellipsoid + * containing 50% of the vertices). + * + * Note that it will not perform fantastic in case two or more eigen values are equal (e.g. a cylinder or + * parallelepiped with a square section give two identical eigenvalues, a sphere or tetrahedron give + * three identical ones, etc.), since you cannot really define all axes in those cases. We default to dummy + * generated orthogonal vectors in this case, instead of using eigen vectors. + */ +static void mesh_calc_eigen_matrix( + const MVert *verts, const float (*vcos)[3], const int numverts, float r_mat[4][4]) +{ + float center[3], covmat[3][3]; + float eigen_val[3], eigen_vec[3][3]; + float (*cos)[3] = (float (*)[3])vcos; + + bool eigen_success; + int i; + + if (verts) { + const MVert *mv; + float (*co)[3]; + + cos = MEM_mallocN(sizeof(*cos) * (size_t)numverts, __func__); + for (i = 0, co = cos, mv = verts; i < numverts; i++, co++, mv++) { + copy_v3_v3(*co, mv->co); + } + } + unit_m4(r_mat); + + /* Note: here we apply sample correction to covariance matrix, since we consider the vertices as a sample + * of the whole 'surface' population of our mesh... */ + BLI_covariance_m3_v3n(cos, numverts, true, covmat, center); + + eigen_success = BLI_eigen_solve_selfadjoint_m3(covmat, eigen_val, eigen_vec); + BLI_assert(eigen_success); + UNUSED_VARS_NDEBUG(eigen_success); + + /* Special handling of cases where some eigen values are (nearly) identical. */ + if (compare_ff_relative(eigen_val[0], eigen_val[1], FLT_EPSILON, 64)) { + if (compare_ff_relative(eigen_val[0], eigen_val[2], FLT_EPSILON, 64)) { + /* No preferred direction, that set of vertices has a spherical average, + * so we simply returned scaled/translated identity matrix (with no rotation). */ + unit_m3(eigen_vec); + } + else { + /* Ellipsoid defined by eigen values/vectors has a spherical section, + * we can only define one axis from eigen_vec[2] (two others computed eigen vecs + * are not so nice for us here, they tend to 'randomly' rotate around valid one). + * Note that eigen vectors as returned by BLI_eigen_solve_selfadjoint_m3() are normalized. */ + ortho_basis_v3v3_v3(eigen_vec[0], eigen_vec[1], eigen_vec[2]); + } + } + else if (compare_ff_relative(eigen_val[0], eigen_val[2], FLT_EPSILON, 64)) { + /* Same as above, but with eigen_vec[1] as valid axis. */ + ortho_basis_v3v3_v3(eigen_vec[2], eigen_vec[0], eigen_vec[1]); + } + else if (compare_ff_relative(eigen_val[1], eigen_val[2], FLT_EPSILON, 64)) { + /* Same as above, but with eigen_vec[0] as valid axis. */ + ortho_basis_v3v3_v3(eigen_vec[1], eigen_vec[2], eigen_vec[0]); + } + + for (i = 0; i < 3; i++) { + float evi = eigen_val[i]; + + /* Protect against 1D/2D degenerated cases! */ + /* Note: not sure why we need square root of eigen values here (which are equivalent to singular values, + * as far as I have understood), but it seems to heavily reduce (if not completly nullify) + * the error due to non-uniform scalings... */ + evi = (evi < 1e-6f && evi > -1e-6f) ? ((evi < 0.0f) ? -1e-3f : 1e-3f) : sqrtf_signed(evi); + mul_v3_fl(eigen_vec[i], evi); + } + + copy_m4_m3(r_mat, eigen_vec); + copy_v3_v3(r_mat[3], center); + + if (verts) { + MEM_freeN(cos); + } +} + +/** + * Set r_space_transform so that best bbox of dst matches best bbox of src. + */ +void BKE_mesh_remap_find_best_match_from_dm( + const MVert *verts_dst, const int numverts_dst, DerivedMesh *dm_src, SpaceTransform *r_space_transform) +{ + /* Note that those are done so that we successively get actual mirror matrix (by multiplication of columns)... */ + const float mirrors[][3] = { + {-1.0f, 1.0f, 1.0f}, /* -> -1, 1, 1 */ + { 1.0f, -1.0f, 1.0f}, /* -> -1, -1, 1 */ + { 1.0f, 1.0f, -1.0f}, /* -> -1, -1, -1 */ + { 1.0f, -1.0f, 1.0f}, /* -> -1, 1, -1 */ + {-1.0f, 1.0f, 1.0f}, /* -> 1, 1, -1 */ + { 1.0f, -1.0f, 1.0f}, /* -> 1, -1, -1 */ + { 1.0f, 1.0f, -1.0f}, /* -> 1, -1, 1 */ + {0.0f, 0.0f, 0.0f}, + }; + const float (*mirr)[3]; + + float mat_src[4][4], mat_dst[4][4], best_mat_dst[4][4]; + float best_match = FLT_MAX, match; + + const int numverts_src = dm_src->getNumVerts(dm_src); + float (*vcos_src)[3] = MEM_mallocN(sizeof(*vcos_src) * (size_t)numverts_src, __func__); + dm_src->getVertCos(dm_src, vcos_src); + + mesh_calc_eigen_matrix(NULL, vcos_src, numverts_src, mat_src); + mesh_calc_eigen_matrix(verts_dst, NULL, numverts_dst, mat_dst); + + BLI_space_transform_global_from_matrices(r_space_transform, mat_dst, mat_src); + match = BKE_mesh_remap_calc_difference_from_dm(r_space_transform, verts_dst, numverts_dst, dm_src); + best_match = match; + copy_m4_m4(best_mat_dst, mat_dst); + + /* And now, we have to check the otehr sixth possible mirrored versions... */ + for (mirr = mirrors; (*mirr)[0]; mirr++) { + mul_v3_fl(mat_dst[0], (*mirr)[0]); + mul_v3_fl(mat_dst[1], (*mirr)[1]); + mul_v3_fl(mat_dst[2], (*mirr)[2]); + + BLI_space_transform_global_from_matrices(r_space_transform, mat_dst, mat_src); + match = BKE_mesh_remap_calc_difference_from_dm(r_space_transform, verts_dst, numverts_dst, dm_src); + if (match < best_match) { + best_match = match; + copy_m4_m4(best_mat_dst, mat_dst); + } + } + + BLI_space_transform_global_from_matrices(r_space_transform, best_mat_dst, mat_src); +} + +/** \} */ + /** \name Mesh to mesh mapping * \{ */ @@ -147,57 +391,6 @@ static int mesh_remap_interp_poly_data_get( return sources_num; } -static bool mesh_remap_bvhtree_query_nearest( - BVHTreeFromMesh *treedata, BVHTreeNearest *nearest, - float co[3], const float max_dist_sq, float *r_hit_dist) -{ - /* Use local proximity heuristics (to reduce the nearest search). */ - if (nearest->index != -1) { - nearest->dist_sq = min_ff(len_squared_v3v3(co, nearest->co), max_dist_sq); - } - else { - nearest->dist_sq = max_dist_sq; - } - /* Compute and store result. If invalid (-1 index), keep FLT_MAX dist. */ - BLI_bvhtree_find_nearest(treedata->tree, co, nearest, treedata->nearest_callback, treedata); - - if ((nearest->index != -1) && (nearest->dist_sq <= max_dist_sq)) { - *r_hit_dist = sqrtf(nearest->dist_sq); - return true; - } - else { - return false; - } -} - -static bool mesh_remap_bvhtree_query_raycast( - BVHTreeFromMesh *treedata, BVHTreeRayHit *rayhit, - const float co[3], const float no[3], const float radius, const float max_dist, float *r_hit_dist) -{ - BVHTreeRayHit rayhit_tmp; - float inv_no[3]; - - rayhit->index = -1; - rayhit->dist = max_dist; - BLI_bvhtree_ray_cast(treedata->tree, co, no, radius, rayhit, treedata->raycast_callback, treedata); - - /* Also cast in the other direction! */ - rayhit_tmp = *rayhit; - negate_v3_v3(inv_no, no); - BLI_bvhtree_ray_cast(treedata->tree, co, inv_no, radius, &rayhit_tmp, treedata->raycast_callback, treedata); - if (rayhit_tmp.dist < rayhit->dist) { - *rayhit = rayhit_tmp; - } - - if ((rayhit->index != -1) && (rayhit->dist <= max_dist)) { - *r_hit_dist = rayhit->dist; - return true; - } - else { - return false; - } -} - /* Little helper when dealing with source islands */ typedef struct IslandResult { float factor; /* A factor, based on which best island for a given set of elements will be selected. */ |