diff options
Diffstat (limited to 'source/blender')
-rw-r--r-- | source/blender/blenkernel/BKE_shrinkwrap.h | 34 | ||||
-rw-r--r-- | source/blender/blenkernel/intern/DerivedMesh.c | 12 | ||||
-rw-r--r-- | source/blender/blenkernel/intern/constraint.c | 6 | ||||
-rw-r--r-- | source/blender/blenkernel/intern/mesh.c | 1 | ||||
-rw-r--r-- | source/blender/blenkernel/intern/mesh_runtime.c | 2 | ||||
-rw-r--r-- | source/blender/blenkernel/intern/shrinkwrap.c | 542 | ||||
-rw-r--r-- | source/blender/blenlib/BLI_math_solvers.h | 9 | ||||
-rw-r--r-- | source/blender/blenlib/intern/math_solvers.c | 95 | ||||
-rw-r--r-- | source/blender/depsgraph/DEG_depsgraph.h | 6 | ||||
-rw-r--r-- | source/blender/depsgraph/intern/builder/deg_builder_relations.cc | 3 | ||||
-rw-r--r-- | source/blender/makesdna/DNA_mesh_types.h | 3 | ||||
-rw-r--r-- | source/blender/makesdna/DNA_modifier_types.h | 1 | ||||
-rw-r--r-- | source/blender/makesrna/intern/rna_constraint.c | 3 | ||||
-rw-r--r-- | source/blender/makesrna/intern/rna_modifier.c | 3 | ||||
-rw-r--r-- | source/blender/modifiers/intern/MOD_shrinkwrap.c | 6 |
15 files changed, 716 insertions, 10 deletions
diff --git a/source/blender/blenkernel/BKE_shrinkwrap.h b/source/blender/blenkernel/BKE_shrinkwrap.h index d1a14003b04..04c0f1c5d9a 100644 --- a/source/blender/blenkernel/BKE_shrinkwrap.h +++ b/source/blender/blenkernel/BKE_shrinkwrap.h @@ -33,6 +33,7 @@ /* Shrinkwrap stuff */ #include "BKE_bvhutils.h" +#include "BLI_bitmap.h" /* * Shrinkwrap is composed by a set of functions and options that define the type of shrink. @@ -55,6 +56,34 @@ struct ShrinkwrapModifierData; struct BVHTree; struct SpaceTransform; +/* Information about boundary edges in the mesh. */ +typedef struct ShrinkwrapBoundaryVertData { + /* Average direction of edges that meet here. */ + float direction[3]; + + /* Closest vector to direction that is orthogonal to vertex normal. */ + float normal_plane[3]; +} ShrinkwrapBoundaryVertData; + +typedef struct ShrinkwrapBoundaryData { + /* True if the edge belongs to exactly one face. */ + const BLI_bitmap *edge_is_boundary; + /* True if the looptri has any boundary edges. */ + const BLI_bitmap *looptri_has_boundary; + + /* Mapping from vertex index to boundary vertex index, or -1. + * Used for compact storage of data about boundary vertices. */ + const int *vert_boundary_id; + unsigned int num_boundary_verts; + + /* Direction data about boundary vertices. */ + const ShrinkwrapBoundaryVertData *boundary_verts; +} ShrinkwrapBoundaryData; + +void BKE_shrinkwrap_discard_boundary_data(struct Mesh *mesh); +void BKE_shrinkwrap_compute_boundary_data(struct Mesh *mesh); + +/* Information about a mesh and BVH tree. */ typedef struct ShrinkwrapTreeData { Mesh *mesh; @@ -62,6 +91,7 @@ typedef struct ShrinkwrapTreeData { BVHTreeFromMesh treeData; float (*clnors)[3]; + ShrinkwrapBoundaryData *boundary; } ShrinkwrapTreeData; /* Checks if the modifier needs target normals with these settings. */ @@ -91,6 +121,10 @@ bool BKE_shrinkwrap_project_normal( char options, const float vert[3], const float dir[3], const float ray_radius, const struct SpaceTransform *transf, struct ShrinkwrapTreeData *tree, BVHTreeRayHit *hit); +/* Maps the point to the nearest surface, either by simple nearest, or by target normal projection. */ +void BKE_shrinkwrap_find_nearest_surface( + struct ShrinkwrapTreeData *tree, struct BVHTreeNearest *nearest, float co[3], int type); + /* Computes a smooth normal of the target (if applicable) at the hit location. */ void BKE_shrinkwrap_compute_smooth_normal( const struct ShrinkwrapTreeData *tree, const struct SpaceTransform *transform, diff --git a/source/blender/blenkernel/intern/DerivedMesh.c b/source/blender/blenkernel/intern/DerivedMesh.c index c84bd378b04..1d787388774 100644 --- a/source/blender/blenkernel/intern/DerivedMesh.c +++ b/source/blender/blenkernel/intern/DerivedMesh.c @@ -75,6 +75,7 @@ #include "DEG_depsgraph.h" #include "DEG_depsgraph_query.h" +#include "BKE_shrinkwrap.h" #ifdef WITH_OPENSUBDIV # include "DNA_userdef_types.h" @@ -1980,6 +1981,15 @@ static void mesh_finalize_eval(Object *object) } } +static void mesh_build_extra_data(struct Depsgraph *depsgraph, Object *ob) +{ + uint32_t eval_flags = DEG_get_eval_flags_for_id(depsgraph, &ob->id); + + if (eval_flags & DAG_EVAL_NEED_SHRINKWRAP_BOUNDARY) { + BKE_shrinkwrap_compute_boundary_data(ob->runtime.mesh_eval); + } +} + static void mesh_build_data( struct Depsgraph *depsgraph, Scene *scene, Object *ob, CustomDataMask dataMask, const bool build_shapekey_layers, const bool need_mapping) @@ -2016,6 +2026,8 @@ static void mesh_build_data( } BLI_assert(!(ob->derivedFinal->dirty & DM_DIRTY_NORMALS)); + + mesh_build_extra_data(depsgraph, ob); } static void editbmesh_build_data( diff --git a/source/blender/blenkernel/intern/constraint.c b/source/blender/blenkernel/intern/constraint.c index 41b07d73dc9..cc94e87c747 100644 --- a/source/blender/blenkernel/intern/constraint.c +++ b/source/blender/blenkernel/intern/constraint.c @@ -3636,6 +3636,7 @@ static void shrinkwrap_get_tarmat(struct Depsgraph *UNUSED(depsgraph), bConstrai switch (scon->shrinkType) { case MOD_SHRINKWRAP_NEAREST_SURFACE: case MOD_SHRINKWRAP_NEAREST_VERTEX: + case MOD_SHRINKWRAP_TARGET_PROJECT: { BVHTreeNearest nearest; @@ -3644,15 +3645,14 @@ static void shrinkwrap_get_tarmat(struct Depsgraph *UNUSED(depsgraph), bConstrai BLI_space_transform_apply(&transform, co); - BVHTreeFromMesh *treeData = &tree.treeData; - BLI_bvhtree_find_nearest(treeData->tree, co, &nearest, treeData->nearest_callback, treeData); + BKE_shrinkwrap_find_nearest_surface(&tree, &nearest, co, scon->shrinkType); if (nearest.index < 0) { fail = true; break; } - if (scon->shrinkType == MOD_SHRINKWRAP_NEAREST_SURFACE) { + if (scon->shrinkType != MOD_SHRINKWRAP_NEAREST_VERTEX) { if (do_track_normal) { track_normal = true; BKE_shrinkwrap_compute_smooth_normal(&tree, NULL, nearest.index, nearest.co, nearest.no, track_no); diff --git a/source/blender/blenkernel/intern/mesh.c b/source/blender/blenkernel/intern/mesh.c index ed01889d200..6a00aaf576b 100644 --- a/source/blender/blenkernel/intern/mesh.c +++ b/source/blender/blenkernel/intern/mesh.c @@ -563,6 +563,7 @@ void BKE_mesh_copy_data(Main *bmain, Mesh *me_dst, const Mesh *me_src, const int me_dst->runtime.batch_cache = NULL; me_dst->runtime.looptris.array = NULL; me_dst->runtime.bvh_cache = NULL; + me_dst->runtime.shrinkwrap_data = NULL; if (me_src->id.tag & LIB_TAG_NO_MAIN) { me_dst->runtime.deformed_only = me_src->runtime.deformed_only; diff --git a/source/blender/blenkernel/intern/mesh_runtime.c b/source/blender/blenkernel/intern/mesh_runtime.c index 68880cacfb9..0931318f17a 100644 --- a/source/blender/blenkernel/intern/mesh_runtime.c +++ b/source/blender/blenkernel/intern/mesh_runtime.c @@ -44,6 +44,7 @@ #include "BKE_mesh.h" #include "BKE_mesh_runtime.h" #include "BKE_subdiv_ccg.h" +#include "BKE_shrinkwrap.h" /* -------------------------------------------------------------------- */ /** \name Mesh Runtime Struct Utils @@ -202,6 +203,7 @@ void BKE_mesh_runtime_clear_geometry(Mesh *mesh) BKE_subdiv_ccg_destroy(mesh->runtime.subdiv_ccg); mesh->runtime.subdiv_ccg = NULL; } + BKE_shrinkwrap_discard_boundary_data(mesh); } /** \} */ diff --git a/source/blender/blenkernel/intern/shrinkwrap.c b/source/blender/blenkernel/intern/shrinkwrap.c index 6942c8f1da3..78045e5fd56 100644 --- a/source/blender/blenkernel/intern/shrinkwrap.c +++ b/source/blender/blenkernel/intern/shrinkwrap.c @@ -45,6 +45,7 @@ #include "BLI_math.h" #include "BLI_utildefines.h" #include "BLI_task.h" +#include "BLI_math_solvers.h" #include "BKE_shrinkwrap.h" #include "BKE_cdderivedmesh.h" @@ -57,6 +58,9 @@ #include "BKE_editmesh.h" #include "BKE_mesh.h" /* for OMP limits. */ #include "BKE_subsurf.h" +#include "BKE_mesh_runtime.h" + +#include "MEM_guardedalloc.h" #include "BLI_strict_flags.h" @@ -70,7 +74,6 @@ /* Util macros */ #define OUT_OF_MEMORY() ((void)printf("Shrinkwrap: Out of memory\n")) - typedef struct ShrinkwrapCalcData { ShrinkwrapModifierData *smd; //shrinkwrap modifier data @@ -106,7 +109,8 @@ typedef struct ShrinkwrapCalcCBData { /* Checks if the modifier needs target normals with these settings. */ bool BKE_shrinkwrap_needs_normals(int shrinkType, int shrinkMode) { - return shrinkType != MOD_SHRINKWRAP_NEAREST_VERTEX && shrinkMode == MOD_SHRINKWRAP_ABOVE_SURFACE; + return (shrinkType == MOD_SHRINKWRAP_TARGET_PROJECT) || + (shrinkType != MOD_SHRINKWRAP_NEAREST_VERTEX && shrinkMode == MOD_SHRINKWRAP_ABOVE_SURFACE); } /* Initializes the mesh data structure from the given mesh and settings. */ @@ -142,6 +146,10 @@ bool BKE_shrinkwrap_init_tree(ShrinkwrapTreeData *data, Mesh *mesh, int shrinkTy } } + if (shrinkType == MOD_SHRINKWRAP_TARGET_PROJECT) { + data->boundary = mesh->runtime.shrinkwrap_data; + } + return true; } } @@ -152,6 +160,176 @@ void BKE_shrinkwrap_free_tree(ShrinkwrapTreeData *data) free_bvhtree_from_mesh(&data->treeData); } +/* Free boundary data for target project */ +void BKE_shrinkwrap_discard_boundary_data(struct Mesh *mesh) +{ + struct ShrinkwrapBoundaryData *data = mesh->runtime.shrinkwrap_data; + + if (data != NULL) { + MEM_freeN((void*)data->edge_is_boundary); + MEM_freeN((void*)data->looptri_has_boundary); + MEM_freeN((void*)data->vert_boundary_id); + MEM_freeN((void*)data->boundary_verts); + + MEM_freeN(data); + } + + mesh->runtime.shrinkwrap_data = NULL; +} + +/* Accumulate edge for average boundary edge direction. */ +static void merge_vert_dir(ShrinkwrapBoundaryVertData *vdata, signed char *status, int index, const float edge_dir[3], signed char side) +{ + BLI_assert(index >= 0); + float *direction = vdata[index].direction; + + /* Invert the direction vector if either: + * - This is the second edge and both edges have the vertex as start or end. + * - For third and above, if it points in the wrong direction. + */ + if (status[index] >= 0 ? status[index] == side : dot_v3v3(direction, edge_dir) < 0) { + sub_v3_v3(direction, edge_dir); + } + else { + add_v3_v3(direction, edge_dir); + } + + status[index] = (status[index] == 0) ? side : -1; +} + +static ShrinkwrapBoundaryData *shrinkwrap_build_boundary_data(struct Mesh *mesh) +{ + const MLoop *mloop = mesh->mloop; + const MEdge *medge = mesh->medge; + const MVert *mvert = mesh->mvert; + + /* Count faces per edge (up to 2). */ + char *edge_mode = MEM_calloc_arrayN((size_t)mesh->totedge, sizeof(char), __func__); + + for (int i = 0; i < mesh->totloop; i++) { + unsigned int eidx = mloop[i].e; + + if (edge_mode[eidx] < 2) { + edge_mode[eidx]++; + } + } + + /* Build the boundary edge bitmask. */ + BLI_bitmap *edge_is_boundary = BLI_BITMAP_NEW(mesh->totedge, "ShrinkwrapBoundaryData::edge_is_boundary"); + unsigned int num_boundary_edges = 0; + + for (int i = 0; i < mesh->totedge; i++) { + edge_mode[i] = (edge_mode[i] == 1); + + if (edge_mode[i]) { + BLI_BITMAP_ENABLE(edge_is_boundary, i); + num_boundary_edges++; + } + } + + /* If no boundary, return NULL. */ + if (num_boundary_edges == 0) { + MEM_freeN(edge_is_boundary); + MEM_freeN(edge_mode); + return NULL; + } + + /* Allocate the data object. */ + ShrinkwrapBoundaryData *data = MEM_callocN(sizeof(ShrinkwrapBoundaryData), "ShrinkwrapBoundaryData"); + + data->edge_is_boundary = edge_is_boundary; + + /* Build the boundary looptri bitmask. */ + const MLoopTri *mlooptri = BKE_mesh_runtime_looptri_ensure(mesh); + int totlooptri = BKE_mesh_runtime_looptri_len(mesh); + + BLI_bitmap *looptri_has_boundary = BLI_BITMAP_NEW(totlooptri, "ShrinkwrapBoundaryData::looptri_is_boundary"); + + for (int i = 0; i < totlooptri; i++) { + int edges[3]; + BKE_mesh_looptri_get_real_edges(mesh, &mlooptri[i], edges); + + for (int j = 0; j < 3; j++) { + if (edges[j] >= 0 && edge_mode[edges[j]]) { + BLI_BITMAP_ENABLE(looptri_has_boundary, i); + break; + } + } + } + + data->looptri_has_boundary = looptri_has_boundary; + + /* Find boundary vertices and build a mapping table for compact storage of data. */ + int *vert_boundary_id = MEM_calloc_arrayN((size_t)mesh->totvert, sizeof(int), "ShrinkwrapBoundaryData::vert_boundary_id"); + + for (int i = 0; i < mesh->totedge; i++) { + if (edge_mode[i]) { + const MEdge *edge = &medge[i]; + + vert_boundary_id[edge->v1] = 1; + vert_boundary_id[edge->v2] = 1; + } + } + + unsigned int num_boundary_verts = 0; + + for (int i = 0; i < mesh->totvert; i++) { + vert_boundary_id[i] = (vert_boundary_id[i] != 0) ? (int)num_boundary_verts++ : -1; + } + + data->vert_boundary_id = vert_boundary_id; + data->num_boundary_verts = num_boundary_verts; + + /* Compute average directions. */ + ShrinkwrapBoundaryVertData *boundary_verts = MEM_calloc_arrayN(num_boundary_verts, sizeof(*boundary_verts), "ShrinkwrapBoundaryData::boundary_verts"); + + signed char *vert_status = MEM_calloc_arrayN(num_boundary_verts, sizeof(char), __func__); + + for (int i = 0; i < mesh->totedge; i++) { + if (edge_mode[i]) { + const MEdge *edge = &medge[i]; + + float dir[3]; + sub_v3_v3v3(dir, mvert[edge->v2].co, mvert[edge->v1].co); + normalize_v3(dir); + + merge_vert_dir(boundary_verts, vert_status, vert_boundary_id[edge->v1], dir, 1); + merge_vert_dir(boundary_verts, vert_status, vert_boundary_id[edge->v2], dir, 2); + } + } + + MEM_freeN(vert_status); + + /* Finalize average direction and compute normal. */ + for (int i = 0; i < mesh->totvert; i++) { + int bidx = vert_boundary_id[i]; + + if (bidx >= 0) { + ShrinkwrapBoundaryVertData *vdata = &boundary_verts[bidx]; + float no[3], tmp[3]; + + normalize_v3(vdata->direction); + + normal_short_to_float_v3(no, mesh->mvert[i].no); + cross_v3_v3v3(tmp, no, vdata->direction); + cross_v3_v3v3(vdata->normal_plane, tmp, no); + normalize_v3(vdata->normal_plane); + } + } + + data->boundary_verts = boundary_verts; + + MEM_freeN(edge_mode); + return data; +} + +void BKE_shrinkwrap_compute_boundary_data(struct Mesh *mesh) +{ + BKE_shrinkwrap_discard_boundary_data(mesh); + + mesh->runtime.shrinkwrap_data = shrinkwrap_build_boundary_data(mesh); +} + /* * Shrinkwrap to the nearest vertex * @@ -523,6 +701,350 @@ static void shrinkwrap_calc_normal_projection(ShrinkwrapCalcData *calc) } /* + * Shrinkwrap Target Surface Project mode + * + * It uses Newton's method to find a surface location with its + * smooth normal pointing at the original point. + * + * The equation system on barycentric weights and normal multiplier: + * + * (w0*V0 + w1*V1 + w2*V2) + l * (w0*N0 + w1*N1 + w2*N2) - CO = 0 + * w0 + w1 + w2 = 1 + * + * The actual solution vector is [ w0, w1, l ], with w2 eliminated. + */ + +//#define TRACE_TARGET_PROJECT + +typedef struct TargetProjectTriData { + const float **vtri_co; + const float (*vtri_no)[3]; + const float *point_co; + + float n0_minus_n2[3], n1_minus_n2[3]; + float c0_minus_c2[3], c1_minus_c2[3]; + + /* Current interpolated position and normal. */ + float co_interp[3], no_interp[3]; +} TargetProjectTriData; + +/* Computes the deviation of the equation system from goal. */ +static void target_project_tri_deviation(void *userdata, const float x[3], float r_delta[3]) +{ + TargetProjectTriData *data = userdata; + + float w[3] = { x[0], x[1], 1.0f - x[0] - x[1] }; + interp_v3_v3v3v3(data->co_interp, data->vtri_co[0], data->vtri_co[1], data->vtri_co[2], w); + interp_v3_v3v3v3(data->no_interp, data->vtri_no[0], data->vtri_no[1], data->vtri_no[2], w); + + madd_v3_v3v3fl(r_delta, data->co_interp, data->no_interp, x[2]); + sub_v3_v3(r_delta, data->point_co); +} + +/* Computes the Jacobian matrix of the equation system. */ +static void target_project_tri_jacobian(void *userdata, const float x[3], float r_jacobian[3][3]) +{ + TargetProjectTriData *data = userdata; + + madd_v3_v3v3fl(r_jacobian[0], data->c0_minus_c2, data->n0_minus_n2, x[2]); + madd_v3_v3v3fl(r_jacobian[1], data->c1_minus_c2, data->n1_minus_n2, x[2]); + + copy_v3_v3(r_jacobian[2], data->vtri_no[2]); + madd_v3_v3fl(r_jacobian[2], data->n0_minus_n2, x[0]); + madd_v3_v3fl(r_jacobian[2], data->n1_minus_n2, x[1]); +} + +/* Clamp barycentric weights to the triangle. */ +static void target_project_tri_clamp(float x[3]) +{ + if (x[0] < 0.0f) { + x[0] = 0.0f; + } + if (x[1] < 0.0f) { + x[1] = 0.0f; + } + if (x[0] + x[1] > 1.0f) { + x[0] = x[0] / (x[0] + x[1]); + x[1] = 1.0f - x[0]; + } +} + +/* Correct the Newton's method step to keep the coordinates within the triangle. */ +static bool target_project_tri_correct(void *UNUSED(userdata), const float x[3], float step[3], float x_next[3]) +{ + /* Insignificant correction threshold */ + const float epsilon = 1e-6f; + const float dir_epsilon = 0.05f; + bool fixed = false, locked = false; + + /* Weight 0 and 1 boundary check. */ + for (int i = 0; i < 2; i++) { + if (step[i] > x[i]) { + if (step[i] > dir_epsilon * fabsf(step[1 - i])) { + /* Abort if the solution is clearly outside the domain. */ + if (x[i] < epsilon) { + return false; + } + + /* Scale a significant step down to arrive at the boundary. */ + mul_v3_fl(step, x[i] / step[i]); + fixed = true; + } + else { + /* Reset precision errors to stay at the boundary. */ + step[i] = x[i]; + fixed = locked = true; + } + } + } + + /* Weight 2 boundary check. */ + float sum = x[0] + x[1]; + float sstep = step[0] + step[1]; + + if (sum - sstep > 1.0f) { + if (sstep < -dir_epsilon * (fabsf(step[0]) + fabsf(step[1]))) { + /* Abort if the solution is clearly outside the domain. */ + if (sum > 1.0f - epsilon) { + return false; + } + + /* Scale a significant step down to arrive at the boundary. */ + mul_v3_fl(step, (1.0f - sum) / -sstep); + fixed = true; + } + else { + /* Reset precision errors to stay at the boundary. */ + if (locked) { + step[0] = step[1] = 0.0f; + } + else { + step[0] -= 0.5f * sstep; + step[1] = -step[0]; + fixed = true; + } + } + } + + /* Recompute and clamp the new coordinates after step correction. */ + if (fixed) { + sub_v3_v3v3(x_next, x, step); + target_project_tri_clamp(x_next); + } + + return true; +} + +static bool target_project_solve_point_tri( + const float *vtri_co[3], const float vtri_no[3][3], const float point_co[3], + const float hit_co[3], float hit_dist_sq, + float r_hit_co[3], float r_hit_no[3]) +{ + float x[3], tmp[3]; + float dist = sqrtf(hit_dist_sq); + float epsilon = dist * 1.0e-5f; + + CLAMP_MIN(epsilon, 1.0e-5f); + + /* Initial solution vector: barycentric weights plus distance along normal. */ + interp_weights_tri_v3(x, UNPACK3(vtri_co), hit_co); + + interp_v3_v3v3v3(r_hit_no, UNPACK3(vtri_no), x); + sub_v3_v3v3(tmp, point_co, hit_co); + + x[2] = (dot_v3v3(tmp, r_hit_no) < 0) ? -dist : dist; + + /* Solve the equations iteratively. */ + TargetProjectTriData tri_data = { .vtri_co = vtri_co, .vtri_no = vtri_no, .point_co = point_co }; + + sub_v3_v3v3(tri_data.n0_minus_n2, vtri_no[0], vtri_no[2]); + sub_v3_v3v3(tri_data.n1_minus_n2, vtri_no[1], vtri_no[2]); + sub_v3_v3v3(tri_data.c0_minus_c2, vtri_co[0], vtri_co[2]); + sub_v3_v3v3(tri_data.c1_minus_c2, vtri_co[1], vtri_co[2]); + + target_project_tri_clamp(x); + +#ifdef TRACE_TARGET_PROJECT + const bool trace = true; +#else + const bool trace = false; +#endif + + bool ok = BLI_newton3d_solve(target_project_tri_deviation, target_project_tri_jacobian, target_project_tri_correct, &tri_data, epsilon, 20, trace, x, x); + + if (ok) { + copy_v3_v3(r_hit_co, tri_data.co_interp); + copy_v3_v3(r_hit_no, tri_data.no_interp); + + return true; + } + + return false; +} + +static bool update_hit(BVHTreeNearest *nearest, int index, const float co[3], const float hit_co[3], const float hit_no[3]) +{ + float dist_sq = len_squared_v3v3(hit_co, co); + + if (dist_sq < nearest->dist_sq) { +#ifdef TRACE_TARGET_PROJECT + printf("===> %d (%.3f,%.3f,%.3f) %g < %g\n", index, UNPACK3(hit_co), dist_sq, nearest->dist_sq); +#endif + nearest->index = index; + nearest->dist_sq = dist_sq; + copy_v3_v3(nearest->co, hit_co); + normalize_v3_v3(nearest->no, hit_no); + return true; + } + + return false; +} + +/* Target projection on a non-manifold boundary edge - treats it like an infinitely thin cylinder. */ +static void target_project_edge(const ShrinkwrapTreeData *tree, int index, const float co[3], BVHTreeNearest *nearest, int eidx) +{ + const BVHTreeFromMesh *data = &tree->treeData; + const MEdge *edge = &tree->mesh->medge[eidx]; + const float *vedge_co[2] = { data->vert[edge->v1].co, data->vert[edge->v2].co }; + +#ifdef TRACE_TARGET_PROJECT + printf("EDGE %d (%.3f,%.3f,%.3f) (%.3f,%.3f,%.3f)\n", eidx, UNPACK3(vedge_co[0]), UNPACK3(vedge_co[1])); +#endif + + /* Retrieve boundary vertex IDs */ + const int *vert_boundary_id = tree->boundary->vert_boundary_id; + int bid1 = vert_boundary_id[edge->v1], bid2 = vert_boundary_id[edge->v2]; + + if (bid1 < 0 || bid2 < 0) { + return; + } + + /* Retrieve boundary vertex normals and align them to direction. */ + const ShrinkwrapBoundaryVertData *boundary_verts = tree->boundary->boundary_verts; + float vedge_dir[2][3], dir[3]; + + copy_v3_v3(vedge_dir[0], boundary_verts[bid1].normal_plane); + copy_v3_v3(vedge_dir[1], boundary_verts[bid2].normal_plane); + + sub_v3_v3v3(dir, vedge_co[1], vedge_co[0]); + + if (dot_v3v3(boundary_verts[bid1].direction, dir) < 0) { + negate_v3(vedge_dir[0]); + } + if (dot_v3v3(boundary_verts[bid2].direction, dir) < 0) { + negate_v3(vedge_dir[1]); + } + + /* Solve a quadratic equation: lerp(d0,d1,x) * (co - lerp(v0,v1,x)) = 0 */ + float d0v0 = dot_v3v3(vedge_dir[0], vedge_co[0]), d0v1 = dot_v3v3(vedge_dir[0], vedge_co[1]); + float d1v0 = dot_v3v3(vedge_dir[1], vedge_co[0]), d1v1 = dot_v3v3(vedge_dir[1], vedge_co[1]); + float d0co = dot_v3v3(vedge_dir[0], co); + + float a = d0v1 - d0v0 + d1v0 - d1v1; + float b = 2 * d0v0 - d0v1 - d0co - d1v0 + dot_v3v3(vedge_dir[1], co); + float c = d0co - d0v0; + float det = b * b - 4 * a * c; + + if (det >= 0) { + const float epsilon = 1e-6f; + float sdet = sqrtf(det); + float hit_co[3], hit_no[3]; + + for (int i = (det > 0 ? 2 : 0); i >= 0; i -= 2) { + float x = (- b + ((float)i - 1) * sdet) / (2 * a); + + if (x >= -epsilon && x <= 1.0f + epsilon) { + CLAMP(x, 0, 1); + + float vedge_no[2][3]; + normal_short_to_float_v3(vedge_no[0], data->vert[edge->v1].no); + normal_short_to_float_v3(vedge_no[1], data->vert[edge->v2].no); + + interp_v3_v3v3(hit_co, vedge_co[0], vedge_co[1], x); + interp_v3_v3v3(hit_no, vedge_no[0], vedge_no[1], x); + + update_hit(nearest, index, co, hit_co, hit_no); + } + } + } +} + +/* Target normal projection BVH callback - based on mesh_looptri_nearest_point. */ +static void mesh_looptri_target_project(void *userdata, int index, const float co[3], BVHTreeNearest *nearest) +{ + const ShrinkwrapTreeData *tree = (ShrinkwrapTreeData *)userdata; + const BVHTreeFromMesh *data = &tree->treeData; + const MLoopTri *lt = &data->looptri[index]; + const MLoop *loop[3] = { &data->loop[lt->tri[0]], &data->loop[lt->tri[1]], &data->loop[lt->tri[2]] }; + const MVert *vtri[3] = { &data->vert[loop[0]->v], &data->vert[loop[1]->v], &data->vert[loop[2]->v] }; + const float *vtri_co[3] = { vtri[0]->co, vtri[1]->co, vtri[2]->co }; + float raw_hit_co[3], hit_co[3], hit_no[3], dist_sq, vtri_no[3][3]; + + /* First find the closest point and bail out if it's worse than the current solution. */ + closest_on_tri_to_point_v3(raw_hit_co, co, UNPACK3(vtri_co)); + dist_sq = len_squared_v3v3(co, raw_hit_co); + +#ifdef TRACE_TARGET_PROJECT + printf("TRIANGLE %d (%.3f,%.3f,%.3f) (%.3f,%.3f,%.3f) (%.3f,%.3f,%.3f) %g %g\n", + index, UNPACK3(vtri_co[0]), UNPACK3(vtri_co[1]), UNPACK3(vtri_co[2]), dist_sq, nearest->dist_sq); +#endif + + if (dist_sq >= nearest->dist_sq) + return; + + /* Decode normals */ + normal_short_to_float_v3(vtri_no[0], vtri[0]->no); + normal_short_to_float_v3(vtri_no[1], vtri[1]->no); + normal_short_to_float_v3(vtri_no[2], vtri[2]->no); + + /* Solve the equations for the triangle */ + if (target_project_solve_point_tri(vtri_co, vtri_no, co, raw_hit_co, dist_sq, hit_co, hit_no)) { + update_hit(nearest, index, co, hit_co, hit_no); + } + /* Boundary edges */ + else if (tree->boundary && BLI_BITMAP_TEST(tree->boundary->looptri_has_boundary, index)) { + const BLI_bitmap *is_boundary = tree->boundary->edge_is_boundary; + int edges[3]; + + BKE_mesh_looptri_get_real_edges(tree->mesh, lt, edges); + + for (int i = 0; i < 3; i++) { + if (edges[i] >= 0 && BLI_BITMAP_TEST(is_boundary, edges[i])) { + target_project_edge(tree, index, co, nearest, edges[i]); + } + } + } +} + +/* + * Maps the point to the nearest surface, either by simple nearest, or by target normal projection. + */ +void BKE_shrinkwrap_find_nearest_surface(struct ShrinkwrapTreeData *tree, BVHTreeNearest *nearest, float co[3], int type) +{ + BVHTreeFromMesh *treeData = &tree->treeData; + + if (type == MOD_SHRINKWRAP_TARGET_PROJECT) { +#ifdef TRACE_TARGET_PROJECT + printf("====== TARGET PROJECT START ======\n"); +#endif + + BLI_bvhtree_find_nearest_ex(tree->bvh, co, nearest, mesh_looptri_target_project, tree, BVH_NEAREST_OPTIMAL_ORDER); + +#ifdef TRACE_TARGET_PROJECT + printf("====== TARGET PROJECT END: %d %g ======\n", nearest->index, nearest->dist_sq); +#endif + + if (nearest->index < 0) { + /* fallback to simple nearest */ + BLI_bvhtree_find_nearest(tree->bvh, co, nearest, treeData->nearest_callback, treeData); + } + } + else { + BLI_bvhtree_find_nearest(tree->bvh, co, nearest, treeData->nearest_callback, treeData); + } +} + +/* * Shrinkwrap moving vertexs to the nearest surface point on the target * * it builds a BVHTree from the target mesh and then performs a @@ -536,7 +1058,6 @@ static void shrinkwrap_calc_nearest_surface_point_cb_ex( ShrinkwrapCalcCBData *data = userdata; ShrinkwrapCalcData *calc = data->calc; - BVHTreeFromMesh *treeData = &data->tree->treeData; BVHTreeNearest *nearest = tls->userdata_chunk; float *co = calc->vertexCos[i]; @@ -565,12 +1086,20 @@ static void shrinkwrap_calc_nearest_surface_point_cb_ex( * If we already had an hit before.. we assume this vertex is going to have a close hit to that other vertex * so we can initiate the "nearest.dist" with the expected value to that last hit. * This will lead in pruning of the search tree. */ - if (nearest->index != -1) - nearest->dist_sq = len_squared_v3v3(tmp_co, nearest->co); + if (nearest->index != -1) { + if (calc->smd->shrinkType == MOD_SHRINKWRAP_TARGET_PROJECT) { + /* Heuristic doesn't work because of additional restrictions. */ + nearest->index = -1; + nearest->dist_sq = FLT_MAX; + } + else { + nearest->dist_sq = len_squared_v3v3(tmp_co, nearest->co); + } + } else nearest->dist_sq = FLT_MAX; - BLI_bvhtree_find_nearest(treeData->tree, tmp_co, nearest, treeData->nearest_callback, treeData); + BKE_shrinkwrap_find_nearest_surface(data->tree, nearest, tmp_co, calc->smd->shrinkType); /* Found the nearest vertex */ if (nearest->index != -1) { @@ -838,6 +1367,7 @@ void shrinkwrapModifier_deform(ShrinkwrapModifierData *smd, struct Scene *scene, switch (smd->shrinkType) { case MOD_SHRINKWRAP_NEAREST_SURFACE: + case MOD_SHRINKWRAP_TARGET_PROJECT: TIMEIT_BENCH(shrinkwrap_calc_nearest_surface_point(&calc), deform_surface); break; diff --git a/source/blender/blenlib/BLI_math_solvers.h b/source/blender/blenlib/BLI_math_solvers.h index 4352779d9ee..845d2100133 100644 --- a/source/blender/blenlib/BLI_math_solvers.h +++ b/source/blender/blenlib/BLI_math_solvers.h @@ -53,6 +53,15 @@ void BLI_svd_m3(const float m3[3][3], float r_U[3][3], float r_S[], float r_V[3] bool BLI_tridiagonal_solve(const float *a, const float *b, const float *c, const float *d, float *r_x, const int count); bool BLI_tridiagonal_solve_cyclic(const float *a, const float *b, const float *c, const float *d, float *r_x, const int count); +/* Generic 3 variable Newton's method solver. */ +typedef void (*Newton3D_DeltaFunc)(void *userdata, const float x[3], float r_delta[3]); +typedef void (*Newton3D_JacobianFunc)(void *userdata, const float x[3], float r_jacobian[3][3]); +typedef bool (*Newton3D_CorrectionFunc)(void *userdata, const float x[3], float step[3], float x_next[3]); + +bool BLI_newton3d_solve( + Newton3D_DeltaFunc func_delta, Newton3D_JacobianFunc func_jacobian, Newton3D_CorrectionFunc func_correction, void *userdata, + float epsilon, int max_iterations, bool trace, const float x_init[3], float result[3]); + #ifdef BLI_MATH_GCC_WARN_PRAGMA # pragma GCC diagnostic pop #endif diff --git a/source/blender/blenlib/intern/math_solvers.c b/source/blender/blenlib/intern/math_solvers.c index e3174d8340a..7b9727ead8e 100644 --- a/source/blender/blenlib/intern/math_solvers.c +++ b/source/blender/blenlib/intern/math_solvers.c @@ -179,3 +179,98 @@ bool BLI_tridiagonal_solve_cyclic(const float *a, const float *b, const float *c return success; } + +/** + * \brief Solve a generic f(x) = 0 equation using Newton's method. + * + * \param func_delta Callback computing the value of f(x). + * \param func_jacobian Callback computing the Jacobian matrix of the function at x. + * \param func_correction Callback for forcing the search into an arbitrary custom domain. May be NULL. + * \param userdata Data for the callbacks. + * \param epsilon Desired precision. + * \param max_iterations Limit on the iterations. + * \param max_corrections Limit on the number of times the correction callback can fire before giving up. + * \param trace Enables logging to console. + * \param x_init Initial solution vector. + * \param result Final result. + * \return true if success + */ +bool BLI_newton3d_solve( + Newton3D_DeltaFunc func_delta, Newton3D_JacobianFunc func_jacobian, Newton3D_CorrectionFunc func_correction, void *userdata, + float epsilon, int max_iterations, bool trace, const float x_init[3], float result[3]) +{ + float fdelta[3], fdeltav, next_fdeltav; + float jacobian[3][3], step[3], x[3], x_next[3]; + + epsilon *= epsilon; + + copy_v3_v3(x, x_init); + + func_delta(userdata, x, fdelta); + fdeltav = len_squared_v3(fdelta); + + if (trace) { + printf("START (%g, %g, %g) %g\n", x[0], x[1], x[2], fdeltav); + } + + for (int i = 0; i < max_iterations && fdeltav > epsilon; i++) { + /* Newton's method step. */ + func_jacobian(userdata, x, jacobian); + + if (!invert_m3(jacobian)) { + return false; + } + + mul_v3_m3v3(step, jacobian, fdelta); + sub_v3_v3v3(x_next, x, step); + + /* Custom out-of-bounds value correction. */ + if (func_correction) { + if (trace) { + printf("%3d * (%g, %g, %g)\n", i, x_next[0], x_next[1], x_next[2]); + } + + if (!func_correction(userdata, x, step, x_next)) { + return false; + } + } + + func_delta(userdata, x_next, fdelta); + next_fdeltav = len_squared_v3(fdelta); + + if (trace) { + printf("%3d ? (%g, %g, %g) %g\n", i, x_next[0], x_next[1], x_next[2], next_fdeltav); + } + + /* Line search correction. */ + while (next_fdeltav > fdeltav) { + float g0 = sqrtf(fdeltav), g1 = sqrtf(next_fdeltav); + float g01 = -g0 / len_v3(step); + float det = 2.0f * (g1 - g0 - g01); + float l = (det == 0.0f) ? 0.1f : -g01 / det; + CLAMP_MIN(l, 0.1f); + + mul_v3_fl(step, l); + sub_v3_v3v3(x_next, x, step); + + func_delta(userdata, x_next, fdelta); + next_fdeltav = len_squared_v3(fdelta); + + if (trace) { + printf("%3d . (%g, %g, %g) %g\n", i, x_next[0], x_next[1], x_next[2], next_fdeltav); + } + } + + copy_v3_v3(x, x_next); + fdeltav = next_fdeltav; + } + + bool success = (fdeltav <= epsilon); + + if (trace) { + printf("%s (%g, %g, %g) %g\n", success ? "OK " : "FAIL", x[0], x[1], x[2], fdeltav); + } + + copy_v3_v3(result, x); + return success; +} diff --git a/source/blender/depsgraph/DEG_depsgraph.h b/source/blender/depsgraph/DEG_depsgraph.h index d3900bc6350..75572a91440 100644 --- a/source/blender/depsgraph/DEG_depsgraph.h +++ b/source/blender/depsgraph/DEG_depsgraph.h @@ -80,7 +80,11 @@ enum { * to meet dependencies with such a things as curve modifier and other guys * who're using curve deform, where_on_path and so. */ - DAG_EVAL_NEED_CURVE_PATH = 1, + DAG_EVAL_NEED_CURVE_PATH = (1 << 0), + /* A shrinkwrap modifier or constraint targeting this mesh needs information + * about non-manifold boundary edges for the Target Normal Project mode. + */ + DAG_EVAL_NEED_SHRINKWRAP_BOUNDARY = (1 << 1), }; #ifdef __cplusplus diff --git a/source/blender/depsgraph/intern/builder/deg_builder_relations.cc b/source/blender/depsgraph/intern/builder/deg_builder_relations.cc index 122f4f7d317..31c16f013f9 100644 --- a/source/blender/depsgraph/intern/builder/deg_builder_relations.cc +++ b/source/blender/depsgraph/intern/builder/deg_builder_relations.cc @@ -1019,6 +1019,9 @@ void DepsgraphRelationBuilder::build_constraints(ID *id, if (track || BKE_shrinkwrap_needs_normals(scon->shrinkType, scon->shrinkMode)) { add_customdata_mask(target_key, CD_MASK_NORMAL | CD_MASK_CUSTOMLOOPNORMAL); } + if (scon->shrinkType == MOD_SHRINKWRAP_TARGET_PROJECT) { + add_special_eval_flag(&ct->tar->id, DAG_EVAL_NEED_SHRINKWRAP_BOUNDARY); + } } /* NOTE: obdata eval now doesn't necessarily depend on the diff --git a/source/blender/makesdna/DNA_mesh_types.h b/source/blender/makesdna/DNA_mesh_types.h index c591b57002a..fda18e5cbd6 100644 --- a/source/blender/makesdna/DNA_mesh_types.h +++ b/source/blender/makesdna/DNA_mesh_types.h @@ -100,6 +100,9 @@ typedef struct Mesh_Runtime { /** 'BVHCache', for 'BKE_bvhutil.c' */ struct LinkNode *bvh_cache; + /** Non-manifold boundary data for Shrinkwrap Target Project. */ + struct ShrinkwrapBoundaryData *shrinkwrap_data; + /** Set by modifier stack if only deformed from original. */ char deformed_only; /** diff --git a/source/blender/makesdna/DNA_modifier_types.h b/source/blender/makesdna/DNA_modifier_types.h index 6f8793bf0f7..4eb9f890be0 100644 --- a/source/blender/makesdna/DNA_modifier_types.h +++ b/source/blender/makesdna/DNA_modifier_types.h @@ -888,6 +888,7 @@ enum { MOD_SHRINKWRAP_NEAREST_SURFACE = 0, MOD_SHRINKWRAP_PROJECT = 1, MOD_SHRINKWRAP_NEAREST_VERTEX = 2, + MOD_SHRINKWRAP_TARGET_PROJECT = 3, }; /* Shrinkwrap->shrinkMode */ diff --git a/source/blender/makesrna/intern/rna_constraint.c b/source/blender/makesrna/intern/rna_constraint.c index e2d4093b753..c809cd7b1a4 100644 --- a/source/blender/makesrna/intern/rna_constraint.c +++ b/source/blender/makesrna/intern/rna_constraint.c @@ -2134,6 +2134,9 @@ static void rna_def_constraint_shrinkwrap(BlenderRNA *brna) "Shrink the location to the nearest target surface along a given axis"}, {MOD_SHRINKWRAP_NEAREST_VERTEX, "NEAREST_VERTEX", 0, "Nearest Vertex", "Shrink the location to the nearest target vertex"}, + {MOD_SHRINKWRAP_TARGET_PROJECT, "TARGET_PROJECT", 0, "Target Normal Project", + "Shrink the location to the nearest target surface " + "along the interpolated vertex normals of the target"}, {0, NULL, 0, NULL, NULL} }; diff --git a/source/blender/makesrna/intern/rna_modifier.c b/source/blender/makesrna/intern/rna_modifier.c index 44f14be259a..d450e95d397 100644 --- a/source/blender/makesrna/intern/rna_modifier.c +++ b/source/blender/makesrna/intern/rna_modifier.c @@ -3167,6 +3167,9 @@ static void rna_def_modifier_shrinkwrap(BlenderRNA *brna) "Shrink the mesh to the nearest target surface along a given axis"}, {MOD_SHRINKWRAP_NEAREST_VERTEX, "NEAREST_VERTEX", 0, "Nearest Vertex", "Shrink the mesh to the nearest target vertex"}, + {MOD_SHRINKWRAP_TARGET_PROJECT, "TARGET_PROJECT", 0, "Target Normal Project", + "Shrink the mesh to the nearest target surface " + "along the interpolated vertex normals of the target"}, {0, NULL, 0, NULL, NULL} }; diff --git a/source/blender/modifiers/intern/MOD_shrinkwrap.c b/source/blender/modifiers/intern/MOD_shrinkwrap.c index 5c110359f5a..477c11039db 100644 --- a/source/blender/modifiers/intern/MOD_shrinkwrap.c +++ b/source/blender/modifiers/intern/MOD_shrinkwrap.c @@ -150,10 +150,16 @@ static void updateDepsgraph(ModifierData *md, const ModifierUpdateDepsgraphConte if (smd->target != NULL) { DEG_add_object_relation(ctx->node, smd->target, DEG_OB_COMP_TRANSFORM, "Shrinkwrap Modifier"); DEG_add_object_relation_with_customdata(ctx->node, smd->target, DEG_OB_COMP_GEOMETRY, mask, "Shrinkwrap Modifier"); + if (smd->shrinkType == MOD_SHRINKWRAP_TARGET_PROJECT) { + DEG_add_special_eval_flag(ctx->node, &smd->target->id, DAG_EVAL_NEED_SHRINKWRAP_BOUNDARY); + } } if (smd->auxTarget != NULL) { DEG_add_object_relation(ctx->node, smd->auxTarget, DEG_OB_COMP_TRANSFORM, "Shrinkwrap Modifier"); DEG_add_object_relation_with_customdata(ctx->node, smd->auxTarget, DEG_OB_COMP_GEOMETRY, mask, "Shrinkwrap Modifier"); + if (smd->shrinkType == MOD_SHRINKWRAP_TARGET_PROJECT) { + DEG_add_special_eval_flag(ctx->node, &smd->auxTarget->id, DAG_EVAL_NEED_SHRINKWRAP_BOUNDARY); + } } DEG_add_object_relation(ctx->node, ctx->object, DEG_OB_COMP_TRANSFORM, "Shrinkwrap Modifier"); } |