From d4f77c22661c316bdaada2ff053f1484412ccb83 Mon Sep 17 00:00:00 2001 From: Hans Goudey Date: Thu, 6 Oct 2022 15:15:49 -0500 Subject: Cleanup: Move shrinkwrap.c to C++ --- source/blender/blenkernel/CMakeLists.txt | 2 +- source/blender/blenkernel/intern/shrinkwrap.c | 1591 ------------------------ source/blender/blenkernel/intern/shrinkwrap.cc | 1590 +++++++++++++++++++++++ 3 files changed, 1591 insertions(+), 1592 deletions(-) delete mode 100644 source/blender/blenkernel/intern/shrinkwrap.c create mode 100644 source/blender/blenkernel/intern/shrinkwrap.cc (limited to 'source') diff --git a/source/blender/blenkernel/CMakeLists.txt b/source/blender/blenkernel/CMakeLists.txt index 192acbf8338..a58fc22b880 100644 --- a/source/blender/blenkernel/CMakeLists.txt +++ b/source/blender/blenkernel/CMakeLists.txt @@ -259,7 +259,7 @@ set(SRC intern/scene.cc intern/screen.c intern/shader_fx.c - intern/shrinkwrap.c + intern/shrinkwrap.cc intern/simulation.cc intern/softbody.c intern/sound.c diff --git a/source/blender/blenkernel/intern/shrinkwrap.c b/source/blender/blenkernel/intern/shrinkwrap.c deleted file mode 100644 index f94836551bb..00000000000 --- a/source/blender/blenkernel/intern/shrinkwrap.c +++ /dev/null @@ -1,1591 +0,0 @@ -/* SPDX-License-Identifier: GPL-2.0-or-later - * Copyright Blender Foundation. All rights reserved. */ - -/** \file - * \ingroup bke - */ - -#include -#include -#include -#include -#include -#include - -#include "DNA_gpencil_modifier_types.h" -#include "DNA_mesh_types.h" -#include "DNA_meshdata_types.h" -#include "DNA_modifier_types.h" -#include "DNA_object_types.h" - -#include "BLI_math.h" -#include "BLI_math_solvers.h" -#include "BLI_task.h" -#include "BLI_utildefines.h" - -#include "BKE_DerivedMesh.h" -#include "BKE_cdderivedmesh.h" -#include "BKE_context.h" -#include "BKE_lattice.h" -#include "BKE_lib_id.h" -#include "BKE_modifier.h" -#include "BKE_shrinkwrap.h" - -#include "BKE_deform.h" -#include "BKE_editmesh.h" -#include "BKE_mesh.h" /* for OMP limits. */ -#include "BKE_mesh_runtime.h" -#include "BKE_mesh_wrapper.h" -#include "BKE_subsurf.h" - -#include "DEG_depsgraph_query.h" - -#include "MEM_guardedalloc.h" - -#include "BLI_strict_flags.h" - -/* for timing... */ -#if 0 -# include "PIL_time_utildefines.h" -#else -# define TIMEIT_BENCH(expr, id) (expr) -#endif - -/* Util macros */ -#define OUT_OF_MEMORY() ((void)printf("Shrinkwrap: Out of memory\n")) - -typedef struct ShrinkwrapCalcData { - ShrinkwrapModifierData *smd; /* shrinkwrap modifier data */ - - struct Object *ob; /* object we are applying shrinkwrap to */ - - struct MVert *vert; /* Array of verts being projected. */ - const float (*vert_normals)[3]; - float (*vertexCos)[3]; /* vertexs being shrinkwraped */ - int numVerts; - - const struct MDeformVert *dvert; /* Pointer to mdeform array */ - int vgroup; /* Vertex group num */ - bool invert_vgroup; /* invert vertex group influence */ - - struct Mesh *target; /* mesh we are shrinking to */ - struct SpaceTransform local2target; /* transform to move between local and target space */ - struct ShrinkwrapTreeData *tree; /* mesh BVH tree data */ - - struct Object *aux_target; - - float keepDist; /* Distance to keep above target surface (units are in local space) */ -} ShrinkwrapCalcData; - -typedef struct ShrinkwrapCalcCBData { - ShrinkwrapCalcData *calc; - - ShrinkwrapTreeData *tree; - ShrinkwrapTreeData *aux_tree; - - float *proj_axis; - SpaceTransform *local2aux; -} ShrinkwrapCalcCBData; - -bool BKE_shrinkwrap_needs_normals(int shrinkType, int shrinkMode) -{ - return (shrinkType == MOD_SHRINKWRAP_TARGET_PROJECT) || - (shrinkType != MOD_SHRINKWRAP_NEAREST_VERTEX && - shrinkMode == MOD_SHRINKWRAP_ABOVE_SURFACE); -} - -bool BKE_shrinkwrap_init_tree( - ShrinkwrapTreeData *data, Mesh *mesh, int shrinkType, int shrinkMode, bool force_normals) -{ - memset(data, 0, sizeof(*data)); - - if (mesh == NULL) { - return false; - } - - /* We could create a BVH tree from the edit mesh, - * however accessing normals from the face/loop normals gets more involved. - * Convert mesh data since this isn't typically used in edit-mode. */ - BKE_mesh_wrapper_ensure_mdata(mesh); - - if (mesh->totvert <= 0) { - return false; - } - - data->mesh = mesh; - data->polys = BKE_mesh_polys(mesh); - - if (shrinkType == MOD_SHRINKWRAP_NEAREST_VERTEX) { - data->bvh = BKE_bvhtree_from_mesh_get(&data->treeData, mesh, BVHTREE_FROM_VERTS, 2); - - return data->bvh != NULL; - } - - if (mesh->totpoly <= 0) { - return false; - } - - data->bvh = BKE_bvhtree_from_mesh_get(&data->treeData, mesh, BVHTREE_FROM_LOOPTRI, 4); - - if (data->bvh == NULL) { - return false; - } - - if (force_normals || BKE_shrinkwrap_needs_normals(shrinkType, shrinkMode)) { - data->pnors = BKE_mesh_poly_normals_ensure(mesh); - if ((mesh->flag & ME_AUTOSMOOTH) != 0) { - data->clnors = CustomData_get_layer(&mesh->ldata, CD_NORMAL); - } - } - - if (shrinkType == MOD_SHRINKWRAP_TARGET_PROJECT) { - data->boundary = mesh->runtime.shrinkwrap_data; - } - - return true; -} - -void BKE_shrinkwrap_free_tree(ShrinkwrapTreeData *data) -{ - free_bvhtree_from_mesh(&data->treeData); -} - -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 MVert *mvert = BKE_mesh_verts(mesh); - const MEdge *medge = BKE_mesh_edges(mesh); - const MLoop *mloop = BKE_mesh_loops(mesh); - - /* 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++) { - uint 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"); - uint 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; - } - } - - uint 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. */ - const float(*vert_normals)[3] = BKE_mesh_vertex_normals_ensure(mesh); - for (int i = 0; i < mesh->totvert; i++) { - int bidx = vert_boundary_id[i]; - - if (bidx >= 0) { - ShrinkwrapBoundaryVertData *vdata = &boundary_verts[bidx]; - float tmp[3]; - - normalize_v3(vdata->direction); - - cross_v3_v3v3(tmp, vert_normals[i], vdata->direction); - cross_v3_v3v3(vdata->normal_plane, tmp, vert_normals[i]); - 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); -} - -/** - * Shrink-wrap to the nearest vertex - * - * it builds a #BVHTree of vertices we can attach to and then - * for each vertex performs a nearest vertex search on the tree - */ -static void shrinkwrap_calc_nearest_vertex_cb_ex(void *__restrict userdata, - const int i, - const TaskParallelTLS *__restrict tls) -{ - ShrinkwrapCalcCBData *data = userdata; - - ShrinkwrapCalcData *calc = data->calc; - BVHTreeFromMesh *treeData = &data->tree->treeData; - BVHTreeNearest *nearest = tls->userdata_chunk; - - float *co = calc->vertexCos[i]; - float tmp_co[3]; - float weight = BKE_defvert_array_find_weight_safe(calc->dvert, i, calc->vgroup); - - if (calc->invert_vgroup) { - weight = 1.0f - weight; - } - - if (weight == 0.0f) { - return; - } - - /* Convert the vertex to tree coordinates */ - if (calc->vert) { - copy_v3_v3(tmp_co, calc->vert[i].co); - } - else { - copy_v3_v3(tmp_co, co); - } - BLI_space_transform_apply(&calc->local2target, tmp_co); - - /* Use local proximity heuristics (to reduce the nearest search) - * - * 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); - } - else { - nearest->dist_sq = FLT_MAX; - } - - BLI_bvhtree_find_nearest(treeData->tree, tmp_co, nearest, treeData->nearest_callback, treeData); - - /* Found the nearest vertex */ - if (nearest->index != -1) { - /* Adjusting the vertex weight, - * so that after interpolating it keeps a certain distance from the nearest position */ - if (nearest->dist_sq > FLT_EPSILON) { - const float dist = sqrtf(nearest->dist_sq); - weight *= (dist - calc->keepDist) / dist; - } - - /* Convert the coordinates back to mesh coordinates */ - copy_v3_v3(tmp_co, nearest->co); - BLI_space_transform_invert(&calc->local2target, tmp_co); - - interp_v3_v3v3(co, co, tmp_co, weight); /* linear interpolation */ - } -} - -static void shrinkwrap_calc_nearest_vertex(ShrinkwrapCalcData *calc) -{ - BVHTreeNearest nearest = NULL_BVHTreeNearest; - - /* Setup nearest */ - nearest.index = -1; - nearest.dist_sq = FLT_MAX; - - ShrinkwrapCalcCBData data = { - .calc = calc, - .tree = calc->tree, - }; - TaskParallelSettings settings; - BLI_parallel_range_settings_defaults(&settings); - settings.use_threading = (calc->numVerts > BKE_MESH_OMP_LIMIT); - settings.userdata_chunk = &nearest; - settings.userdata_chunk_size = sizeof(nearest); - BLI_task_parallel_range( - 0, calc->numVerts, &data, shrinkwrap_calc_nearest_vertex_cb_ex, &settings); -} - -bool BKE_shrinkwrap_project_normal(char options, - const float vert[3], - const float dir[3], - const float ray_radius, - const SpaceTransform *transf, - ShrinkwrapTreeData *tree, - BVHTreeRayHit *hit) -{ - /* don't use this because this dist value could be incompatible - * this value used by the callback for comparing previous/new dist values. - * also, at the moment there is no need to have a corrected 'dist' value */ - // #define USE_DIST_CORRECT - - float tmp_co[3], tmp_no[3]; - const float *co, *no; - BVHTreeRayHit hit_tmp; - - /* Copy from hit (we need to convert hit rays from one space coordinates to the other */ - memcpy(&hit_tmp, hit, sizeof(hit_tmp)); - - /* Apply space transform (TODO readjust dist) */ - if (transf) { - copy_v3_v3(tmp_co, vert); - BLI_space_transform_apply(transf, tmp_co); - co = tmp_co; - - copy_v3_v3(tmp_no, dir); - BLI_space_transform_apply_normal(transf, tmp_no); - no = tmp_no; - -#ifdef USE_DIST_CORRECT - hit_tmp.dist *= mat4_to_scale(((SpaceTransform *)transf)->local2target); -#endif - } - else { - co = vert; - no = dir; - } - - hit_tmp.index = -1; - - BLI_bvhtree_ray_cast( - tree->bvh, co, no, ray_radius, &hit_tmp, tree->treeData.raycast_callback, &tree->treeData); - - if (hit_tmp.index != -1) { - /* invert the normal first so face culling works on rotated objects */ - if (transf) { - BLI_space_transform_invert_normal(transf, hit_tmp.no); - } - - if (options & MOD_SHRINKWRAP_CULL_TARGET_MASK) { - /* Apply back-face. */ - const float dot = dot_v3v3(dir, hit_tmp.no); - if (((options & MOD_SHRINKWRAP_CULL_TARGET_FRONTFACE) && dot <= 0.0f) || - ((options & MOD_SHRINKWRAP_CULL_TARGET_BACKFACE) && dot >= 0.0f)) { - return false; /* Ignore hit */ - } - } - - if (transf) { - /* Inverting space transform (TODO: make coherent with the initial dist readjust). */ - BLI_space_transform_invert(transf, hit_tmp.co); -#ifdef USE_DIST_CORRECT - hit_tmp.dist = len_v3v3(vert, hit_tmp.co); -#endif - } - - BLI_assert(hit_tmp.dist <= hit->dist); - - memcpy(hit, &hit_tmp, sizeof(hit_tmp)); - return true; - } - return false; -} - -static void shrinkwrap_calc_normal_projection_cb_ex(void *__restrict userdata, - const int i, - const TaskParallelTLS *__restrict tls) -{ - ShrinkwrapCalcCBData *data = userdata; - - ShrinkwrapCalcData *calc = data->calc; - ShrinkwrapTreeData *tree = data->tree; - ShrinkwrapTreeData *aux_tree = data->aux_tree; - - float *proj_axis = data->proj_axis; - SpaceTransform *local2aux = data->local2aux; - - BVHTreeRayHit *hit = tls->userdata_chunk; - - const float proj_limit_squared = calc->smd->projLimit * calc->smd->projLimit; - float *co = calc->vertexCos[i]; - float tmp_co[3], tmp_no[3]; - float weight = BKE_defvert_array_find_weight_safe(calc->dvert, i, calc->vgroup); - - if (calc->invert_vgroup) { - weight = 1.0f - weight; - } - - if (weight == 0.0f) { - return; - } - - if (calc->vert != NULL && calc->smd->projAxis == MOD_SHRINKWRAP_PROJECT_OVER_NORMAL) { - /* calc->vert contains verts from evaluated mesh. */ - /* These coordinates are deformed by vertexCos only for normal projection - * (to get correct normals) for other cases calc->verts contains undeformed coordinates and - * vertexCos should be used */ - copy_v3_v3(tmp_co, calc->vert[i].co); - copy_v3_v3(tmp_no, calc->vert_normals[i]); - } - else { - copy_v3_v3(tmp_co, co); - copy_v3_v3(tmp_no, proj_axis); - } - - hit->index = -1; - - /* TODO: we should use FLT_MAX here, but sweepsphere code isn't prepared for that */ - hit->dist = BVH_RAYCAST_DIST_MAX; - - bool is_aux = false; - - /* Project over positive direction of axis */ - if (calc->smd->shrinkOpts & MOD_SHRINKWRAP_PROJECT_ALLOW_POS_DIR) { - if (aux_tree) { - if (BKE_shrinkwrap_project_normal(0, tmp_co, tmp_no, 0.0, local2aux, aux_tree, hit)) { - is_aux = true; - } - } - - if (BKE_shrinkwrap_project_normal( - calc->smd->shrinkOpts, tmp_co, tmp_no, 0.0, &calc->local2target, tree, hit)) { - is_aux = false; - } - } - - /* Project over negative direction of axis */ - if (calc->smd->shrinkOpts & MOD_SHRINKWRAP_PROJECT_ALLOW_NEG_DIR) { - float inv_no[3]; - negate_v3_v3(inv_no, tmp_no); - - char options = calc->smd->shrinkOpts; - - if ((options & MOD_SHRINKWRAP_INVERT_CULL_TARGET) && - (options & MOD_SHRINKWRAP_CULL_TARGET_MASK)) { - options ^= MOD_SHRINKWRAP_CULL_TARGET_MASK; - } - - if (aux_tree) { - if (BKE_shrinkwrap_project_normal(0, tmp_co, inv_no, 0.0, local2aux, aux_tree, hit)) { - is_aux = true; - } - } - - if (BKE_shrinkwrap_project_normal( - options, tmp_co, inv_no, 0.0, &calc->local2target, tree, hit)) { - is_aux = false; - } - } - - /* don't set the initial dist (which is more efficient), - * because its calculated in the targets space, we want the dist in our own space */ - if (proj_limit_squared != 0.0f) { - if (hit->index != -1 && len_squared_v3v3(hit->co, co) > proj_limit_squared) { - hit->index = -1; - } - } - - if (hit->index != -1) { - if (is_aux) { - BKE_shrinkwrap_snap_point_to_surface(aux_tree, - local2aux, - calc->smd->shrinkMode, - hit->index, - hit->co, - hit->no, - calc->keepDist, - tmp_co, - hit->co); - } - else { - BKE_shrinkwrap_snap_point_to_surface(tree, - &calc->local2target, - calc->smd->shrinkMode, - hit->index, - hit->co, - hit->no, - calc->keepDist, - tmp_co, - hit->co); - } - - interp_v3_v3v3(co, co, hit->co, weight); - } -} - -static void shrinkwrap_calc_normal_projection(ShrinkwrapCalcData *calc) -{ - /* Options about projection direction */ - float proj_axis[3] = {0.0f, 0.0f, 0.0f}; - - /* Ray-cast and tree stuff. */ - - /** \note 'hit.dist' is kept in the targets space, this is only used - * for finding the best hit, to get the real dist, - * measure the len_v3v3() from the input coord to hit.co */ - BVHTreeRayHit hit; - - /* auxiliary target */ - Mesh *auxMesh = NULL; - ShrinkwrapTreeData *aux_tree = NULL; - ShrinkwrapTreeData aux_tree_stack; - SpaceTransform local2aux; - - /* If the user doesn't allows to project in any direction of projection axis - * then there's nothing todo. */ - if ((calc->smd->shrinkOpts & - (MOD_SHRINKWRAP_PROJECT_ALLOW_POS_DIR | MOD_SHRINKWRAP_PROJECT_ALLOW_NEG_DIR)) == 0) { - return; - } - - /* Prepare data to retrieve the direction in which we should project each vertex */ - if (calc->smd->projAxis == MOD_SHRINKWRAP_PROJECT_OVER_NORMAL) { - if (calc->vert == NULL) { - return; - } - } - else { - /* The code supports any axis that is a combination of X,Y,Z - * although currently UI only allows to set the 3 different axis */ - if (calc->smd->projAxis & MOD_SHRINKWRAP_PROJECT_OVER_X_AXIS) { - proj_axis[0] = 1.0f; - } - if (calc->smd->projAxis & MOD_SHRINKWRAP_PROJECT_OVER_Y_AXIS) { - proj_axis[1] = 1.0f; - } - if (calc->smd->projAxis & MOD_SHRINKWRAP_PROJECT_OVER_Z_AXIS) { - proj_axis[2] = 1.0f; - } - - normalize_v3(proj_axis); - - /* Invalid projection direction */ - if (len_squared_v3(proj_axis) < FLT_EPSILON) { - return; - } - } - - if (calc->aux_target) { - auxMesh = BKE_modifier_get_evaluated_mesh_from_evaluated_object(calc->aux_target); - if (!auxMesh) { - return; - } - BLI_SPACE_TRANSFORM_SETUP(&local2aux, calc->ob, calc->aux_target); - } - - if (BKE_shrinkwrap_init_tree( - &aux_tree_stack, auxMesh, calc->smd->shrinkType, calc->smd->shrinkMode, false)) { - aux_tree = &aux_tree_stack; - } - - /* After successfully build the trees, start projection vertices. */ - ShrinkwrapCalcCBData data = { - .calc = calc, - .tree = calc->tree, - .aux_tree = aux_tree, - .proj_axis = proj_axis, - .local2aux = &local2aux, - }; - TaskParallelSettings settings; - BLI_parallel_range_settings_defaults(&settings); - settings.use_threading = (calc->numVerts > BKE_MESH_OMP_LIMIT); - settings.userdata_chunk = &hit; - settings.userdata_chunk_size = sizeof(hit); - BLI_task_parallel_range( - 0, calc->numVerts, &data, shrinkwrap_calc_normal_projection_cb_ex, &settings); - - /* free data structures */ - if (aux_tree) { - BKE_shrinkwrap_free_tree(aux_tree); - } -} - -/* - * 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; - - const 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-5f; - /* Dot product threshold for checking if step is 'clearly' pointing outside. */ - const float dir_epsilon = 0.5f; - bool fixed = false, locked = false; - - /* The barycentric coordinate domain is a triangle bounded by - * the X and Y axes, plus the x+y=1 diagonal. First, clamp the - * movement against the diagonal. Note that step is subtracted. */ - float sum = x[0] + x[1]; - float sstep = -(step[0] + step[1]); - - if (sum + sstep > 1.0f) { - float ldist = 1.0f - sum; - - /* If already at the boundary, slide along it. */ - if (ldist < epsilon * (float)M_SQRT2) { - float step_len = len_v2(step); - - /* Abort if the solution is clearly outside the domain. */ - if (step_len > epsilon && sstep > step_len * dir_epsilon * (float)M_SQRT2) { - return false; - } - - /* Project the new position onto the diagonal. */ - add_v2_fl(step, (sum + sstep - 1.0f) * 0.5f); - fixed = locked = true; - } - else { - /* Scale a significant step down to arrive at the boundary. */ - mul_v3_fl(step, ldist / sstep); - fixed = true; - } - } - - /* Weight 0 and 1 boundary checks - along axis. */ - for (int i = 0; i < 2; i++) { - if (step[i] > x[i]) { - /* If already at the boundary, slide along it. */ - if (x[i] < epsilon) { - float step_len = len_v2(step); - - /* Abort if the solution is clearly outside the domain. */ - if (step_len > epsilon && (locked || step[i] > step_len * dir_epsilon)) { - return false; - } - - /* Reset precision errors to stay at the boundary. */ - step[i] = x[i]; - fixed = true; - } - else { - /* Scale a significant step down to arrive at the boundary. */ - mul_v3_fl(step, x[i] / step[i]); - 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 magnitude_estimate = dist + len_manhattan_v3(vtri_co[0]) + len_manhattan_v3(vtri_co[1]) + - len_manhattan_v3(vtri_co[2]); - float epsilon = magnitude_estimate * 1.0e-6f; - - /* 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 = &data->edge[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]; - copy_v3_v3(vedge_no[0], data->vert_normals[edge->v1]); - copy_v3_v3(vedge_no[1], data->vert_normals[edge->v2]); - - 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 */ - copy_v3_v3(vtri_no[0], tree->treeData.vert_normals[loop[0]->v]); - copy_v3_v3(vtri_no[1], tree->treeData.vert_normals[loop[1]->v]); - copy_v3_v3(vtri_no[2], tree->treeData.vert_normals[loop[2]->v]); - - /* 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]); - } - } - } -} - -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("\n====== 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\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 - * NN matches for each vertex - */ -static void shrinkwrap_calc_nearest_surface_point_cb_ex(void *__restrict userdata, - const int i, - const TaskParallelTLS *__restrict tls) -{ - ShrinkwrapCalcCBData *data = userdata; - - ShrinkwrapCalcData *calc = data->calc; - BVHTreeNearest *nearest = tls->userdata_chunk; - - float *co = calc->vertexCos[i]; - float tmp_co[3]; - float weight = BKE_defvert_array_find_weight_safe(calc->dvert, i, calc->vgroup); - - if (calc->invert_vgroup) { - weight = 1.0f - weight; - } - - if (weight == 0.0f) { - return; - } - - /* Convert the vertex to tree coordinates */ - if (calc->vert) { - copy_v3_v3(tmp_co, calc->vert[i].co); - } - else { - copy_v3_v3(tmp_co, co); - } - BLI_space_transform_apply(&calc->local2target, tmp_co); - - /* Use local proximity heuristics (to reduce the nearest search) - * - * 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) { - 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; - } - - BKE_shrinkwrap_find_nearest_surface(data->tree, nearest, tmp_co, calc->smd->shrinkType); - - /* Found the nearest vertex */ - if (nearest->index != -1) { - BKE_shrinkwrap_snap_point_to_surface(data->tree, - NULL, - calc->smd->shrinkMode, - nearest->index, - nearest->co, - nearest->no, - calc->keepDist, - tmp_co, - tmp_co); - - /* Convert the coordinates back to mesh coordinates */ - BLI_space_transform_invert(&calc->local2target, tmp_co); - interp_v3_v3v3(co, co, tmp_co, weight); /* linear interpolation */ - } -} - -void BKE_shrinkwrap_compute_smooth_normal(const struct ShrinkwrapTreeData *tree, - const struct SpaceTransform *transform, - int looptri_idx, - const float hit_co[3], - const float hit_no[3], - float r_no[3]) -{ - const BVHTreeFromMesh *treeData = &tree->treeData; - const MLoopTri *tri = &treeData->looptri[looptri_idx]; - const float(*vert_normals)[3] = tree->treeData.vert_normals; - - /* Interpolate smooth normals if enabled. */ - if ((tree->polys[tri->poly].flag & ME_SMOOTH) != 0) { - const uint32_t vert_indices[3] = {treeData->loop[tri->tri[0]].v, - treeData->loop[tri->tri[1]].v, - treeData->loop[tri->tri[2]].v}; - float w[3], no[3][3], tmp_co[3]; - - /* Custom and auto smooth split normals. */ - if (tree->clnors) { - copy_v3_v3(no[0], tree->clnors[tri->tri[0]]); - copy_v3_v3(no[1], tree->clnors[tri->tri[1]]); - copy_v3_v3(no[2], tree->clnors[tri->tri[2]]); - } - /* Ordinary vertex normals. */ - else { - copy_v3_v3(no[0], vert_normals[vert_indices[0]]); - copy_v3_v3(no[1], vert_normals[vert_indices[1]]); - copy_v3_v3(no[2], vert_normals[vert_indices[2]]); - } - - /* Barycentric weights from hit point. */ - copy_v3_v3(tmp_co, hit_co); - - if (transform) { - BLI_space_transform_apply(transform, tmp_co); - } - - interp_weights_tri_v3(w, - treeData->vert[vert_indices[0]].co, - treeData->vert[vert_indices[1]].co, - treeData->vert[vert_indices[2]].co, - tmp_co); - - /* Interpolate using weights. */ - interp_v3_v3v3v3(r_no, no[0], no[1], no[2], w); - - if (transform) { - BLI_space_transform_invert_normal(transform, r_no); - } - else { - normalize_v3(r_no); - } - } - /* Use the polygon normal if flat. */ - else if (tree->pnors != NULL) { - copy_v3_v3(r_no, tree->pnors[tri->poly]); - } - /* Finally fallback to the looptri normal. */ - else { - copy_v3_v3(r_no, hit_no); - } -} - -/* Helper for MOD_SHRINKWRAP_INSIDE, MOD_SHRINKWRAP_OUTSIDE and MOD_SHRINKWRAP_OUTSIDE_SURFACE. */ -static void shrinkwrap_snap_with_side(float r_point_co[3], - const float point_co[3], - const float hit_co[3], - const float hit_no[3], - float goal_dist, - float forcesign, - bool forcesnap) -{ - float delta[3]; - sub_v3_v3v3(delta, point_co, hit_co); - - float dist = len_v3(delta); - - /* If exactly on the surface, push out along normal */ - if (dist < FLT_EPSILON) { - if (forcesnap || goal_dist > 0) { - madd_v3_v3v3fl(r_point_co, hit_co, hit_no, goal_dist * forcesign); - } - else { - copy_v3_v3(r_point_co, hit_co); - } - } - /* Move to the correct side if needed */ - else { - float dsign = signf(dot_v3v3(delta, hit_no)); - - if (forcesign == 0.0f) { - forcesign = dsign; - } - - /* If on the wrong side or too close, move to correct */ - if (forcesnap || dsign * dist * forcesign < goal_dist) { - mul_v3_fl(delta, dsign / dist); - - /* At very small distance, blend in the hit normal to stabilize math. */ - float dist_epsilon = (fabsf(goal_dist) + len_manhattan_v3(hit_co)) * 1e-4f; - - if (dist < dist_epsilon) { -#ifdef TRACE_TARGET_PROJECT - printf("zero_factor %g = %g / %g\n", dist / dist_epsilon, dist, dist_epsilon); -#endif - - interp_v3_v3v3(delta, hit_no, delta, dist / dist_epsilon); - } - - madd_v3_v3v3fl(r_point_co, hit_co, delta, goal_dist * forcesign); - } - else { - copy_v3_v3(r_point_co, point_co); - } - } -} - -void BKE_shrinkwrap_snap_point_to_surface(const struct ShrinkwrapTreeData *tree, - const struct SpaceTransform *transform, - int mode, - int hit_idx, - const float hit_co[3], - const float hit_no[3], - float goal_dist, - const float point_co[3], - float r_point_co[3]) -{ - float tmp[3]; - - switch (mode) { - /* Offsets along the line between point_co and hit_co. */ - case MOD_SHRINKWRAP_ON_SURFACE: - if (goal_dist != 0) { - shrinkwrap_snap_with_side(r_point_co, point_co, hit_co, hit_no, goal_dist, 0, true); - } - else { - copy_v3_v3(r_point_co, hit_co); - } - break; - - case MOD_SHRINKWRAP_INSIDE: - shrinkwrap_snap_with_side(r_point_co, point_co, hit_co, hit_no, goal_dist, -1, false); - break; - - case MOD_SHRINKWRAP_OUTSIDE: - shrinkwrap_snap_with_side(r_point_co, point_co, hit_co, hit_no, goal_dist, +1, false); - break; - - case MOD_SHRINKWRAP_OUTSIDE_SURFACE: - if (goal_dist != 0) { - shrinkwrap_snap_with_side(r_point_co, point_co, hit_co, hit_no, goal_dist, +1, true); - } - else { - copy_v3_v3(r_point_co, hit_co); - } - break; - - /* Offsets along the normal */ - case MOD_SHRINKWRAP_ABOVE_SURFACE: - if (goal_dist != 0) { - BKE_shrinkwrap_compute_smooth_normal(tree, transform, hit_idx, hit_co, hit_no, tmp); - madd_v3_v3v3fl(r_point_co, hit_co, tmp, goal_dist); - } - else { - copy_v3_v3(r_point_co, hit_co); - } - break; - - default: - printf("Unknown Shrinkwrap surface snap mode: %d\n", mode); - copy_v3_v3(r_point_co, hit_co); - } -} - -static void shrinkwrap_calc_nearest_surface_point(ShrinkwrapCalcData *calc) -{ - BVHTreeNearest nearest = NULL_BVHTreeNearest; - - /* Setup nearest */ - nearest.index = -1; - nearest.dist_sq = FLT_MAX; - - /* Find the nearest vertex */ - ShrinkwrapCalcCBData data = { - .calc = calc, - .tree = calc->tree, - }; - TaskParallelSettings settings; - BLI_parallel_range_settings_defaults(&settings); - settings.use_threading = (calc->numVerts > BKE_MESH_OMP_LIMIT); - settings.userdata_chunk = &nearest; - settings.userdata_chunk_size = sizeof(nearest); - BLI_task_parallel_range( - 0, calc->numVerts, &data, shrinkwrap_calc_nearest_surface_point_cb_ex, &settings); -} - -void shrinkwrapModifier_deform(ShrinkwrapModifierData *smd, - const ModifierEvalContext *ctx, - struct Scene *scene, - Object *ob, - Mesh *mesh, - const MDeformVert *dvert, - const int defgrp_index, - float (*vertexCos)[3], - int numVerts) -{ - - DerivedMesh *ss_mesh = NULL; - ShrinkwrapCalcData calc = NULL_ShrinkwrapCalcData; - - /* remove loop dependencies on derived meshes (TODO should this be done elsewhere?) */ - if (smd->target == ob) { - smd->target = NULL; - } - if (smd->auxTarget == ob) { - smd->auxTarget = NULL; - } - - /* Configure Shrinkwrap calc data */ - calc.smd = smd; - calc.ob = ob; - calc.numVerts = numVerts; - calc.vertexCos = vertexCos; - calc.dvert = dvert; - calc.vgroup = defgrp_index; - calc.invert_vgroup = (smd->shrinkOpts & MOD_SHRINKWRAP_INVERT_VGROUP) != 0; - - if (smd->target != NULL) { - Object *ob_target = DEG_get_evaluated_object(ctx->depsgraph, smd->target); - calc.target = BKE_modifier_get_evaluated_mesh_from_evaluated_object(ob_target); - - /* TODO: there might be several "bugs" with non-uniform scales matrices - * because it will no longer be nearest surface, not sphere projection - * because space has been deformed */ - BLI_SPACE_TRANSFORM_SETUP(&calc.local2target, ob, ob_target); - - /* TODO: smd->keepDist is in global units.. must change to local */ - calc.keepDist = smd->keepDist; - } - calc.aux_target = DEG_get_evaluated_object(ctx->depsgraph, smd->auxTarget); - - if (mesh != NULL && smd->shrinkType == MOD_SHRINKWRAP_PROJECT) { - /* Setup arrays to get vertexs positions, normals and deform weights */ - calc.vert = BKE_mesh_verts_for_write(mesh); - calc.vert_normals = BKE_mesh_vertex_normals_ensure(mesh); - - /* Using vertexs positions/normals as if a subsurface was applied */ - if (smd->subsurfLevels) { - SubsurfModifierData ssmd = {{NULL}}; - ssmd.subdivType = ME_CC_SUBSURF; /* catmull clark */ - ssmd.levels = smd->subsurfLevels; /* levels */ - - /* TODO: to be moved to Mesh once we are done with changes in subsurf code. */ - DerivedMesh *dm = CDDM_from_mesh(mesh); - - ss_mesh = subsurf_make_derived_from_derived( - dm, &ssmd, scene, NULL, (ob->mode & OB_MODE_EDIT) ? SUBSURF_IN_EDIT_MODE : 0); - - if (ss_mesh) { - calc.vert = ss_mesh->getVertDataArray(ss_mesh, CD_MVERT); - if (calc.vert) { - /* TRICKY: this code assumes subsurface will have the transformed original vertices - * in their original order at the end of the vert array. */ - calc.vert = calc.vert + ss_mesh->getNumVerts(ss_mesh) - dm->getNumVerts(dm); - } - } - - /* Just to make sure we are not leaving any memory behind */ - BLI_assert(ssmd.emCache == NULL); - BLI_assert(ssmd.mCache == NULL); - - dm->release(dm); - } - } - - /* Projecting target defined - lets work! */ - ShrinkwrapTreeData tree; - - if (BKE_shrinkwrap_init_tree(&tree, calc.target, smd->shrinkType, smd->shrinkMode, false)) { - calc.tree = &tree; - - switch (smd->shrinkType) { - case MOD_SHRINKWRAP_NEAREST_SURFACE: - case MOD_SHRINKWRAP_TARGET_PROJECT: - TIMEIT_BENCH(shrinkwrap_calc_nearest_surface_point(&calc), deform_surface); - break; - - case MOD_SHRINKWRAP_PROJECT: - TIMEIT_BENCH(shrinkwrap_calc_normal_projection(&calc), deform_project); - break; - - case MOD_SHRINKWRAP_NEAREST_VERTEX: - TIMEIT_BENCH(shrinkwrap_calc_nearest_vertex(&calc), deform_vertex); - break; - } - - BKE_shrinkwrap_free_tree(&tree); - } - - /* free memory */ - if (ss_mesh) { - ss_mesh->release(ss_mesh); - } -} - -void shrinkwrapGpencilModifier_deform(ShrinkwrapGpencilModifierData *mmd, - Object *ob, - MDeformVert *dvert, - const int defgrp_index, - float (*vertexCos)[3], - int numVerts) -{ - - ShrinkwrapCalcData calc = NULL_ShrinkwrapCalcData; - /* Convert gpencil struct to use the same struct and function used with meshes. */ - ShrinkwrapModifierData smd; - smd.target = mmd->target; - smd.auxTarget = mmd->aux_target; - smd.keepDist = mmd->keep_dist; - smd.shrinkType = mmd->shrink_type; - smd.shrinkOpts = mmd->shrink_opts; - smd.shrinkMode = mmd->shrink_mode; - smd.projLimit = mmd->proj_limit; - smd.projAxis = mmd->proj_axis; - - /* Configure Shrinkwrap calc data. */ - calc.smd = &smd; - calc.ob = ob; - calc.numVerts = numVerts; - calc.vertexCos = vertexCos; - calc.dvert = dvert; - calc.vgroup = defgrp_index; - calc.invert_vgroup = (mmd->flag & GP_SHRINKWRAP_INVERT_VGROUP) != 0; - - BLI_SPACE_TRANSFORM_SETUP(&calc.local2target, ob, mmd->target); - calc.keepDist = mmd->keep_dist; - calc.tree = mmd->cache_data; - - switch (mmd->shrink_type) { - case MOD_SHRINKWRAP_NEAREST_SURFACE: - case MOD_SHRINKWRAP_TARGET_PROJECT: - TIMEIT_BENCH(shrinkwrap_calc_nearest_surface_point(&calc), gpdeform_surface); - break; - - case MOD_SHRINKWRAP_PROJECT: - TIMEIT_BENCH(shrinkwrap_calc_normal_projection(&calc), gpdeform_project); - break; - - case MOD_SHRINKWRAP_NEAREST_VERTEX: - TIMEIT_BENCH(shrinkwrap_calc_nearest_vertex(&calc), gpdeform_vertex); - break; - } -} - -void BKE_shrinkwrap_mesh_nearest_surface_deform(struct bContext *C, - Object *ob_source, - Object *ob_target) -{ - Depsgraph *depsgraph = CTX_data_depsgraph_pointer(C); - struct Scene *sce = CTX_data_scene(C); - ShrinkwrapModifierData ssmd = {{0}}; - ModifierEvalContext ctx = {depsgraph, ob_source, 0}; - int totvert; - - ssmd.target = ob_target; - ssmd.shrinkType = MOD_SHRINKWRAP_NEAREST_SURFACE; - ssmd.shrinkMode = MOD_SHRINKWRAP_ON_SURFACE; - ssmd.keepDist = 0.0f; - - Mesh *src_me = ob_source->data; - float(*vertexCos)[3] = BKE_mesh_vert_coords_alloc(src_me, &totvert); - - shrinkwrapModifier_deform(&ssmd, &ctx, sce, ob_source, src_me, NULL, -1, vertexCos, totvert); - - BKE_mesh_vert_coords_apply(src_me, vertexCos); - - MEM_freeN(vertexCos); -} - -void BKE_shrinkwrap_remesh_target_project(Mesh *src_me, Mesh *target_me, Object *ob_target) -{ - ShrinkwrapModifierData ssmd = {{0}}; - int totvert; - - ssmd.target = ob_target; - ssmd.shrinkType = MOD_SHRINKWRAP_PROJECT; - ssmd.shrinkMode = MOD_SHRINKWRAP_ON_SURFACE; - ssmd.shrinkOpts = MOD_SHRINKWRAP_PROJECT_ALLOW_NEG_DIR | MOD_SHRINKWRAP_PROJECT_ALLOW_POS_DIR; - ssmd.keepDist = 0.0f; - - /* Tolerance value to prevent artifacts on sharp edges of a mesh. - * This constant and based on experimenting with different values. */ - const float projLimitTolerance = 5.0f; - ssmd.projLimit = target_me->remesh_voxel_size * projLimitTolerance; - - float(*vertexCos)[3] = BKE_mesh_vert_coords_alloc(src_me, &totvert); - - ShrinkwrapCalcData calc = NULL_ShrinkwrapCalcData; - - calc.smd = &ssmd; - calc.numVerts = src_me->totvert; - calc.vertexCos = vertexCos; - calc.vert_normals = BKE_mesh_vertex_normals_ensure(src_me); - calc.vgroup = -1; - calc.target = target_me; - calc.keepDist = ssmd.keepDist; - calc.vert = BKE_mesh_verts_for_write(src_me); - BLI_SPACE_TRANSFORM_SETUP(&calc.local2target, ob_target, ob_target); - - ShrinkwrapTreeData tree; - if (BKE_shrinkwrap_init_tree(&tree, calc.target, ssmd.shrinkType, ssmd.shrinkMode, false)) { - calc.tree = &tree; - TIMEIT_BENCH(shrinkwrap_calc_normal_projection(&calc), deform_project); - BKE_shrinkwrap_free_tree(&tree); - } - - BKE_mesh_vert_coords_apply(src_me, vertexCos); - - MEM_freeN(vertexCos); -} diff --git a/source/blender/blenkernel/intern/shrinkwrap.cc b/source/blender/blenkernel/intern/shrinkwrap.cc new file mode 100644 index 00000000000..2b8af00d281 --- /dev/null +++ b/source/blender/blenkernel/intern/shrinkwrap.cc @@ -0,0 +1,1590 @@ +/* SPDX-License-Identifier: GPL-2.0-or-later + * Copyright Blender Foundation. All rights reserved. */ + +/** \file + * \ingroup bke + */ + +#include +#include +#include +#include +#include +#include + +#include "DNA_gpencil_modifier_types.h" +#include "DNA_mesh_types.h" +#include "DNA_meshdata_types.h" +#include "DNA_modifier_types.h" +#include "DNA_object_types.h" + +#include "BLI_math.h" +#include "BLI_math_solvers.h" +#include "BLI_task.h" +#include "BLI_utildefines.h" + +#include "BKE_DerivedMesh.h" +#include "BKE_cdderivedmesh.h" +#include "BKE_context.h" +#include "BKE_lattice.h" +#include "BKE_lib_id.h" +#include "BKE_modifier.h" +#include "BKE_shrinkwrap.h" + +#include "BKE_deform.h" +#include "BKE_editmesh.h" +#include "BKE_mesh.h" /* for OMP limits. */ +#include "BKE_mesh_runtime.h" +#include "BKE_mesh_wrapper.h" +#include "BKE_subsurf.h" + +#include "DEG_depsgraph_query.h" + +#include "MEM_guardedalloc.h" + +#include "BLI_strict_flags.h" + +/* for timing... */ +#if 0 +# include "PIL_time_utildefines.h" +#else +# define TIMEIT_BENCH(expr, id) (expr) +#endif + +/* Util macros */ +#define OUT_OF_MEMORY() ((void)printf("Shrinkwrap: Out of memory\n")) + +struct ShrinkwrapCalcData { + ShrinkwrapModifierData *smd; /* shrinkwrap modifier data */ + + Object *ob; /* object we are applying shrinkwrap to */ + + MVert *vert; /* Array of verts being projected. */ + const float (*vert_normals)[3]; + float (*vertexCos)[3]; /* vertexs being shrinkwraped */ + int numVerts; + + const MDeformVert *dvert; /* Pointer to mdeform array */ + int vgroup; /* Vertex group num */ + bool invert_vgroup; /* invert vertex group influence */ + + Mesh *target; /* mesh we are shrinking to */ + SpaceTransform local2target; /* transform to move between local and target space */ + ShrinkwrapTreeData *tree; /* mesh BVH tree data */ + + Object *aux_target; + + float keepDist; /* Distance to keep above target surface (units are in local space) */ +}; + +struct ShrinkwrapCalcCBData { + ShrinkwrapCalcData *calc; + + ShrinkwrapTreeData *tree; + ShrinkwrapTreeData *aux_tree; + + float *proj_axis; + SpaceTransform *local2aux; +}; + +bool BKE_shrinkwrap_needs_normals(int shrinkType, int shrinkMode) +{ + return (shrinkType == MOD_SHRINKWRAP_TARGET_PROJECT) || + (shrinkType != MOD_SHRINKWRAP_NEAREST_VERTEX && + shrinkMode == MOD_SHRINKWRAP_ABOVE_SURFACE); +} + +bool BKE_shrinkwrap_init_tree( + ShrinkwrapTreeData *data, Mesh *mesh, int shrinkType, int shrinkMode, bool force_normals) +{ + memset(data, 0, sizeof(*data)); + + if (mesh == nullptr) { + return false; + } + + /* We could create a BVH tree from the edit mesh, + * however accessing normals from the face/loop normals gets more involved. + * Convert mesh data since this isn't typically used in edit-mode. */ + BKE_mesh_wrapper_ensure_mdata(mesh); + + if (mesh->totvert <= 0) { + return false; + } + + data->mesh = mesh; + data->polys = BKE_mesh_polys(mesh); + + if (shrinkType == MOD_SHRINKWRAP_NEAREST_VERTEX) { + data->bvh = BKE_bvhtree_from_mesh_get(&data->treeData, mesh, BVHTREE_FROM_VERTS, 2); + + return data->bvh != nullptr; + } + + if (mesh->totpoly <= 0) { + return false; + } + + data->bvh = BKE_bvhtree_from_mesh_get(&data->treeData, mesh, BVHTREE_FROM_LOOPTRI, 4); + + if (data->bvh == nullptr) { + return false; + } + + if (force_normals || BKE_shrinkwrap_needs_normals(shrinkType, shrinkMode)) { + data->pnors = BKE_mesh_poly_normals_ensure(mesh); + if ((mesh->flag & ME_AUTOSMOOTH) != 0) { + data->clnors = static_cast(CustomData_get_layer(&mesh->ldata, CD_NORMAL)); + } + } + + if (shrinkType == MOD_SHRINKWRAP_TARGET_PROJECT) { + data->boundary = mesh->runtime.shrinkwrap_data; + } + + return true; +} + +void BKE_shrinkwrap_free_tree(ShrinkwrapTreeData *data) +{ + free_bvhtree_from_mesh(&data->treeData); +} + +void BKE_shrinkwrap_discard_boundary_data(Mesh *mesh) +{ + ShrinkwrapBoundaryData *data = mesh->runtime.shrinkwrap_data; + + if (data != nullptr) { + 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 = nullptr; +} + +/* 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(Mesh *mesh) +{ + const MVert *mvert = BKE_mesh_verts(mesh); + const MEdge *medge = BKE_mesh_edges(mesh); + const MLoop *mloop = BKE_mesh_loops(mesh); + + /* Count faces per edge (up to 2). */ + char *edge_mode = static_cast( + MEM_calloc_arrayN((size_t)mesh->totedge, sizeof(char), __func__)); + + for (int i = 0; i < mesh->totloop; i++) { + uint 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"); + uint 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 nullptr. */ + if (num_boundary_edges == 0) { + MEM_freeN(edge_is_boundary); + MEM_freeN(edge_mode); + return nullptr; + } + + /* Allocate the data object. */ + ShrinkwrapBoundaryData *data = MEM_cnew(__func__); + + 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 = static_cast( + MEM_calloc_arrayN((size_t)mesh->totvert, sizeof(int), __func__)); + + 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; + } + } + + uint 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 = static_cast( + MEM_calloc_arrayN(num_boundary_verts, sizeof(*boundary_verts), __func__)); + + signed char *vert_status = static_cast( + 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. */ + const float(*vert_normals)[3] = BKE_mesh_vertex_normals_ensure(mesh); + for (int i = 0; i < mesh->totvert; i++) { + int bidx = vert_boundary_id[i]; + + if (bidx >= 0) { + ShrinkwrapBoundaryVertData *vdata = &boundary_verts[bidx]; + float tmp[3]; + + normalize_v3(vdata->direction); + + cross_v3_v3v3(tmp, vert_normals[i], vdata->direction); + cross_v3_v3v3(vdata->normal_plane, tmp, vert_normals[i]); + normalize_v3(vdata->normal_plane); + } + } + + data->boundary_verts = boundary_verts; + + MEM_freeN(edge_mode); + return data; +} + +void BKE_shrinkwrap_compute_boundary_data(Mesh *mesh) +{ + BKE_shrinkwrap_discard_boundary_data(mesh); + + mesh->runtime.shrinkwrap_data = shrinkwrap_build_boundary_data(mesh); +} + +/** + * Shrink-wrap to the nearest vertex + * + * it builds a #BVHTree of vertices we can attach to and then + * for each vertex performs a nearest vertex search on the tree + */ +static void shrinkwrap_calc_nearest_vertex_cb_ex(void *__restrict userdata, + const int i, + const TaskParallelTLS *__restrict tls) +{ + ShrinkwrapCalcCBData *data = static_cast(userdata); + + ShrinkwrapCalcData *calc = data->calc; + BVHTreeFromMesh *treeData = &data->tree->treeData; + BVHTreeNearest *nearest = static_cast(tls->userdata_chunk); + + float *co = calc->vertexCos[i]; + float tmp_co[3]; + float weight = BKE_defvert_array_find_weight_safe(calc->dvert, i, calc->vgroup); + + if (calc->invert_vgroup) { + weight = 1.0f - weight; + } + + if (weight == 0.0f) { + return; + } + + /* Convert the vertex to tree coordinates */ + if (calc->vert) { + copy_v3_v3(tmp_co, calc->vert[i].co); + } + else { + copy_v3_v3(tmp_co, co); + } + BLI_space_transform_apply(&calc->local2target, tmp_co); + + /* Use local proximity heuristics (to reduce the nearest search) + * + * 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); + } + else { + nearest->dist_sq = FLT_MAX; + } + + BLI_bvhtree_find_nearest(treeData->tree, tmp_co, nearest, treeData->nearest_callback, treeData); + + /* Found the nearest vertex */ + if (nearest->index != -1) { + /* Adjusting the vertex weight, + * so that after interpolating it keeps a certain distance from the nearest position */ + if (nearest->dist_sq > FLT_EPSILON) { + const float dist = sqrtf(nearest->dist_sq); + weight *= (dist - calc->keepDist) / dist; + } + + /* Convert the coordinates back to mesh coordinates */ + copy_v3_v3(tmp_co, nearest->co); + BLI_space_transform_invert(&calc->local2target, tmp_co); + + interp_v3_v3v3(co, co, tmp_co, weight); /* linear interpolation */ + } +} + +static void shrinkwrap_calc_nearest_vertex(ShrinkwrapCalcData *calc) +{ + BVHTreeNearest nearest = NULL_BVHTreeNearest; + + /* Setup nearest */ + nearest.index = -1; + nearest.dist_sq = FLT_MAX; + + ShrinkwrapCalcCBData data{}; + data.calc = calc; + data.tree = calc->tree; + TaskParallelSettings settings; + BLI_parallel_range_settings_defaults(&settings); + settings.use_threading = (calc->numVerts > BKE_MESH_OMP_LIMIT); + settings.userdata_chunk = &nearest; + settings.userdata_chunk_size = sizeof(nearest); + BLI_task_parallel_range( + 0, calc->numVerts, &data, shrinkwrap_calc_nearest_vertex_cb_ex, &settings); +} + +bool BKE_shrinkwrap_project_normal(char options, + const float vert[3], + const float dir[3], + const float ray_radius, + const SpaceTransform *transf, + ShrinkwrapTreeData *tree, + BVHTreeRayHit *hit) +{ + /* don't use this because this dist value could be incompatible + * this value used by the callback for comparing previous/new dist values. + * also, at the moment there is no need to have a corrected 'dist' value */ + // #define USE_DIST_CORRECT + + float tmp_co[3], tmp_no[3]; + const float *co, *no; + BVHTreeRayHit hit_tmp; + + /* Copy from hit (we need to convert hit rays from one space coordinates to the other */ + memcpy(&hit_tmp, hit, sizeof(hit_tmp)); + + /* Apply space transform (TODO readjust dist) */ + if (transf) { + copy_v3_v3(tmp_co, vert); + BLI_space_transform_apply(transf, tmp_co); + co = tmp_co; + + copy_v3_v3(tmp_no, dir); + BLI_space_transform_apply_normal(transf, tmp_no); + no = tmp_no; + +#ifdef USE_DIST_CORRECT + hit_tmp.dist *= mat4_to_scale(((SpaceTransform *)transf)->local2target); +#endif + } + else { + co = vert; + no = dir; + } + + hit_tmp.index = -1; + + BLI_bvhtree_ray_cast( + tree->bvh, co, no, ray_radius, &hit_tmp, tree->treeData.raycast_callback, &tree->treeData); + + if (hit_tmp.index != -1) { + /* invert the normal first so face culling works on rotated objects */ + if (transf) { + BLI_space_transform_invert_normal(transf, hit_tmp.no); + } + + if (options & MOD_SHRINKWRAP_CULL_TARGET_MASK) { + /* Apply back-face. */ + const float dot = dot_v3v3(dir, hit_tmp.no); + if (((options & MOD_SHRINKWRAP_CULL_TARGET_FRONTFACE) && dot <= 0.0f) || + ((options & MOD_SHRINKWRAP_CULL_TARGET_BACKFACE) && dot >= 0.0f)) { + return false; /* Ignore hit */ + } + } + + if (transf) { + /* Inverting space transform (TODO: make coherent with the initial dist readjust). */ + BLI_space_transform_invert(transf, hit_tmp.co); +#ifdef USE_DIST_CORRECT + hit_tmp.dist = len_v3v3(vert, hit_tmp.co); +#endif + } + + BLI_assert(hit_tmp.dist <= hit->dist); + + memcpy(hit, &hit_tmp, sizeof(hit_tmp)); + return true; + } + return false; +} + +static void shrinkwrap_calc_normal_projection_cb_ex(void *__restrict userdata, + const int i, + const TaskParallelTLS *__restrict tls) +{ + ShrinkwrapCalcCBData *data = static_cast(userdata); + + ShrinkwrapCalcData *calc = data->calc; + ShrinkwrapTreeData *tree = data->tree; + ShrinkwrapTreeData *aux_tree = data->aux_tree; + + float *proj_axis = data->proj_axis; + SpaceTransform *local2aux = data->local2aux; + + BVHTreeRayHit *hit = static_cast(tls->userdata_chunk); + + const float proj_limit_squared = calc->smd->projLimit * calc->smd->projLimit; + float *co = calc->vertexCos[i]; + float tmp_co[3], tmp_no[3]; + float weight = BKE_defvert_array_find_weight_safe(calc->dvert, i, calc->vgroup); + + if (calc->invert_vgroup) { + weight = 1.0f - weight; + } + + if (weight == 0.0f) { + return; + } + + if (calc->vert != nullptr && calc->smd->projAxis == MOD_SHRINKWRAP_PROJECT_OVER_NORMAL) { + /* calc->vert contains verts from evaluated mesh. */ + /* These coordinates are deformed by vertexCos only for normal projection + * (to get correct normals) for other cases calc->verts contains undeformed coordinates and + * vertexCos should be used */ + copy_v3_v3(tmp_co, calc->vert[i].co); + copy_v3_v3(tmp_no, calc->vert_normals[i]); + } + else { + copy_v3_v3(tmp_co, co); + copy_v3_v3(tmp_no, proj_axis); + } + + hit->index = -1; + + /* TODO: we should use FLT_MAX here, but sweepsphere code isn't prepared for that */ + hit->dist = BVH_RAYCAST_DIST_MAX; + + bool is_aux = false; + + /* Project over positive direction of axis */ + if (calc->smd->shrinkOpts & MOD_SHRINKWRAP_PROJECT_ALLOW_POS_DIR) { + if (aux_tree) { + if (BKE_shrinkwrap_project_normal(0, tmp_co, tmp_no, 0.0, local2aux, aux_tree, hit)) { + is_aux = true; + } + } + + if (BKE_shrinkwrap_project_normal( + calc->smd->shrinkOpts, tmp_co, tmp_no, 0.0, &calc->local2target, tree, hit)) { + is_aux = false; + } + } + + /* Project over negative direction of axis */ + if (calc->smd->shrinkOpts & MOD_SHRINKWRAP_PROJECT_ALLOW_NEG_DIR) { + float inv_no[3]; + negate_v3_v3(inv_no, tmp_no); + + char options = calc->smd->shrinkOpts; + + if ((options & MOD_SHRINKWRAP_INVERT_CULL_TARGET) && + (options & MOD_SHRINKWRAP_CULL_TARGET_MASK)) { + options ^= MOD_SHRINKWRAP_CULL_TARGET_MASK; + } + + if (aux_tree) { + if (BKE_shrinkwrap_project_normal(0, tmp_co, inv_no, 0.0, local2aux, aux_tree, hit)) { + is_aux = true; + } + } + + if (BKE_shrinkwrap_project_normal( + options, tmp_co, inv_no, 0.0, &calc->local2target, tree, hit)) { + is_aux = false; + } + } + + /* don't set the initial dist (which is more efficient), + * because its calculated in the targets space, we want the dist in our own space */ + if (proj_limit_squared != 0.0f) { + if (hit->index != -1 && len_squared_v3v3(hit->co, co) > proj_limit_squared) { + hit->index = -1; + } + } + + if (hit->index != -1) { + if (is_aux) { + BKE_shrinkwrap_snap_point_to_surface(aux_tree, + local2aux, + calc->smd->shrinkMode, + hit->index, + hit->co, + hit->no, + calc->keepDist, + tmp_co, + hit->co); + } + else { + BKE_shrinkwrap_snap_point_to_surface(tree, + &calc->local2target, + calc->smd->shrinkMode, + hit->index, + hit->co, + hit->no, + calc->keepDist, + tmp_co, + hit->co); + } + + interp_v3_v3v3(co, co, hit->co, weight); + } +} + +static void shrinkwrap_calc_normal_projection(ShrinkwrapCalcData *calc) +{ + /* Options about projection direction */ + float proj_axis[3] = {0.0f, 0.0f, 0.0f}; + + /* Ray-cast and tree stuff. */ + + /** \note 'hit.dist' is kept in the targets space, this is only used + * for finding the best hit, to get the real dist, + * measure the len_v3v3() from the input coord to hit.co */ + BVHTreeRayHit hit; + + /* auxiliary target */ + Mesh *auxMesh = nullptr; + ShrinkwrapTreeData *aux_tree = nullptr; + ShrinkwrapTreeData aux_tree_stack; + SpaceTransform local2aux; + + /* If the user doesn't allows to project in any direction of projection axis + * then there's nothing todo. */ + if ((calc->smd->shrinkOpts & + (MOD_SHRINKWRAP_PROJECT_ALLOW_POS_DIR | MOD_SHRINKWRAP_PROJECT_ALLOW_NEG_DIR)) == 0) { + return; + } + + /* Prepare data to retrieve the direction in which we should project each vertex */ + if (calc->smd->projAxis == MOD_SHRINKWRAP_PROJECT_OVER_NORMAL) { + if (calc->vert == nullptr) { + return; + } + } + else { + /* The code supports any axis that is a combination of X,Y,Z + * although currently UI only allows to set the 3 different axis */ + if (calc->smd->projAxis & MOD_SHRINKWRAP_PROJECT_OVER_X_AXIS) { + proj_axis[0] = 1.0f; + } + if (calc->smd->projAxis & MOD_SHRINKWRAP_PROJECT_OVER_Y_AXIS) { + proj_axis[1] = 1.0f; + } + if (calc->smd->projAxis & MOD_SHRINKWRAP_PROJECT_OVER_Z_AXIS) { + proj_axis[2] = 1.0f; + } + + normalize_v3(proj_axis); + + /* Invalid projection direction */ + if (len_squared_v3(proj_axis) < FLT_EPSILON) { + return; + } + } + + if (calc->aux_target) { + auxMesh = BKE_modifier_get_evaluated_mesh_from_evaluated_object(calc->aux_target); + if (!auxMesh) { + return; + } + BLI_SPACE_TRANSFORM_SETUP(&local2aux, calc->ob, calc->aux_target); + } + + if (BKE_shrinkwrap_init_tree( + &aux_tree_stack, auxMesh, calc->smd->shrinkType, calc->smd->shrinkMode, false)) { + aux_tree = &aux_tree_stack; + } + + /* After successfully build the trees, start projection vertices. */ + ShrinkwrapCalcCBData data{}; + data.calc = calc; + data.tree = calc->tree; + data.aux_tree = aux_tree; + data.proj_axis = proj_axis; + data.local2aux = &local2aux; + TaskParallelSettings settings; + BLI_parallel_range_settings_defaults(&settings); + settings.use_threading = (calc->numVerts > BKE_MESH_OMP_LIMIT); + settings.userdata_chunk = &hit; + settings.userdata_chunk_size = sizeof(hit); + BLI_task_parallel_range( + 0, calc->numVerts, &data, shrinkwrap_calc_normal_projection_cb_ex, &settings); + + /* free data structures */ + if (aux_tree) { + BKE_shrinkwrap_free_tree(aux_tree); + } +} + +/* + * 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 + +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]; +}; + +/* 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 = static_cast(userdata); + + const 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 = static_cast(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 * /*userdata*/, + const float x[3], + float step[3], + float x_next[3]) +{ + /* Insignificant correction threshold */ + const float epsilon = 1e-5f; + /* Dot product threshold for checking if step is 'clearly' pointing outside. */ + const float dir_epsilon = 0.5f; + bool fixed = false, locked = false; + + /* The barycentric coordinate domain is a triangle bounded by + * the X and Y axes, plus the x+y=1 diagonal. First, clamp the + * movement against the diagonal. Note that step is subtracted. */ + float sum = x[0] + x[1]; + float sstep = -(step[0] + step[1]); + + if (sum + sstep > 1.0f) { + float ldist = 1.0f - sum; + + /* If already at the boundary, slide along it. */ + if (ldist < epsilon * float(M_SQRT2)) { + float step_len = len_v2(step); + + /* Abort if the solution is clearly outside the domain. */ + if (step_len > epsilon && sstep > step_len * dir_epsilon * float(M_SQRT2)) { + return false; + } + + /* Project the new position onto the diagonal. */ + add_v2_fl(step, (sum + sstep - 1.0f) * 0.5f); + fixed = locked = true; + } + else { + /* Scale a significant step down to arrive at the boundary. */ + mul_v3_fl(step, ldist / sstep); + fixed = true; + } + } + + /* Weight 0 and 1 boundary checks - along axis. */ + for (int i = 0; i < 2; i++) { + if (step[i] > x[i]) { + /* If already at the boundary, slide along it. */ + if (x[i] < epsilon) { + float step_len = len_v2(step); + + /* Abort if the solution is clearly outside the domain. */ + if (step_len > epsilon && (locked || step[i] > step_len * dir_epsilon)) { + return false; + } + + /* Reset precision errors to stay at the boundary. */ + step[i] = x[i]; + fixed = true; + } + else { + /* Scale a significant step down to arrive at the boundary. */ + mul_v3_fl(step, x[i] / step[i]); + 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 magnitude_estimate = dist + len_manhattan_v3(vtri_co[0]) + len_manhattan_v3(vtri_co[1]) + + len_manhattan_v3(vtri_co[2]); + float epsilon = magnitude_estimate * 1.0e-6f; + + /* 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{}; + tri_data.vtri_co = vtri_co; + tri_data.vtri_no = vtri_no; + tri_data.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 = &data->edge[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]; + copy_v3_v3(vedge_no[0], data->vert_normals[edge->v1]); + copy_v3_v3(vedge_no[1], data->vert_normals[edge->v2]); + + 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 */ + copy_v3_v3(vtri_no[0], tree->treeData.vert_normals[loop[0]->v]); + copy_v3_v3(vtri_no[1], tree->treeData.vert_normals[loop[1]->v]); + copy_v3_v3(vtri_no[2], tree->treeData.vert_normals[loop[2]->v]); + + /* 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]); + } + } + } +} + +void BKE_shrinkwrap_find_nearest_surface(ShrinkwrapTreeData *tree, + BVHTreeNearest *nearest, + float co[3], + int type) +{ + BVHTreeFromMesh *treeData = &tree->treeData; + + if (type == MOD_SHRINKWRAP_TARGET_PROJECT) { +#ifdef TRACE_TARGET_PROJECT + printf("\n====== 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\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 + * NN matches for each vertex + */ +static void shrinkwrap_calc_nearest_surface_point_cb_ex(void *__restrict userdata, + const int i, + const TaskParallelTLS *__restrict tls) +{ + ShrinkwrapCalcCBData *data = static_cast(userdata); + + ShrinkwrapCalcData *calc = data->calc; + BVHTreeNearest *nearest = static_cast(tls->userdata_chunk); + + float *co = calc->vertexCos[i]; + float tmp_co[3]; + float weight = BKE_defvert_array_find_weight_safe(calc->dvert, i, calc->vgroup); + + if (calc->invert_vgroup) { + weight = 1.0f - weight; + } + + if (weight == 0.0f) { + return; + } + + /* Convert the vertex to tree coordinates */ + if (calc->vert) { + copy_v3_v3(tmp_co, calc->vert[i].co); + } + else { + copy_v3_v3(tmp_co, co); + } + BLI_space_transform_apply(&calc->local2target, tmp_co); + + /* Use local proximity heuristics (to reduce the nearest search) + * + * 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) { + 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; + } + + BKE_shrinkwrap_find_nearest_surface(data->tree, nearest, tmp_co, calc->smd->shrinkType); + + /* Found the nearest vertex */ + if (nearest->index != -1) { + BKE_shrinkwrap_snap_point_to_surface(data->tree, + nullptr, + calc->smd->shrinkMode, + nearest->index, + nearest->co, + nearest->no, + calc->keepDist, + tmp_co, + tmp_co); + + /* Convert the coordinates back to mesh coordinates */ + BLI_space_transform_invert(&calc->local2target, tmp_co); + interp_v3_v3v3(co, co, tmp_co, weight); /* linear interpolation */ + } +} + +void BKE_shrinkwrap_compute_smooth_normal(const ShrinkwrapTreeData *tree, + const SpaceTransform *transform, + int looptri_idx, + const float hit_co[3], + const float hit_no[3], + float r_no[3]) +{ + const BVHTreeFromMesh *treeData = &tree->treeData; + const MLoopTri *tri = &treeData->looptri[looptri_idx]; + const float(*vert_normals)[3] = tree->treeData.vert_normals; + + /* Interpolate smooth normals if enabled. */ + if ((tree->polys[tri->poly].flag & ME_SMOOTH) != 0) { + const uint32_t vert_indices[3] = {treeData->loop[tri->tri[0]].v, + treeData->loop[tri->tri[1]].v, + treeData->loop[tri->tri[2]].v}; + float w[3], no[3][3], tmp_co[3]; + + /* Custom and auto smooth split normals. */ + if (tree->clnors) { + copy_v3_v3(no[0], tree->clnors[tri->tri[0]]); + copy_v3_v3(no[1], tree->clnors[tri->tri[1]]); + copy_v3_v3(no[2], tree->clnors[tri->tri[2]]); + } + /* Ordinary vertex normals. */ + else { + copy_v3_v3(no[0], vert_normals[vert_indices[0]]); + copy_v3_v3(no[1], vert_normals[vert_indices[1]]); + copy_v3_v3(no[2], vert_normals[vert_indices[2]]); + } + + /* Barycentric weights from hit point. */ + copy_v3_v3(tmp_co, hit_co); + + if (transform) { + BLI_space_transform_apply(transform, tmp_co); + } + + interp_weights_tri_v3(w, + treeData->vert[vert_indices[0]].co, + treeData->vert[vert_indices[1]].co, + treeData->vert[vert_indices[2]].co, + tmp_co); + + /* Interpolate using weights. */ + interp_v3_v3v3v3(r_no, no[0], no[1], no[2], w); + + if (transform) { + BLI_space_transform_invert_normal(transform, r_no); + } + else { + normalize_v3(r_no); + } + } + /* Use the polygon normal if flat. */ + else if (tree->pnors != nullptr) { + copy_v3_v3(r_no, tree->pnors[tri->poly]); + } + /* Finally fallback to the looptri normal. */ + else { + copy_v3_v3(r_no, hit_no); + } +} + +/* Helper for MOD_SHRINKWRAP_INSIDE, MOD_SHRINKWRAP_OUTSIDE and MOD_SHRINKWRAP_OUTSIDE_SURFACE. */ +static void shrinkwrap_snap_with_side(float r_point_co[3], + const float point_co[3], + const float hit_co[3], + const float hit_no[3], + float goal_dist, + float forcesign, + bool forcesnap) +{ + float delta[3]; + sub_v3_v3v3(delta, point_co, hit_co); + + float dist = len_v3(delta); + + /* If exactly on the surface, push out along normal */ + if (dist < FLT_EPSILON) { + if (forcesnap || goal_dist > 0) { + madd_v3_v3v3fl(r_point_co, hit_co, hit_no, goal_dist * forcesign); + } + else { + copy_v3_v3(r_point_co, hit_co); + } + } + /* Move to the correct side if needed */ + else { + float dsign = signf(dot_v3v3(delta, hit_no)); + + if (forcesign == 0.0f) { + forcesign = dsign; + } + + /* If on the wrong side or too close, move to correct */ + if (forcesnap || dsign * dist * forcesign < goal_dist) { + mul_v3_fl(delta, dsign / dist); + + /* At very small distance, blend in the hit normal to stabilize math. */ + float dist_epsilon = (fabsf(goal_dist) + len_manhattan_v3(hit_co)) * 1e-4f; + + if (dist < dist_epsilon) { +#ifdef TRACE_TARGET_PROJECT + printf("zero_factor %g = %g / %g\n", dist / dist_epsilon, dist, dist_epsilon); +#endif + + interp_v3_v3v3(delta, hit_no, delta, dist / dist_epsilon); + } + + madd_v3_v3v3fl(r_point_co, hit_co, delta, goal_dist * forcesign); + } + else { + copy_v3_v3(r_point_co, point_co); + } + } +} + +void BKE_shrinkwrap_snap_point_to_surface(const ShrinkwrapTreeData *tree, + const SpaceTransform *transform, + int mode, + int hit_idx, + const float hit_co[3], + const float hit_no[3], + float goal_dist, + const float point_co[3], + float r_point_co[3]) +{ + float tmp[3]; + + switch (mode) { + /* Offsets along the line between point_co and hit_co. */ + case MOD_SHRINKWRAP_ON_SURFACE: + if (goal_dist != 0) { + shrinkwrap_snap_with_side(r_point_co, point_co, hit_co, hit_no, goal_dist, 0, true); + } + else { + copy_v3_v3(r_point_co, hit_co); + } + break; + + case MOD_SHRINKWRAP_INSIDE: + shrinkwrap_snap_with_side(r_point_co, point_co, hit_co, hit_no, goal_dist, -1, false); + break; + + case MOD_SHRINKWRAP_OUTSIDE: + shrinkwrap_snap_with_side(r_point_co, point_co, hit_co, hit_no, goal_dist, +1, false); + break; + + case MOD_SHRINKWRAP_OUTSIDE_SURFACE: + if (goal_dist != 0) { + shrinkwrap_snap_with_side(r_point_co, point_co, hit_co, hit_no, goal_dist, +1, true); + } + else { + copy_v3_v3(r_point_co, hit_co); + } + break; + + /* Offsets along the normal */ + case MOD_SHRINKWRAP_ABOVE_SURFACE: + if (goal_dist != 0) { + BKE_shrinkwrap_compute_smooth_normal(tree, transform, hit_idx, hit_co, hit_no, tmp); + madd_v3_v3v3fl(r_point_co, hit_co, tmp, goal_dist); + } + else { + copy_v3_v3(r_point_co, hit_co); + } + break; + + default: + printf("Unknown Shrinkwrap surface snap mode: %d\n", mode); + copy_v3_v3(r_point_co, hit_co); + } +} + +static void shrinkwrap_calc_nearest_surface_point(ShrinkwrapCalcData *calc) +{ + BVHTreeNearest nearest = NULL_BVHTreeNearest; + + /* Setup nearest */ + nearest.index = -1; + nearest.dist_sq = FLT_MAX; + + /* Find the nearest vertex */ + ShrinkwrapCalcCBData data{}; + data.calc = calc; + data.tree = calc->tree; + TaskParallelSettings settings; + BLI_parallel_range_settings_defaults(&settings); + settings.use_threading = (calc->numVerts > BKE_MESH_OMP_LIMIT); + settings.userdata_chunk = &nearest; + settings.userdata_chunk_size = sizeof(nearest); + BLI_task_parallel_range( + 0, calc->numVerts, &data, shrinkwrap_calc_nearest_surface_point_cb_ex, &settings); +} + +void shrinkwrapModifier_deform(ShrinkwrapModifierData *smd, + const ModifierEvalContext *ctx, + Scene *scene, + Object *ob, + Mesh *mesh, + const MDeformVert *dvert, + const int defgrp_index, + float (*vertexCos)[3], + int numVerts) +{ + + DerivedMesh *ss_mesh = nullptr; + ShrinkwrapCalcData calc = NULL_ShrinkwrapCalcData; + + /* remove loop dependencies on derived meshes (TODO should this be done elsewhere?) */ + if (smd->target == ob) { + smd->target = nullptr; + } + if (smd->auxTarget == ob) { + smd->auxTarget = nullptr; + } + + /* Configure Shrinkwrap calc data */ + calc.smd = smd; + calc.ob = ob; + calc.numVerts = numVerts; + calc.vertexCos = vertexCos; + calc.dvert = dvert; + calc.vgroup = defgrp_index; + calc.invert_vgroup = (smd->shrinkOpts & MOD_SHRINKWRAP_INVERT_VGROUP) != 0; + + if (smd->target != nullptr) { + Object *ob_target = DEG_get_evaluated_object(ctx->depsgraph, smd->target); + calc.target = BKE_modifier_get_evaluated_mesh_from_evaluated_object(ob_target); + + /* TODO: there might be several "bugs" with non-uniform scales matrices + * because it will no longer be nearest surface, not sphere projection + * because space has been deformed */ + BLI_SPACE_TRANSFORM_SETUP(&calc.local2target, ob, ob_target); + + /* TODO: smd->keepDist is in global units.. must change to local */ + calc.keepDist = smd->keepDist; + } + calc.aux_target = DEG_get_evaluated_object(ctx->depsgraph, smd->auxTarget); + + if (mesh != nullptr && smd->shrinkType == MOD_SHRINKWRAP_PROJECT) { + /* Setup arrays to get vertexs positions, normals and deform weights */ + calc.vert = BKE_mesh_verts_for_write(mesh); + calc.vert_normals = BKE_mesh_vertex_normals_ensure(mesh); + + /* Using vertexs positions/normals as if a subsurface was applied */ + if (smd->subsurfLevels) { + SubsurfModifierData ssmd = {{nullptr}}; + ssmd.subdivType = ME_CC_SUBSURF; /* catmull clark */ + ssmd.levels = smd->subsurfLevels; /* levels */ + + /* TODO: to be moved to Mesh once we are done with changes in subsurf code. */ + DerivedMesh *dm = CDDM_from_mesh(mesh); + + ss_mesh = subsurf_make_derived_from_derived( + dm, + &ssmd, + scene, + nullptr, + (ob->mode & OB_MODE_EDIT) ? SUBSURF_IN_EDIT_MODE : SubsurfFlags(0)); + + if (ss_mesh) { + calc.vert = static_cast(ss_mesh->getVertDataArray(ss_mesh, CD_MVERT)); + if (calc.vert) { + /* TRICKY: this code assumes subsurface will have the transformed original vertices + * in their original order at the end of the vert array. */ + calc.vert = calc.vert + ss_mesh->getNumVerts(ss_mesh) - dm->getNumVerts(dm); + } + } + + /* Just to make sure we are not leaving any memory behind */ + BLI_assert(ssmd.emCache == nullptr); + BLI_assert(ssmd.mCache == nullptr); + + dm->release(dm); + } + } + + /* Projecting target defined - lets work! */ + ShrinkwrapTreeData tree; + + if (BKE_shrinkwrap_init_tree(&tree, calc.target, smd->shrinkType, smd->shrinkMode, false)) { + calc.tree = &tree; + + switch (smd->shrinkType) { + case MOD_SHRINKWRAP_NEAREST_SURFACE: + case MOD_SHRINKWRAP_TARGET_PROJECT: + TIMEIT_BENCH(shrinkwrap_calc_nearest_surface_point(&calc), deform_surface); + break; + + case MOD_SHRINKWRAP_PROJECT: + TIMEIT_BENCH(shrinkwrap_calc_normal_projection(&calc), deform_project); + break; + + case MOD_SHRINKWRAP_NEAREST_VERTEX: + TIMEIT_BENCH(shrinkwrap_calc_nearest_vertex(&calc), deform_vertex); + break; + } + + BKE_shrinkwrap_free_tree(&tree); + } + + /* free memory */ + if (ss_mesh) { + ss_mesh->release(ss_mesh); + } +} + +void shrinkwrapGpencilModifier_deform(ShrinkwrapGpencilModifierData *mmd, + Object *ob, + MDeformVert *dvert, + const int defgrp_index, + float (*vertexCos)[3], + int numVerts) +{ + + ShrinkwrapCalcData calc = NULL_ShrinkwrapCalcData; + /* Convert gpencil struct to use the same struct and function used with meshes. */ + ShrinkwrapModifierData smd; + smd.target = mmd->target; + smd.auxTarget = mmd->aux_target; + smd.keepDist = mmd->keep_dist; + smd.shrinkType = mmd->shrink_type; + smd.shrinkOpts = mmd->shrink_opts; + smd.shrinkMode = mmd->shrink_mode; + smd.projLimit = mmd->proj_limit; + smd.projAxis = mmd->proj_axis; + + /* Configure Shrinkwrap calc data. */ + calc.smd = &smd; + calc.ob = ob; + calc.numVerts = numVerts; + calc.vertexCos = vertexCos; + calc.dvert = dvert; + calc.vgroup = defgrp_index; + calc.invert_vgroup = (mmd->flag & GP_SHRINKWRAP_INVERT_VGROUP) != 0; + + BLI_SPACE_TRANSFORM_SETUP(&calc.local2target, ob, mmd->target); + calc.keepDist = mmd->keep_dist; + calc.tree = mmd->cache_data; + + switch (mmd->shrink_type) { + case MOD_SHRINKWRAP_NEAREST_SURFACE: + case MOD_SHRINKWRAP_TARGET_PROJECT: + TIMEIT_BENCH(shrinkwrap_calc_nearest_surface_point(&calc), gpdeform_surface); + break; + + case MOD_SHRINKWRAP_PROJECT: + TIMEIT_BENCH(shrinkwrap_calc_normal_projection(&calc), gpdeform_project); + break; + + case MOD_SHRINKWRAP_NEAREST_VERTEX: + TIMEIT_BENCH(shrinkwrap_calc_nearest_vertex(&calc), gpdeform_vertex); + break; + } +} + +void BKE_shrinkwrap_mesh_nearest_surface_deform(bContext *C, Object *ob_source, Object *ob_target) +{ + Depsgraph *depsgraph = CTX_data_depsgraph_pointer(C); + Scene *sce = CTX_data_scene(C); + ShrinkwrapModifierData ssmd = {{0}}; + ModifierEvalContext ctx = {depsgraph, ob_source, ModifierApplyFlag(0)}; + int totvert; + + ssmd.target = ob_target; + ssmd.shrinkType = MOD_SHRINKWRAP_NEAREST_SURFACE; + ssmd.shrinkMode = MOD_SHRINKWRAP_ON_SURFACE; + ssmd.keepDist = 0.0f; + + Mesh *src_me = static_cast(ob_source->data); + float(*vertexCos)[3] = BKE_mesh_vert_coords_alloc(src_me, &totvert); + + shrinkwrapModifier_deform(&ssmd, &ctx, sce, ob_source, src_me, nullptr, -1, vertexCos, totvert); + + BKE_mesh_vert_coords_apply(src_me, vertexCos); + + MEM_freeN(vertexCos); +} + +void BKE_shrinkwrap_remesh_target_project(Mesh *src_me, Mesh *target_me, Object *ob_target) +{ + ShrinkwrapModifierData ssmd = {{0}}; + int totvert; + + ssmd.target = ob_target; + ssmd.shrinkType = MOD_SHRINKWRAP_PROJECT; + ssmd.shrinkMode = MOD_SHRINKWRAP_ON_SURFACE; + ssmd.shrinkOpts = MOD_SHRINKWRAP_PROJECT_ALLOW_NEG_DIR | MOD_SHRINKWRAP_PROJECT_ALLOW_POS_DIR; + ssmd.keepDist = 0.0f; + + /* Tolerance value to prevent artifacts on sharp edges of a mesh. + * This constant and based on experimenting with different values. */ + const float projLimitTolerance = 5.0f; + ssmd.projLimit = target_me->remesh_voxel_size * projLimitTolerance; + + float(*vertexCos)[3] = BKE_mesh_vert_coords_alloc(src_me, &totvert); + + ShrinkwrapCalcData calc = NULL_ShrinkwrapCalcData; + + calc.smd = &ssmd; + calc.numVerts = src_me->totvert; + calc.vertexCos = vertexCos; + calc.vert_normals = BKE_mesh_vertex_normals_ensure(src_me); + calc.vgroup = -1; + calc.target = target_me; + calc.keepDist = ssmd.keepDist; + calc.vert = BKE_mesh_verts_for_write(src_me); + BLI_SPACE_TRANSFORM_SETUP(&calc.local2target, ob_target, ob_target); + + ShrinkwrapTreeData tree; + if (BKE_shrinkwrap_init_tree(&tree, calc.target, ssmd.shrinkType, ssmd.shrinkMode, false)) { + calc.tree = &tree; + TIMEIT_BENCH(shrinkwrap_calc_normal_projection(&calc), deform_project); + BKE_shrinkwrap_free_tree(&tree); + } + + BKE_mesh_vert_coords_apply(src_me, vertexCos); + + MEM_freeN(vertexCos); +} -- cgit v1.2.3