archive/blender-archive
Archived
1
1
Fork 0
Code Pull Requests Packages Activity
This repository has been archived on 2023-10-09. You can view files and clone it. You cannot open issues or pull requests or push a commit.
Files
517afcc858b09ecb51e73b2cfbcaffeba75a9933
blender-archive/source/blender/blenkernel/intern/shrinkwrap.c
Hans Goudey cfa53e0fbe Refactor: Move normals out of MVert, lazy calculation 
As described in T91186, this commit moves mesh vertex normals into a
contiguous array of float vectors in a custom data layer, how face
normals are currently stored.

The main interface is documented in `BKE_mesh.h`. Vertex and face
normals are now calculated on-demand and cached, retrieved with an
"ensure" function. Since the logical state of a mesh is now "has
normals when necessary", they can be retrieved from a `const` mesh.

The goal is to use on-demand calculation for all derived data, but
leave room for eager calculation for performance purposes (modifier
evaluation is threaded, but viewport data generation is not).

**Benefits**
This moves us closer to a SoA approach rather than the current AoS
paradigm. Accessing a contiguous `float3` is much more efficient than
retrieving data from a larger struct. The memory requirements for
accessing only normals or vertex locations are smaller, and at the
cost of more memory usage for just normals, they now don't have to
be converted between float and short, which also simplifies code

In the future, the remaining items can be removed from `MVert`,
leaving only `float3`, which has similar benefits (see T93602).

Removing the combination of derived and original data makes it
conceptually simpler to only calculate normals when necessary.
This is especially important now that we have more opportunities
for temporary meshes in geometry nodes.

**Performance**
In addition to the theoretical future performance improvements by
making `MVert == float3`, I've done some basic performance testing
on this patch directly. The data is fairly rough, but it gives an idea
about where things stand generally.
 - Mesh line primitive 4m Verts: 1.16x faster (36 -> 31 ms),
   showing that accessing just `MVert` is now more efficient.
 - Spring Splash Screen: 1.03-1.06 -> 1.06-1.11 FPS, a very slight
   change that at least shows there is no regression.
 - Sprite Fright Snail Smoosh: 3.30-3.40 -> 3.42-3.50 FPS, a small
   but observable speedup.
 - Set Position Node with Scaled Normal: 1.36x faster (53 -> 39 ms),
   shows that using normals in geometry nodes is faster.
 - Normal Calculation 1.6m Vert Cube: 1.19x faster (25 -> 21 ms),
   shows that calculating normals is slightly faster now.
 - File Size of 1.6m Vert Cube: 1.03x smaller (214.7 -> 208.4 MB),
   Normals are not saved in files, which can help with large meshes.

As for memory usage, it may be slightly more in some cases, but
I didn't observe any difference in the production files I tested.

**Tests**
Some modifiers and cycles test results need to be updated with this
commit, for two reasons:
 - The subdivision surface modifier is not responsible for calculating
   normals anymore. In master, the modifier creates different normals
   than the result of the `Mesh` normal calculation, so this is a bug
   fix.
 - There are small differences in the results of some modifiers that
   use normals because they are not converted to and from `short`
   anymore.

**Future improvements**
 - Remove `ModifierTypeInfo::dependsOnNormals`. Code in each modifier
   already retrieves normals if they are needed anyway.
 - Copy normals as part of a better CoW system for attributes.
 - Make more areas use lazy instead of eager normal calculation.
 - Remove `BKE_mesh_normals_tag_dirty` in more places since that is
   now the default state of a new mesh.
 - Possibly apply a similar change to derived face corner normals.

Differential Revision: https://developer.blender.org/D12770
2022-01-13 14:38:25 -06:00

1607 lines
50 KiB
C
Raw Blame History

/*
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License
* as published by the Free Software Foundation; either version 2
* of the License, or (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software Foundation,
* Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*
* The Original Code is Copyright (C) Blender Foundation.
* All rights reserved.
*/
/** \file
* \ingroup bke
*/
#include <float.h>
#include <math.h>
#include <memory.h>
#include <stdio.h>
#include <string.h>
#include <time.h>
#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;
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;
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 MLoop *mloop = mesh->mloop;
const MEdge *medge = mesh->medge;
const MVert *mvert = mesh->mvert;
/* Count faces per edge (up to 2). */
char *edge_mode = MEM_calloc_arrayN((size_t)mesh->totedge, sizeof(char), __func__);
for (int i = 0; i < mesh->totloop; i++) {
unsigned int eidx = mloop[i].e;
if (edge_mode[eidx] < 2) {
edge_mode[eidx]++;
}
}
/* Build the boundary edge bitmask. */
BLI_bitmap *edge_is_boundary = BLI_BITMAP_NEW(mesh->totedge,
"ShrinkwrapBoundaryData::edge_is_boundary");
unsigned int num_boundary_edges = 0;
for (int i = 0; i < mesh->totedge; i++) {
edge_mode[i] = (edge_mode[i] == 1);
if (edge_mode[i]) {
BLI_BITMAP_ENABLE(edge_is_boundary, i);
num_boundary_edges++;
}
}
/* If no boundary, return NULL. */
if (num_boundary_edges == 0) {
MEM_freeN(edge_is_boundary);
MEM_freeN(edge_mode);
return NULL;
}
/* Allocate the data object. */
ShrinkwrapBoundaryData *data = MEM_callocN(sizeof(ShrinkwrapBoundaryData),
"ShrinkwrapBoundaryData");
data->edge_is_boundary = edge_is_boundary;
/* Build the boundary looptri bitmask. */
const MLoopTri *mlooptri = BKE_mesh_runtime_looptri_ensure(mesh);
int totlooptri = BKE_mesh_runtime_looptri_len(mesh);
BLI_bitmap *looptri_has_boundary = BLI_BITMAP_NEW(totlooptri,
"ShrinkwrapBoundaryData::looptri_is_boundary");
for (int i = 0; i < totlooptri; i++) {
int edges[3];
BKE_mesh_looptri_get_real_edges(mesh, &mlooptri[i], edges);
for (int j = 0; j < 3; j++) {
if (edges[j] >= 0 && edge_mode[edges[j]]) {
BLI_BITMAP_ENABLE(looptri_has_boundary, i);
break;
}
}
}
data->looptri_has_boundary = looptri_has_boundary;
/* Find boundary vertices and build a mapping table for compact storage of data. */
int *vert_boundary_id = MEM_calloc_arrayN(
(size_t)mesh->totvert, sizeof(int), "ShrinkwrapBoundaryData::vert_boundary_id");
for (int i = 0; i < mesh->totedge; i++) {
if (edge_mode[i]) {
const MEdge *edge = &medge[i];
vert_boundary_id[edge->v1] = 1;
vert_boundary_id[edge->v2] = 1;
}
}
unsigned int num_boundary_verts = 0;
for (int i = 0; i < mesh->totvert; i++) {
vert_boundary_id[i] = (vert_boundary_id[i] != 0) ? (int)num_boundary_verts++ : -1;
}
data->vert_boundary_id = vert_boundary_id;
data->num_boundary_verts = num_boundary_verts;
/* Compute average directions. */
ShrinkwrapBoundaryVertData *boundary_verts = MEM_calloc_arrayN(
num_boundary_verts, sizeof(*boundary_verts), "ShrinkwrapBoundaryData::boundary_verts");
signed char *vert_status = MEM_calloc_arrayN(num_boundary_verts, sizeof(char), __func__);
for (int i = 0; i < mesh->totedge; i++) {
if (edge_mode[i]) {
const MEdge *edge = &medge[i];
float dir[3];
sub_v3_v3v3(dir, mvert[edge->v2].co, mvert[edge->v1].co);
normalize_v3(dir);
merge_vert_dir(boundary_verts, vert_status, vert_boundary_id[edge->v1], dir, 1);
merge_vert_dir(boundary_verts, vert_status, vert_boundary_id[edge->v2], dir, 2);
}
}
MEM_freeN(vert_status);
/* Finalize average direction and compute normal. */
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, false);
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 = &tree->mesh->medge[eidx];
const float *vedge_co[2] = {data->vert[edge->v1].co, data->vert[edge->v2].co};
#ifdef TRACE_TARGET_PROJECT
printf("EDGE %d (%.3f,%.3f,%.3f) (%.3f,%.3f,%.3f)\n",
eidx,
UNPACK3(vedge_co[0]),
UNPACK3(vedge_co[1]));
#endif
/* Retrieve boundary vertex IDs */
const int *vert_boundary_id = tree->boundary->vert_boundary_id;
int bid1 = vert_boundary_id[edge->v1], bid2 = vert_boundary_id[edge->v2];
if (bid1 < 0 || bid2 < 0) {
return;
}
/* Retrieve boundary vertex normals and align them to direction. */
const ShrinkwrapBoundaryVertData *boundary_verts = tree->boundary->boundary_verts;
float vedge_dir[2][3], dir[3];
copy_v3_v3(vedge_dir[0], boundary_verts[bid1].normal_plane);
copy_v3_v3(vedge_dir[1], boundary_verts[bid2].normal_plane);
sub_v3_v3v3(dir, vedge_co[1], vedge_co[0]);
if (dot_v3v3(boundary_verts[bid1].direction, dir) < 0) {
negate_v3(vedge_dir[0]);
}
if (dot_v3v3(boundary_verts[bid2].direction, dir) < 0) {
negate_v3(vedge_dir[1]);
}
/* Solve a quadratic equation: lerp(d0,d1,x) * (co - lerp(v0,v1,x)) = 0 */
float d0v0 = dot_v3v3(vedge_dir[0], vedge_co[0]), d0v1 = dot_v3v3(vedge_dir[0], vedge_co[1]);
float d1v0 = dot_v3v3(vedge_dir[1], vedge_co[0]), d1v1 = dot_v3v3(vedge_dir[1], vedge_co[1]);
float d0co = dot_v3v3(vedge_dir[0], co);
float a = d0v1 - d0v0 + d1v0 - d1v1;
float b = 2 * d0v0 - d0v1 - d0co - d1v0 + dot_v3v3(vedge_dir[1], co);
float c = d0co - d0v0;
float det = b * b - 4 * a * c;
if (det >= 0) {
const float epsilon = 1e-6f;
float sdet = sqrtf(det);
float hit_co[3], hit_no[3];
for (int i = (det > 0 ? 2 : 0); i >= 0; i -= 2) {
float x = (-b + ((float)i - 1) * sdet) / (2 * a);
if (x >= -epsilon && x <= 1.0f + epsilon) {
CLAMP(x, 0, 1);
float vedge_no[2][3];
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->mesh->mpoly[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,
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, false);
/* 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 = mesh->mvert;
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 = src_me->mvert;
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);
}
View Git Blame Copy Permalink