Shrinkwrap: new mode that projects along the target normal.

The Nearest Surface Point shrink method, while fast, is neither
smooth nor continuous: as the source point moves, the projected
point can both stop and jump. This causes distortions in the
deformation of the shrinkwrap modifier, and the motion of an
animated object with a shrinkwrap constraint.

This patch implements a new mode, which, instead of using the simple
nearest point search, iteratively solves an equation for each triangle
to find a point which has its interpolated normal point to or from the
original vertex. Non-manifold boundary edges are treated as infinitely
thin cylinders that cast normals in all perpendicular directions.

Since this is useful for the constraint, and having multiple
objects with constraints targeting the same guide mesh is a quite
reasonable use case, rather than calculating the mesh boundary edge
data over and over again, it is precomputed and cached in the mesh.

Reviewers: mont29

Differential Revision: https://developer.blender.org/D3836
This commit is contained in:
Alexander Gavrilov 2018-11-06 21:04:53 +03:00
parent 0709fac41e
commit f600b4bc67
17 changed files with 719 additions and 13 deletions

View File

@ -748,7 +748,7 @@ class ConstraintButtonsPanel:
layout.prop(con, "distance")
layout.prop(con, "shrinkwrap_type")
if con.shrinkwrap_type in {'PROJECT', 'NEAREST_SURFACE'}:
if con.shrinkwrap_type in {'PROJECT', 'NEAREST_SURFACE', 'TARGET_PROJECT'}:
layout.prop(con, 'wrap_mode', text="Snap Mode")
if con.shrinkwrap_type == 'PROJECT':
@ -769,7 +769,7 @@ class ConstraintButtonsPanel:
rowsub.prop(con, "use_invert_cull")
layout.prop(con, "project_limit")
if con.shrinkwrap_type in {'PROJECT', 'NEAREST_SURFACE'}:
if con.shrinkwrap_type in {'PROJECT', 'NEAREST_SURFACE', 'TARGET_PROJECT'}:
layout.prop(con, "use_track_normal")
row = layout.row(align=True)

View File

@ -828,7 +828,7 @@ class DATA_PT_modifiers(ModifierButtonsPanel, Panel):
col.label(text="Mode:")
col.prop(md, "wrap_method", text="")
if md.wrap_method in {'PROJECT', 'NEAREST_SURFACEPOINT'}:
if md.wrap_method in {'PROJECT', 'NEAREST_SURFACEPOINT', 'TARGET_PROJECT'}:
col.prop(md, "wrap_mode", text="")
if md.wrap_method == 'PROJECT':

View File

@ -33,6 +33,7 @@
/* Shrinkwrap stuff */
#include "BKE_bvhutils.h"
#include "BLI_bitmap.h"
/*
* Shrinkwrap is composed by a set of functions and options that define the type of shrink.
@ -55,6 +56,34 @@ struct ShrinkwrapModifierData;
struct BVHTree;
struct SpaceTransform;
/* Information about boundary edges in the mesh. */
typedef struct ShrinkwrapBoundaryVertData {
/* Average direction of edges that meet here. */
float direction[3];
/* Closest vector to direction that is orthogonal to vertex normal. */
float normal_plane[3];
} ShrinkwrapBoundaryVertData;
typedef struct ShrinkwrapBoundaryData {
/* True if the edge belongs to exactly one face. */
const BLI_bitmap *edge_is_boundary;
/* True if the looptri has any boundary edges. */
const BLI_bitmap *looptri_has_boundary;
/* Mapping from vertex index to boundary vertex index, or -1.
* Used for compact storage of data about boundary vertices. */
const int *vert_boundary_id;
unsigned int num_boundary_verts;
/* Direction data about boundary vertices. */
const ShrinkwrapBoundaryVertData *boundary_verts;
} ShrinkwrapBoundaryData;
void BKE_shrinkwrap_discard_boundary_data(struct Mesh *mesh);
void BKE_shrinkwrap_compute_boundary_data(struct Mesh *mesh);
/* Information about a mesh and BVH tree. */
typedef struct ShrinkwrapTreeData {
Mesh *mesh;
@ -62,6 +91,7 @@ typedef struct ShrinkwrapTreeData {
BVHTreeFromMesh treeData;
float (*clnors)[3];
ShrinkwrapBoundaryData *boundary;
} ShrinkwrapTreeData;
/* Checks if the modifier needs target normals with these settings. */
@ -91,6 +121,10 @@ bool BKE_shrinkwrap_project_normal(
char options, const float vert[3], const float dir[3], const float ray_radius,
const struct SpaceTransform *transf, struct ShrinkwrapTreeData *tree, BVHTreeRayHit *hit);
/* Maps the point to the nearest surface, either by simple nearest, or by target normal projection. */
void BKE_shrinkwrap_find_nearest_surface(
struct ShrinkwrapTreeData *tree, struct BVHTreeNearest *nearest, float co[3], int type);
/* Computes a smooth normal of the target (if applicable) at the hit location. */
void BKE_shrinkwrap_compute_smooth_normal(
const struct ShrinkwrapTreeData *tree, const struct SpaceTransform *transform,

View File

@ -75,6 +75,7 @@
#include "DEG_depsgraph.h"
#include "DEG_depsgraph_query.h"
#include "BKE_shrinkwrap.h"
#ifdef WITH_OPENSUBDIV
# include "DNA_userdef_types.h"
@ -1980,6 +1981,15 @@ static void mesh_finalize_eval(Object *object)
}
}
static void mesh_build_extra_data(struct Depsgraph *depsgraph, Object *ob)
{
uint32_t eval_flags = DEG_get_eval_flags_for_id(depsgraph, &ob->id);
if (eval_flags & DAG_EVAL_NEED_SHRINKWRAP_BOUNDARY) {
BKE_shrinkwrap_compute_boundary_data(ob->runtime.mesh_eval);
}
}
static void mesh_build_data(
struct Depsgraph *depsgraph, Scene *scene, Object *ob, CustomDataMask dataMask,
const bool build_shapekey_layers, const bool need_mapping)
@ -2016,6 +2026,8 @@ static void mesh_build_data(
}
BLI_assert(!(ob->derivedFinal->dirty & DM_DIRTY_NORMALS));
mesh_build_extra_data(depsgraph, ob);
}
static void editbmesh_build_data(

View File

@ -3636,6 +3636,7 @@ static void shrinkwrap_get_tarmat(struct Depsgraph *UNUSED(depsgraph), bConstrai
switch (scon->shrinkType) {
case MOD_SHRINKWRAP_NEAREST_SURFACE:
case MOD_SHRINKWRAP_NEAREST_VERTEX:
case MOD_SHRINKWRAP_TARGET_PROJECT:
{
BVHTreeNearest nearest;
@ -3644,15 +3645,14 @@ static void shrinkwrap_get_tarmat(struct Depsgraph *UNUSED(depsgraph), bConstrai
BLI_space_transform_apply(&transform, co);
BVHTreeFromMesh *treeData = &tree.treeData;
BLI_bvhtree_find_nearest(treeData->tree, co, &nearest, treeData->nearest_callback, treeData);
BKE_shrinkwrap_find_nearest_surface(&tree, &nearest, co, scon->shrinkType);
if (nearest.index < 0) {
fail = true;
break;
}
if (scon->shrinkType == MOD_SHRINKWRAP_NEAREST_SURFACE) {
if (scon->shrinkType != MOD_SHRINKWRAP_NEAREST_VERTEX) {
if (do_track_normal) {
track_normal = true;
BKE_shrinkwrap_compute_smooth_normal(&tree, NULL, nearest.index, nearest.co, nearest.no, track_no);

View File

@ -563,6 +563,7 @@ void BKE_mesh_copy_data(Main *bmain, Mesh *me_dst, const Mesh *me_src, const int
me_dst->runtime.batch_cache = NULL;
me_dst->runtime.looptris.array = NULL;
me_dst->runtime.bvh_cache = NULL;
me_dst->runtime.shrinkwrap_data = NULL;
if (me_src->id.tag & LIB_TAG_NO_MAIN) {
me_dst->runtime.deformed_only = me_src->runtime.deformed_only;

View File

@ -44,6 +44,7 @@
#include "BKE_mesh.h"
#include "BKE_mesh_runtime.h"
#include "BKE_subdiv_ccg.h"
#include "BKE_shrinkwrap.h"
/* -------------------------------------------------------------------- */
/** \name Mesh Runtime Struct Utils
@ -202,6 +203,7 @@ void BKE_mesh_runtime_clear_geometry(Mesh *mesh)
BKE_subdiv_ccg_destroy(mesh->runtime.subdiv_ccg);
mesh->runtime.subdiv_ccg = NULL;
}
BKE_shrinkwrap_discard_boundary_data(mesh);
}
/** \} */

View File

@ -45,6 +45,7 @@
#include "BLI_math.h"
#include "BLI_utildefines.h"
#include "BLI_task.h"
#include "BLI_math_solvers.h"
#include "BKE_shrinkwrap.h"
#include "BKE_cdderivedmesh.h"
@ -57,6 +58,9 @@
#include "BKE_editmesh.h"
#include "BKE_mesh.h" /* for OMP limits. */
#include "BKE_subsurf.h"
#include "BKE_mesh_runtime.h"
#include "MEM_guardedalloc.h"
#include "BLI_strict_flags.h"
@ -70,7 +74,6 @@
/* Util macros */
#define OUT_OF_MEMORY() ((void)printf("Shrinkwrap: Out of memory\n"))
typedef struct ShrinkwrapCalcData {
ShrinkwrapModifierData *smd; //shrinkwrap modifier data
@ -106,7 +109,8 @@ typedef struct ShrinkwrapCalcCBData {
/* Checks if the modifier needs target normals with these settings. */
bool BKE_shrinkwrap_needs_normals(int shrinkType, int shrinkMode)
{
return shrinkType != MOD_SHRINKWRAP_NEAREST_VERTEX && shrinkMode == MOD_SHRINKWRAP_ABOVE_SURFACE;
return (shrinkType == MOD_SHRINKWRAP_TARGET_PROJECT) ||
(shrinkType != MOD_SHRINKWRAP_NEAREST_VERTEX && shrinkMode == MOD_SHRINKWRAP_ABOVE_SURFACE);
}
/* Initializes the mesh data structure from the given mesh and settings. */
@ -142,6 +146,10 @@ bool BKE_shrinkwrap_init_tree(ShrinkwrapTreeData *data, Mesh *mesh, int shrinkTy
}
}
if (shrinkType == MOD_SHRINKWRAP_TARGET_PROJECT) {
data->boundary = mesh->runtime.shrinkwrap_data;
}
return true;
}
}
@ -152,6 +160,176 @@ void BKE_shrinkwrap_free_tree(ShrinkwrapTreeData *data)
free_bvhtree_from_mesh(&data->treeData);
}
/* Free boundary data for target project */
void BKE_shrinkwrap_discard_boundary_data(struct Mesh *mesh)
{
struct ShrinkwrapBoundaryData *data = mesh->runtime.shrinkwrap_data;
if (data != NULL) {
MEM_freeN((void*)data->edge_is_boundary);
MEM_freeN((void*)data->looptri_has_boundary);
MEM_freeN((void*)data->vert_boundary_id);
MEM_freeN((void*)data->boundary_verts);
MEM_freeN(data);
}
mesh->runtime.shrinkwrap_data = NULL;
}
/* Accumulate edge for average boundary edge direction. */
static void merge_vert_dir(ShrinkwrapBoundaryVertData *vdata, signed char *status, int index, const float edge_dir[3], signed char side)
{
BLI_assert(index >= 0);
float *direction = vdata[index].direction;
/* Invert the direction vector if either:
* - This is the second edge and both edges have the vertex as start or end.
* - For third and above, if it points in the wrong direction.
*/
if (status[index] >= 0 ? status[index] == side : dot_v3v3(direction, edge_dir) < 0) {
sub_v3_v3(direction, edge_dir);
}
else {
add_v3_v3(direction, edge_dir);
}
status[index] = (status[index] == 0) ? side : -1;
}
static ShrinkwrapBoundaryData *shrinkwrap_build_boundary_data(struct Mesh *mesh)
{
const MLoop *mloop = mesh->mloop;
const MEdge *medge = mesh->medge;
const MVert *mvert = mesh->mvert;
/* Count faces per edge (up to 2). */
char *edge_mode = MEM_calloc_arrayN((size_t)mesh->totedge, sizeof(char), __func__);
for (int i = 0; i < mesh->totloop; i++) {
unsigned int eidx = mloop[i].e;
if (edge_mode[eidx] < 2) {
edge_mode[eidx]++;
}
}
/* Build the boundary edge bitmask. */
BLI_bitmap *edge_is_boundary = BLI_BITMAP_NEW(mesh->totedge, "ShrinkwrapBoundaryData::edge_is_boundary");
unsigned int num_boundary_edges = 0;
for (int i = 0; i < mesh->totedge; i++) {
edge_mode[i] = (edge_mode[i] == 1);
if (edge_mode[i]) {
BLI_BITMAP_ENABLE(edge_is_boundary, i);
num_boundary_edges++;
}
}
/* If no boundary, return NULL. */
if (num_boundary_edges == 0) {
MEM_freeN(edge_is_boundary);
MEM_freeN(edge_mode);
return NULL;
}
/* Allocate the data object. */
ShrinkwrapBoundaryData *data = MEM_callocN(sizeof(ShrinkwrapBoundaryData), "ShrinkwrapBoundaryData");
data->edge_is_boundary = edge_is_boundary;
/* Build the boundary looptri bitmask. */
const MLoopTri *mlooptri = BKE_mesh_runtime_looptri_ensure(mesh);
int totlooptri = BKE_mesh_runtime_looptri_len(mesh);
BLI_bitmap *looptri_has_boundary = BLI_BITMAP_NEW(totlooptri, "ShrinkwrapBoundaryData::looptri_is_boundary");
for (int i = 0; i < totlooptri; i++) {
int edges[3];
BKE_mesh_looptri_get_real_edges(mesh, &mlooptri[i], edges);
for (int j = 0; j < 3; j++) {
if (edges[j] >= 0 && edge_mode[edges[j]]) {
BLI_BITMAP_ENABLE(looptri_has_boundary, i);
break;
}
}
}
data->looptri_has_boundary = looptri_has_boundary;
/* Find boundary vertices and build a mapping table for compact storage of data. */
int *vert_boundary_id = MEM_calloc_arrayN((size_t)mesh->totvert, sizeof(int), "ShrinkwrapBoundaryData::vert_boundary_id");
for (int i = 0; i < mesh->totedge; i++) {
if (edge_mode[i]) {
const MEdge *edge = &medge[i];
vert_boundary_id[edge->v1] = 1;
vert_boundary_id[edge->v2] = 1;
}
}
unsigned int num_boundary_verts = 0;
for (int i = 0; i < mesh->totvert; i++) {
vert_boundary_id[i] = (vert_boundary_id[i] != 0) ? (int)num_boundary_verts++ : -1;
}
data->vert_boundary_id = vert_boundary_id;
data->num_boundary_verts = num_boundary_verts;
/* Compute average directions. */
ShrinkwrapBoundaryVertData *boundary_verts = MEM_calloc_arrayN(num_boundary_verts, sizeof(*boundary_verts), "ShrinkwrapBoundaryData::boundary_verts");
signed char *vert_status = MEM_calloc_arrayN(num_boundary_verts, sizeof(char), __func__);
for (int i = 0; i < mesh->totedge; i++) {
if (edge_mode[i]) {
const MEdge *edge = &medge[i];
float dir[3];
sub_v3_v3v3(dir, mvert[edge->v2].co, mvert[edge->v1].co);
normalize_v3(dir);
merge_vert_dir(boundary_verts, vert_status, vert_boundary_id[edge->v1], dir, 1);
merge_vert_dir(boundary_verts, vert_status, vert_boundary_id[edge->v2], dir, 2);
}
}
MEM_freeN(vert_status);
/* Finalize average direction and compute normal. */
for (int i = 0; i < mesh->totvert; i++) {
int bidx = vert_boundary_id[i];
if (bidx >= 0) {
ShrinkwrapBoundaryVertData *vdata = &boundary_verts[bidx];
float no[3], tmp[3];
normalize_v3(vdata->direction);
normal_short_to_float_v3(no, mesh->mvert[i].no);
cross_v3_v3v3(tmp, no, vdata->direction);
cross_v3_v3v3(vdata->normal_plane, tmp, no);
normalize_v3(vdata->normal_plane);
}
}
data->boundary_verts = boundary_verts;
MEM_freeN(edge_mode);
return data;
}
void BKE_shrinkwrap_compute_boundary_data(struct Mesh *mesh)
{
BKE_shrinkwrap_discard_boundary_data(mesh);
mesh->runtime.shrinkwrap_data = shrinkwrap_build_boundary_data(mesh);
}
/*
* Shrinkwrap to the nearest vertex
*
@ -522,6 +700,350 @@ static void shrinkwrap_calc_normal_projection(ShrinkwrapCalcData *calc)
}
}
/*
* Shrinkwrap Target Surface Project mode
*
* It uses Newton's method to find a surface location with its
* smooth normal pointing at the original point.
*
* The equation system on barycentric weights and normal multiplier:
*
* (w0*V0 + w1*V1 + w2*V2) + l * (w0*N0 + w1*N1 + w2*N2) - CO = 0
* w0 + w1 + w2 = 1
*
* The actual solution vector is [ w0, w1, l ], with w2 eliminated.
*/
//#define TRACE_TARGET_PROJECT
typedef struct TargetProjectTriData {
const float **vtri_co;
const float (*vtri_no)[3];
const float *point_co;
float n0_minus_n2[3], n1_minus_n2[3];
float c0_minus_c2[3], c1_minus_c2[3];
/* Current interpolated position and normal. */
float co_interp[3], no_interp[3];
} TargetProjectTriData;
/* Computes the deviation of the equation system from goal. */
static void target_project_tri_deviation(void *userdata, const float x[3], float r_delta[3])
{
TargetProjectTriData *data = userdata;
float w[3] = { x[0], x[1], 1.0f - x[0] - x[1] };
interp_v3_v3v3v3(data->co_interp, data->vtri_co[0], data->vtri_co[1], data->vtri_co[2], w);
interp_v3_v3v3v3(data->no_interp, data->vtri_no[0], data->vtri_no[1], data->vtri_no[2], w);
madd_v3_v3v3fl(r_delta, data->co_interp, data->no_interp, x[2]);
sub_v3_v3(r_delta, data->point_co);
}
/* Computes the Jacobian matrix of the equation system. */
static void target_project_tri_jacobian(void *userdata, const float x[3], float r_jacobian[3][3])
{
TargetProjectTriData *data = userdata;
madd_v3_v3v3fl(r_jacobian[0], data->c0_minus_c2, data->n0_minus_n2, x[2]);
madd_v3_v3v3fl(r_jacobian[1], data->c1_minus_c2, data->n1_minus_n2, x[2]);
copy_v3_v3(r_jacobian[2], data->vtri_no[2]);
madd_v3_v3fl(r_jacobian[2], data->n0_minus_n2, x[0]);
madd_v3_v3fl(r_jacobian[2], data->n1_minus_n2, x[1]);
}
/* Clamp barycentric weights to the triangle. */
static void target_project_tri_clamp(float x[3])
{
if (x[0] < 0.0f) {
x[0] = 0.0f;
}
if (x[1] < 0.0f) {
x[1] = 0.0f;
}
if (x[0] + x[1] > 1.0f) {
x[0] = x[0] / (x[0] + x[1]);
x[1] = 1.0f - x[0];
}
}
/* Correct the Newton's method step to keep the coordinates within the triangle. */
static bool target_project_tri_correct(void *UNUSED(userdata), const float x[3], float step[3], float x_next[3])
{
/* Insignificant correction threshold */
const float epsilon = 1e-6f;
const float dir_epsilon = 0.05f;
bool fixed = false, locked = false;
/* Weight 0 and 1 boundary check. */
for (int i = 0; i < 2; i++) {
if (step[i] > x[i]) {
if (step[i] > dir_epsilon * fabsf(step[1 - i])) {
/* Abort if the solution is clearly outside the domain. */
if (x[i] < epsilon) {
return false;
}
/* Scale a significant step down to arrive at the boundary. */
mul_v3_fl(step, x[i] / step[i]);
fixed = true;
}
else {
/* Reset precision errors to stay at the boundary. */
step[i] = x[i];
fixed = locked = true;
}
}
}
/* Weight 2 boundary check. */
float sum = x[0] + x[1];
float sstep = step[0] + step[1];
if (sum - sstep > 1.0f) {
if (sstep < -dir_epsilon * (fabsf(step[0]) + fabsf(step[1]))) {
/* Abort if the solution is clearly outside the domain. */
if (sum > 1.0f - epsilon) {
return false;
}
/* Scale a significant step down to arrive at the boundary. */
mul_v3_fl(step, (1.0f - sum) / -sstep);
fixed = true;
}
else {
/* Reset precision errors to stay at the boundary. */
if (locked) {
step[0] = step[1] = 0.0f;
}
else {
step[0] -= 0.5f * sstep;
step[1] = -step[0];
fixed = true;
}
}
}
/* Recompute and clamp the new coordinates after step correction. */
if (fixed) {
sub_v3_v3v3(x_next, x, step);
target_project_tri_clamp(x_next);
}
return true;
}
static bool target_project_solve_point_tri(
const float *vtri_co[3], const float vtri_no[3][3], const float point_co[3],
const float hit_co[3], float hit_dist_sq,
float r_hit_co[3], float r_hit_no[3])
{
float x[3], tmp[3];
float dist = sqrtf(hit_dist_sq);
float epsilon = dist * 1.0e-5f;
CLAMP_MIN(epsilon, 1.0e-5f);
/* Initial solution vector: barycentric weights plus distance along normal. */
interp_weights_tri_v3(x, UNPACK3(vtri_co), hit_co);
interp_v3_v3v3v3(r_hit_no, UNPACK3(vtri_no), x);
sub_v3_v3v3(tmp, point_co, hit_co);
x[2] = (dot_v3v3(tmp, r_hit_no) < 0) ? -dist : dist;
/* Solve the equations iteratively. */
TargetProjectTriData tri_data = { .vtri_co = vtri_co, .vtri_no = vtri_no, .point_co = point_co };
sub_v3_v3v3(tri_data.n0_minus_n2, vtri_no[0], vtri_no[2]);
sub_v3_v3v3(tri_data.n1_minus_n2, vtri_no[1], vtri_no[2]);
sub_v3_v3v3(tri_data.c0_minus_c2, vtri_co[0], vtri_co[2]);
sub_v3_v3v3(tri_data.c1_minus_c2, vtri_co[1], vtri_co[2]);
target_project_tri_clamp(x);
#ifdef TRACE_TARGET_PROJECT
const bool trace = true;
#else
const bool trace = false;
#endif
bool ok = BLI_newton3d_solve(target_project_tri_deviation, target_project_tri_jacobian, target_project_tri_correct, &tri_data, epsilon, 20, trace, x, x);
if (ok) {
copy_v3_v3(r_hit_co, tri_data.co_interp);
copy_v3_v3(r_hit_no, tri_data.no_interp);
return true;
}
return false;
}
static bool update_hit(BVHTreeNearest *nearest, int index, const float co[3], const float hit_co[3], const float hit_no[3])
{
float dist_sq = len_squared_v3v3(hit_co, co);
if (dist_sq < nearest->dist_sq) {
#ifdef TRACE_TARGET_PROJECT
printf("===> %d (%.3f,%.3f,%.3f) %g < %g\n", index, UNPACK3(hit_co), dist_sq, nearest->dist_sq);
#endif
nearest->index = index;
nearest->dist_sq = dist_sq;
copy_v3_v3(nearest->co, hit_co);
normalize_v3_v3(nearest->no, hit_no);
return true;
}
return false;
}
/* Target projection on a non-manifold boundary edge - treats it like an infinitely thin cylinder. */
static void target_project_edge(const ShrinkwrapTreeData *tree, int index, const float co[3], BVHTreeNearest *nearest, int eidx)
{
const BVHTreeFromMesh *data = &tree->treeData;
const MEdge *edge = &tree->mesh->medge[eidx];
const float *vedge_co[2] = { data->vert[edge->v1].co, data->vert[edge->v2].co };
#ifdef TRACE_TARGET_PROJECT
printf("EDGE %d (%.3f,%.3f,%.3f) (%.3f,%.3f,%.3f)\n", eidx, UNPACK3(vedge_co[0]), UNPACK3(vedge_co[1]));
#endif
/* Retrieve boundary vertex IDs */
const int *vert_boundary_id = tree->boundary->vert_boundary_id;
int bid1 = vert_boundary_id[edge->v1], bid2 = vert_boundary_id[edge->v2];
if (bid1 < 0 || bid2 < 0) {
return;
}
/* Retrieve boundary vertex normals and align them to direction. */
const ShrinkwrapBoundaryVertData *boundary_verts = tree->boundary->boundary_verts;
float vedge_dir[2][3], dir[3];
copy_v3_v3(vedge_dir[0], boundary_verts[bid1].normal_plane);
copy_v3_v3(vedge_dir[1], boundary_verts[bid2].normal_plane);
sub_v3_v3v3(dir, vedge_co[1], vedge_co[0]);
if (dot_v3v3(boundary_verts[bid1].direction, dir) < 0) {
negate_v3(vedge_dir[0]);
}
if (dot_v3v3(boundary_verts[bid2].direction, dir) < 0) {
negate_v3(vedge_dir[1]);
}
/* Solve a quadratic equation: lerp(d0,d1,x) * (co - lerp(v0,v1,x)) = 0 */
float d0v0 = dot_v3v3(vedge_dir[0], vedge_co[0]), d0v1 = dot_v3v3(vedge_dir[0], vedge_co[1]);
float d1v0 = dot_v3v3(vedge_dir[1], vedge_co[0]), d1v1 = dot_v3v3(vedge_dir[1], vedge_co[1]);
float d0co = dot_v3v3(vedge_dir[0], co);
float a = d0v1 - d0v0 + d1v0 - d1v1;
float b = 2 * d0v0 - d0v1 - d0co - d1v0 + dot_v3v3(vedge_dir[1], co);
float c = d0co - d0v0;
float det = b * b - 4 * a * c;
if (det >= 0) {
const float epsilon = 1e-6f;
float sdet = sqrtf(det);
float hit_co[3], hit_no[3];
for (int i = (det > 0 ? 2 : 0); i >= 0; i -= 2) {
float x = (- b + ((float)i - 1) * sdet) / (2 * a);
if (x >= -epsilon && x <= 1.0f + epsilon) {
CLAMP(x, 0, 1);
float vedge_no[2][3];
normal_short_to_float_v3(vedge_no[0], data->vert[edge->v1].no);
normal_short_to_float_v3(vedge_no[1], data->vert[edge->v2].no);
interp_v3_v3v3(hit_co, vedge_co[0], vedge_co[1], x);
interp_v3_v3v3(hit_no, vedge_no[0], vedge_no[1], x);
update_hit(nearest, index, co, hit_co, hit_no);
}
}
}
}
/* Target normal projection BVH callback - based on mesh_looptri_nearest_point. */
static void mesh_looptri_target_project(void *userdata, int index, const float co[3], BVHTreeNearest *nearest)
{
const ShrinkwrapTreeData *tree = (ShrinkwrapTreeData *)userdata;
const BVHTreeFromMesh *data = &tree->treeData;
const MLoopTri *lt = &data->looptri[index];
const MLoop *loop[3] = { &data->loop[lt->tri[0]], &data->loop[lt->tri[1]], &data->loop[lt->tri[2]] };
const MVert *vtri[3] = { &data->vert[loop[0]->v], &data->vert[loop[1]->v], &data->vert[loop[2]->v] };
const float *vtri_co[3] = { vtri[0]->co, vtri[1]->co, vtri[2]->co };
float raw_hit_co[3], hit_co[3], hit_no[3], dist_sq, vtri_no[3][3];
/* First find the closest point and bail out if it's worse than the current solution. */
closest_on_tri_to_point_v3(raw_hit_co, co, UNPACK3(vtri_co));
dist_sq = len_squared_v3v3(co, raw_hit_co);
#ifdef TRACE_TARGET_PROJECT
printf("TRIANGLE %d (%.3f,%.3f,%.3f) (%.3f,%.3f,%.3f) (%.3f,%.3f,%.3f) %g %g\n",
index, UNPACK3(vtri_co[0]), UNPACK3(vtri_co[1]), UNPACK3(vtri_co[2]), dist_sq, nearest->dist_sq);
#endif
if (dist_sq >= nearest->dist_sq)
return;
/* Decode normals */
normal_short_to_float_v3(vtri_no[0], vtri[0]->no);
normal_short_to_float_v3(vtri_no[1], vtri[1]->no);
normal_short_to_float_v3(vtri_no[2], vtri[2]->no);
/* Solve the equations for the triangle */
if (target_project_solve_point_tri(vtri_co, vtri_no, co, raw_hit_co, dist_sq, hit_co, hit_no)) {
update_hit(nearest, index, co, hit_co, hit_no);
}
/* Boundary edges */
else if (tree->boundary && BLI_BITMAP_TEST(tree->boundary->looptri_has_boundary, index)) {
const BLI_bitmap *is_boundary = tree->boundary->edge_is_boundary;
int edges[3];
BKE_mesh_looptri_get_real_edges(tree->mesh, lt, edges);
for (int i = 0; i < 3; i++) {
if (edges[i] >= 0 && BLI_BITMAP_TEST(is_boundary, edges[i])) {
target_project_edge(tree, index, co, nearest, edges[i]);
}
}
}
}
/*
* Maps the point to the nearest surface, either by simple nearest, or by target normal projection.
*/
void BKE_shrinkwrap_find_nearest_surface(struct ShrinkwrapTreeData *tree, BVHTreeNearest *nearest, float co[3], int type)
{
BVHTreeFromMesh *treeData = &tree->treeData;
if (type == MOD_SHRINKWRAP_TARGET_PROJECT) {
#ifdef TRACE_TARGET_PROJECT
printf("====== TARGET PROJECT START ======\n");
#endif
BLI_bvhtree_find_nearest_ex(tree->bvh, co, nearest, mesh_looptri_target_project, tree, BVH_NEAREST_OPTIMAL_ORDER);
#ifdef TRACE_TARGET_PROJECT
printf("====== TARGET PROJECT END: %d %g ======\n", nearest->index, nearest->dist_sq);
#endif
if (nearest->index < 0) {
/* fallback to simple nearest */
BLI_bvhtree_find_nearest(tree->bvh, co, nearest, treeData->nearest_callback, treeData);
}
}
else {
BLI_bvhtree_find_nearest(tree->bvh, co, nearest, treeData->nearest_callback, treeData);
}
}
/*
* Shrinkwrap moving vertexs to the nearest surface point on the target
*
@ -536,7 +1058,6 @@ static void shrinkwrap_calc_nearest_surface_point_cb_ex(
ShrinkwrapCalcCBData *data = userdata;
ShrinkwrapCalcData *calc = data->calc;
BVHTreeFromMesh *treeData = &data->tree->treeData;
BVHTreeNearest *nearest = tls->userdata_chunk;
float *co = calc->vertexCos[i];
@ -565,12 +1086,20 @@ static void shrinkwrap_calc_nearest_surface_point_cb_ex(
* If we already had an hit before.. we assume this vertex is going to have a close hit to that other vertex
* so we can initiate the "nearest.dist" with the expected value to that last hit.
* This will lead in pruning of the search tree. */
if (nearest->index != -1)
nearest->dist_sq = len_squared_v3v3(tmp_co, nearest->co);
if (nearest->index != -1) {
if (calc->smd->shrinkType == MOD_SHRINKWRAP_TARGET_PROJECT) {
/* Heuristic doesn't work because of additional restrictions. */
nearest->index = -1;
nearest->dist_sq = FLT_MAX;
}
else {
nearest->dist_sq = len_squared_v3v3(tmp_co, nearest->co);
}
}
else
nearest->dist_sq = FLT_MAX;
BLI_bvhtree_find_nearest(treeData->tree, tmp_co, nearest, treeData->nearest_callback, treeData);
BKE_shrinkwrap_find_nearest_surface(data->tree, nearest, tmp_co, calc->smd->shrinkType);
/* Found the nearest vertex */
if (nearest->index != -1) {
@ -838,6 +1367,7 @@ void shrinkwrapModifier_deform(ShrinkwrapModifierData *smd, struct Scene *scene,
switch (smd->shrinkType) {
case MOD_SHRINKWRAP_NEAREST_SURFACE:
case MOD_SHRINKWRAP_TARGET_PROJECT:
TIMEIT_BENCH(shrinkwrap_calc_nearest_surface_point(&calc), deform_surface);
break;

View File

@ -53,6 +53,15 @@ void BLI_svd_m3(const float m3[3][3], float r_U[3][3], float r_S[], float r_V[3]
bool BLI_tridiagonal_solve(const float *a, const float *b, const float *c, const float *d, float *r_x, const int count);
bool BLI_tridiagonal_solve_cyclic(const float *a, const float *b, const float *c, const float *d, float *r_x, const int count);
/* Generic 3 variable Newton's method solver. */
typedef void (*Newton3D_DeltaFunc)(void *userdata, const float x[3], float r_delta[3]);
typedef void (*Newton3D_JacobianFunc)(void *userdata, const float x[3], float r_jacobian[3][3]);
typedef bool (*Newton3D_CorrectionFunc)(void *userdata, const float x[3], float step[3], float x_next[3]);
bool BLI_newton3d_solve(
Newton3D_DeltaFunc func_delta, Newton3D_JacobianFunc func_jacobian, Newton3D_CorrectionFunc func_correction, void *userdata,
float epsilon, int max_iterations, bool trace, const float x_init[3], float result[3]);
#ifdef BLI_MATH_GCC_WARN_PRAGMA
# pragma GCC diagnostic pop
#endif

View File

@ -179,3 +179,98 @@ bool BLI_tridiagonal_solve_cyclic(const float *a, const float *b, const float *c
return success;
}
/**
* \brief Solve a generic f(x) = 0 equation using Newton's method.
*
* \param func_delta Callback computing the value of f(x).
* \param func_jacobian Callback computing the Jacobian matrix of the function at x.
* \param func_correction Callback for forcing the search into an arbitrary custom domain. May be NULL.
* \param userdata Data for the callbacks.
* \param epsilon Desired precision.
* \param max_iterations Limit on the iterations.
* \param max_corrections Limit on the number of times the correction callback can fire before giving up.
* \param trace Enables logging to console.
* \param x_init Initial solution vector.
* \param result Final result.
* \return true if success
*/
bool BLI_newton3d_solve(
Newton3D_DeltaFunc func_delta, Newton3D_JacobianFunc func_jacobian, Newton3D_CorrectionFunc func_correction, void *userdata,
float epsilon, int max_iterations, bool trace, const float x_init[3], float result[3])
{
float fdelta[3], fdeltav, next_fdeltav;
float jacobian[3][3], step[3], x[3], x_next[3];
epsilon *= epsilon;
copy_v3_v3(x, x_init);
func_delta(userdata, x, fdelta);
fdeltav = len_squared_v3(fdelta);
if (trace) {
printf("START (%g, %g, %g) %g\n", x[0], x[1], x[2], fdeltav);
}
for (int i = 0; i < max_iterations && fdeltav > epsilon; i++) {
/* Newton's method step. */
func_jacobian(userdata, x, jacobian);
if (!invert_m3(jacobian)) {
return false;
}
mul_v3_m3v3(step, jacobian, fdelta);
sub_v3_v3v3(x_next, x, step);
/* Custom out-of-bounds value correction. */
if (func_correction) {
if (trace) {
printf("%3d * (%g, %g, %g)\n", i, x_next[0], x_next[1], x_next[2]);
}
if (!func_correction(userdata, x, step, x_next)) {
return false;
}
}
func_delta(userdata, x_next, fdelta);
next_fdeltav = len_squared_v3(fdelta);
if (trace) {
printf("%3d ? (%g, %g, %g) %g\n", i, x_next[0], x_next[1], x_next[2], next_fdeltav);
}
/* Line search correction. */
while (next_fdeltav > fdeltav) {
float g0 = sqrtf(fdeltav), g1 = sqrtf(next_fdeltav);
float g01 = -g0 / len_v3(step);
float det = 2.0f * (g1 - g0 - g01);
float l = (det == 0.0f) ? 0.1f : -g01 / det;
CLAMP_MIN(l, 0.1f);
mul_v3_fl(step, l);
sub_v3_v3v3(x_next, x, step);
func_delta(userdata, x_next, fdelta);
next_fdeltav = len_squared_v3(fdelta);
if (trace) {
printf("%3d . (%g, %g, %g) %g\n", i, x_next[0], x_next[1], x_next[2], next_fdeltav);
}
}
copy_v3_v3(x, x_next);
fdeltav = next_fdeltav;
}
bool success = (fdeltav <= epsilon);
if (trace) {
printf("%s (%g, %g, %g) %g\n", success ? "OK " : "FAIL", x[0], x[1], x[2], fdeltav);
}
copy_v3_v3(result, x);
return success;
}

View File

@ -80,7 +80,11 @@ enum {
* to meet dependencies with such a things as curve modifier and other guys
* who're using curve deform, where_on_path and so.
*/
DAG_EVAL_NEED_CURVE_PATH = 1,
DAG_EVAL_NEED_CURVE_PATH = (1 << 0),
/* A shrinkwrap modifier or constraint targeting this mesh needs information
* about non-manifold boundary edges for the Target Normal Project mode.
*/
DAG_EVAL_NEED_SHRINKWRAP_BOUNDARY = (1 << 1),
};
#ifdef __cplusplus

View File

@ -1019,6 +1019,9 @@ void DepsgraphRelationBuilder::build_constraints(ID *id,
if (track || BKE_shrinkwrap_needs_normals(scon->shrinkType, scon->shrinkMode)) {
add_customdata_mask(target_key, CD_MASK_NORMAL | CD_MASK_CUSTOMLOOPNORMAL);
}
if (scon->shrinkType == MOD_SHRINKWRAP_TARGET_PROJECT) {
add_special_eval_flag(&ct->tar->id, DAG_EVAL_NEED_SHRINKWRAP_BOUNDARY);
}
}
/* NOTE: obdata eval now doesn't necessarily depend on the

View File

@ -100,6 +100,9 @@ typedef struct Mesh_Runtime {
/** 'BVHCache', for 'BKE_bvhutil.c' */
struct LinkNode *bvh_cache;
/** Non-manifold boundary data for Shrinkwrap Target Project. */
struct ShrinkwrapBoundaryData *shrinkwrap_data;
/** Set by modifier stack if only deformed from original. */
char deformed_only;
/**

View File

@ -888,6 +888,7 @@ enum {
MOD_SHRINKWRAP_NEAREST_SURFACE = 0,
MOD_SHRINKWRAP_PROJECT = 1,
MOD_SHRINKWRAP_NEAREST_VERTEX = 2,
MOD_SHRINKWRAP_TARGET_PROJECT = 3,
};
/* Shrinkwrap->shrinkMode */

View File

@ -2134,6 +2134,9 @@ static void rna_def_constraint_shrinkwrap(BlenderRNA *brna)
"Shrink the location to the nearest target surface along a given axis"},
{MOD_SHRINKWRAP_NEAREST_VERTEX, "NEAREST_VERTEX", 0, "Nearest Vertex",
"Shrink the location to the nearest target vertex"},
{MOD_SHRINKWRAP_TARGET_PROJECT, "TARGET_PROJECT", 0, "Target Normal Project",
"Shrink the location to the nearest target surface "
"along the interpolated vertex normals of the target"},
{0, NULL, 0, NULL, NULL}
};

View File

@ -3167,6 +3167,9 @@ static void rna_def_modifier_shrinkwrap(BlenderRNA *brna)
"Shrink the mesh to the nearest target surface along a given axis"},
{MOD_SHRINKWRAP_NEAREST_VERTEX, "NEAREST_VERTEX", 0, "Nearest Vertex",
"Shrink the mesh to the nearest target vertex"},
{MOD_SHRINKWRAP_TARGET_PROJECT, "TARGET_PROJECT", 0, "Target Normal Project",
"Shrink the mesh to the nearest target surface "
"along the interpolated vertex normals of the target"},
{0, NULL, 0, NULL, NULL}
};

View File

@ -150,10 +150,16 @@ static void updateDepsgraph(ModifierData *md, const ModifierUpdateDepsgraphConte
if (smd->target != NULL) {
DEG_add_object_relation(ctx->node, smd->target, DEG_OB_COMP_TRANSFORM, "Shrinkwrap Modifier");
DEG_add_object_relation_with_customdata(ctx->node, smd->target, DEG_OB_COMP_GEOMETRY, mask, "Shrinkwrap Modifier");
if (smd->shrinkType == MOD_SHRINKWRAP_TARGET_PROJECT) {
DEG_add_special_eval_flag(ctx->node, &smd->target->id, DAG_EVAL_NEED_SHRINKWRAP_BOUNDARY);
}
}
if (smd->auxTarget != NULL) {
DEG_add_object_relation(ctx->node, smd->auxTarget, DEG_OB_COMP_TRANSFORM, "Shrinkwrap Modifier");
DEG_add_object_relation_with_customdata(ctx->node, smd->auxTarget, DEG_OB_COMP_GEOMETRY, mask, "Shrinkwrap Modifier");
if (smd->shrinkType == MOD_SHRINKWRAP_TARGET_PROJECT) {
DEG_add_special_eval_flag(ctx->node, &smd->auxTarget->id, DAG_EVAL_NEED_SHRINKWRAP_BOUNDARY);
}
}
DEG_add_object_relation(ctx->node, ctx->object, DEG_OB_COMP_TRANSFORM, "Shrinkwrap Modifier");
}