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

6 files changed, 588 insertions, 9 deletions

diff --git a/source/blender/blenkernel/BKE_shrinkwrap.h b/source/blender/blenkernel/BKE_shrinkwrap.h
index d1a14003b04..04c0f1c5d9a 100644
--- a/source/blender/blenkernel/BKE_shrinkwrap.h
+++ b/source/blender/blenkernel/BKE_shrinkwrap.h
@@ -33,6 +33,7 @@
 
 /* Shrinkwrap stuff */
 #include "BKE_bvhutils.h"
+#include "BLI_bitmap.h"
 
 /*
  * Shrinkwrap is composed by a set of functions and options that define the type of shrink.
@@ -55,6 +56,34 @@ struct ShrinkwrapModifierData;
 struct BVHTree;
 struct SpaceTransform;
 
+/* Information about boundary edges in the mesh. */
+typedef struct ShrinkwrapBoundaryVertData {
+	/* Average direction of edges that meet here. */
+	float direction[3];
+
+	/* Closest vector to direction that is orthogonal to vertex normal. */
+	float normal_plane[3];
+} ShrinkwrapBoundaryVertData;
+
+typedef struct ShrinkwrapBoundaryData {
+	/* True if the edge belongs to exactly one face. */
+	const BLI_bitmap *edge_is_boundary;
+	/* True if the looptri has any boundary edges. */
+	const BLI_bitmap *looptri_has_boundary;
+
+	/* Mapping from vertex index to boundary vertex index, or -1.
+	 * Used for compact storage of data about boundary vertices. */
+	const int *vert_boundary_id;
+	unsigned int num_boundary_verts;
+
+	/* Direction data about boundary vertices. */
+	const ShrinkwrapBoundaryVertData *boundary_verts;
+} ShrinkwrapBoundaryData;
+
+void BKE_shrinkwrap_discard_boundary_data(struct Mesh *mesh);
+void BKE_shrinkwrap_compute_boundary_data(struct Mesh *mesh);
+
+/* Information about a mesh and BVH tree. */
 typedef struct ShrinkwrapTreeData {
 	Mesh *mesh;
 
@@ -62,6 +91,7 @@ typedef struct ShrinkwrapTreeData {
 	BVHTreeFromMesh treeData;
 
 	float (*clnors)[3];
+	ShrinkwrapBoundaryData *boundary;
 } ShrinkwrapTreeData;
 
 /* Checks if the modifier needs target normals with these settings. */
@@ -91,6 +121,10 @@ bool BKE_shrinkwrap_project_normal(
         char options, const float vert[3], const float dir[3], const float ray_radius,
         const struct SpaceTransform *transf, struct ShrinkwrapTreeData *tree, BVHTreeRayHit *hit);
 
+/* Maps the point to the nearest surface, either by simple nearest, or by target normal projection. */
+void BKE_shrinkwrap_find_nearest_surface(
+        struct ShrinkwrapTreeData *tree, struct BVHTreeNearest *nearest, float co[3], int type);
+
 /* Computes a smooth normal of the target (if applicable) at the hit location. */
 void BKE_shrinkwrap_compute_smooth_normal(
         const struct ShrinkwrapTreeData *tree, const struct SpaceTransform *transform,
diff --git a/source/blender/blenkernel/intern/DerivedMesh.c b/source/blender/blenkernel/intern/DerivedMesh.c
index c84bd378b04..1d787388774 100644
--- a/source/blender/blenkernel/intern/DerivedMesh.c
+++ b/source/blender/blenkernel/intern/DerivedMesh.c
@@ -75,6 +75,7 @@
 
 #include "DEG_depsgraph.h"
 #include "DEG_depsgraph_query.h"
+#include "BKE_shrinkwrap.h"
 
 #ifdef WITH_OPENSUBDIV
 #  include "DNA_userdef_types.h"
@@ -1980,6 +1981,15 @@ static void mesh_finalize_eval(Object *object)
 	}
 }
 
+static void mesh_build_extra_data(struct Depsgraph *depsgraph, Object *ob)
+{
+	uint32_t eval_flags = DEG_get_eval_flags_for_id(depsgraph, &ob->id);
+
+	if (eval_flags & DAG_EVAL_NEED_SHRINKWRAP_BOUNDARY) {
+		BKE_shrinkwrap_compute_boundary_data(ob->runtime.mesh_eval);
+	}
+}
+
 static void mesh_build_data(
         struct Depsgraph *depsgraph, Scene *scene, Object *ob, CustomDataMask dataMask,
         const bool build_shapekey_layers, const bool need_mapping)
@@ -2016,6 +2026,8 @@ static void mesh_build_data(
 	}
 
 	BLI_assert(!(ob->derivedFinal->dirty & DM_DIRTY_NORMALS));
+
+	mesh_build_extra_data(depsgraph, ob);
 }
 
 static void editbmesh_build_data(
diff --git a/source/blender/blenkernel/intern/constraint.c b/source/blender/blenkernel/intern/constraint.c
index 41b07d73dc9..cc94e87c747 100644
--- a/source/blender/blenkernel/intern/constraint.c
+++ b/source/blender/blenkernel/intern/constraint.c
@@ -3636,6 +3636,7 @@ static void shrinkwrap_get_tarmat(struct Depsgraph *UNUSED(depsgraph), bConstrai
 			switch (scon->shrinkType) {
 				case MOD_SHRINKWRAP_NEAREST_SURFACE:
 				case MOD_SHRINKWRAP_NEAREST_VERTEX:
+				case MOD_SHRINKWRAP_TARGET_PROJECT:
 				{
 					BVHTreeNearest nearest;
 
@@ -3644,15 +3645,14 @@ static void shrinkwrap_get_tarmat(struct Depsgraph *UNUSED(depsgraph), bConstrai
 
 					BLI_space_transform_apply(&transform, co);
 
-					BVHTreeFromMesh *treeData = &tree.treeData;
-					BLI_bvhtree_find_nearest(treeData->tree, co, &nearest, treeData->nearest_callback, treeData);
+					BKE_shrinkwrap_find_nearest_surface(&tree, &nearest, co, scon->shrinkType);
 
 					if (nearest.index < 0) {
 						fail = true;
 						break;
 					}
 
-					if (scon->shrinkType == MOD_SHRINKWRAP_NEAREST_SURFACE) {
+					if (scon->shrinkType != MOD_SHRINKWRAP_NEAREST_VERTEX) {
 						if (do_track_normal) {
 							track_normal = true;
 							BKE_shrinkwrap_compute_smooth_normal(&tree, NULL, nearest.index, nearest.co, nearest.no, track_no);
diff --git a/source/blender/blenkernel/intern/mesh.c b/source/blender/blenkernel/intern/mesh.c
index ed01889d200..6a00aaf576b 100644
--- a/source/blender/blenkernel/intern/mesh.c
+++ b/source/blender/blenkernel/intern/mesh.c
@@ -563,6 +563,7 @@ void BKE_mesh_copy_data(Main *bmain, Mesh *me_dst, const Mesh *me_src, const int
 	me_dst->runtime.batch_cache = NULL;
 	me_dst->runtime.looptris.array = NULL;
 	me_dst->runtime.bvh_cache = NULL;
+	me_dst->runtime.shrinkwrap_data = NULL;
 
 	if (me_src->id.tag & LIB_TAG_NO_MAIN) {
 		me_dst->runtime.deformed_only = me_src->runtime.deformed_only;
diff --git a/source/blender/blenkernel/intern/mesh_runtime.c b/source/blender/blenkernel/intern/mesh_runtime.c
index 68880cacfb9..0931318f17a 100644
--- a/source/blender/blenkernel/intern/mesh_runtime.c
+++ b/source/blender/blenkernel/intern/mesh_runtime.c
@@ -44,6 +44,7 @@
 #include "BKE_mesh.h"
 #include "BKE_mesh_runtime.h"
 #include "BKE_subdiv_ccg.h"
+#include "BKE_shrinkwrap.h"
 
 /* -------------------------------------------------------------------- */
 /** \name Mesh Runtime Struct Utils
@@ -202,6 +203,7 @@ void BKE_mesh_runtime_clear_geometry(Mesh *mesh)
 		BKE_subdiv_ccg_destroy(mesh->runtime.subdiv_ccg);
 		mesh->runtime.subdiv_ccg = NULL;
 	}
+	BKE_shrinkwrap_discard_boundary_data(mesh);
 }
 
 /** \} */
diff --git a/source/blender/blenkernel/intern/shrinkwrap.c b/source/blender/blenkernel/intern/shrinkwrap.c
index 6942c8f1da3..78045e5fd56 100644
--- a/source/blender/blenkernel/intern/shrinkwrap.c
+++ b/source/blender/blenkernel/intern/shrinkwrap.c
@@ -45,6 +45,7 @@
 #include "BLI_math.h"
 #include "BLI_utildefines.h"
 #include "BLI_task.h"
+#include "BLI_math_solvers.h"
 
 #include "BKE_shrinkwrap.h"
 #include "BKE_cdderivedmesh.h"
@@ -57,6 +58,9 @@
 #include "BKE_editmesh.h"
 #include "BKE_mesh.h"  /* for OMP limits. */
 #include "BKE_subsurf.h"
+#include "BKE_mesh_runtime.h"
+
+#include "MEM_guardedalloc.h"
 
 #include "BLI_strict_flags.h"
 
@@ -70,7 +74,6 @@
 /* Util macros */
 #define OUT_OF_MEMORY() ((void)printf("Shrinkwrap: Out of memory\n"))
 
-
 typedef struct ShrinkwrapCalcData {
 	ShrinkwrapModifierData *smd;    //shrinkwrap modifier data
 
@@ -106,7 +109,8 @@ typedef struct ShrinkwrapCalcCBData {
 /* Checks if the modifier needs target normals with these settings. */
 bool BKE_shrinkwrap_needs_normals(int shrinkType, int shrinkMode)
 {
-	return shrinkType != MOD_SHRINKWRAP_NEAREST_VERTEX && shrinkMode == MOD_SHRINKWRAP_ABOVE_SURFACE;
+	return (shrinkType == MOD_SHRINKWRAP_TARGET_PROJECT) ||
+	       (shrinkType != MOD_SHRINKWRAP_NEAREST_VERTEX && shrinkMode == MOD_SHRINKWRAP_ABOVE_SURFACE);
 }
 
 /* Initializes the mesh data structure from the given mesh and settings. */
@@ -142,6 +146,10 @@ bool BKE_shrinkwrap_init_tree(ShrinkwrapTreeData *data, Mesh *mesh, int shrinkTy
 			}
 		}
 
+		if (shrinkType == MOD_SHRINKWRAP_TARGET_PROJECT) {
+			data->boundary = mesh->runtime.shrinkwrap_data;
+		}
+
 		return true;
 	}
 }
@@ -152,6 +160,176 @@ void BKE_shrinkwrap_free_tree(ShrinkwrapTreeData *data)
 	free_bvhtree_from_mesh(&data->treeData);
 }
 
+/* Free boundary data for target project */
+void BKE_shrinkwrap_discard_boundary_data(struct Mesh *mesh)
+{
+	struct ShrinkwrapBoundaryData *data = mesh->runtime.shrinkwrap_data;
+
+	if (data != NULL) {
+		MEM_freeN((void*)data->edge_is_boundary);
+		MEM_freeN((void*)data->looptri_has_boundary);
+		MEM_freeN((void*)data->vert_boundary_id);
+		MEM_freeN((void*)data->boundary_verts);
+
+		MEM_freeN(data);
+	}
+
+	mesh->runtime.shrinkwrap_data = NULL;
+}
+
+/* Accumulate edge for average boundary edge direction. */
+static void merge_vert_dir(ShrinkwrapBoundaryVertData *vdata, signed char *status, int index, const float edge_dir[3], signed char side)
+{
+	BLI_assert(index >= 0);
+	float *direction = vdata[index].direction;
+
+	/* Invert the direction vector if either:
+	 * - This is the second edge and both edges have the vertex as start or end.
+	 * - For third and above, if it points in the wrong direction.
+	 */
+	if (status[index] >= 0 ? status[index] == side : dot_v3v3(direction, edge_dir) < 0) {
+		sub_v3_v3(direction, edge_dir);
+	}
+	else {
+		add_v3_v3(direction, edge_dir);
+	}
+
+	status[index] = (status[index] == 0) ? side : -1;
+}
+
+static ShrinkwrapBoundaryData *shrinkwrap_build_boundary_data(struct Mesh *mesh)
+{
+	const MLoop *mloop = mesh->mloop;
+	const MEdge *medge = mesh->medge;
+	const MVert *mvert = mesh->mvert;
+
+	/* Count faces per edge (up to 2). */
+	char *edge_mode = MEM_calloc_arrayN((size_t)mesh->totedge, sizeof(char), __func__);
+
+	for (int i = 0; i < mesh->totloop; i++) {
+		unsigned int eidx = mloop[i].e;
+
+		if (edge_mode[eidx] < 2) {
+			edge_mode[eidx]++;
+		}
+	}
+
+	/* Build the boundary edge bitmask. */
+	BLI_bitmap *edge_is_boundary = BLI_BITMAP_NEW(mesh->totedge, "ShrinkwrapBoundaryData::edge_is_boundary");
+	unsigned int num_boundary_edges = 0;
+
+	for (int i = 0; i < mesh->totedge; i++) {
+		edge_mode[i] = (edge_mode[i] == 1);
+
+		if (edge_mode[i]) {
+			BLI_BITMAP_ENABLE(edge_is_boundary, i);
+			num_boundary_edges++;
+		}
+	}
+
+	/* If no boundary, return NULL. */
+	if (num_boundary_edges == 0) {
+		MEM_freeN(edge_is_boundary);
+		MEM_freeN(edge_mode);
+		return NULL;
+	}
+
+	/* Allocate the data object. */
+	ShrinkwrapBoundaryData *data = MEM_callocN(sizeof(ShrinkwrapBoundaryData), "ShrinkwrapBoundaryData");
+
+	data->edge_is_boundary = edge_is_boundary;
+
+	/* Build the boundary looptri bitmask. */
+	const MLoopTri *mlooptri = BKE_mesh_runtime_looptri_ensure(mesh);
+	int totlooptri = BKE_mesh_runtime_looptri_len(mesh);
+
+	BLI_bitmap *looptri_has_boundary = BLI_BITMAP_NEW(totlooptri, "ShrinkwrapBoundaryData::looptri_is_boundary");
+
+	for (int i = 0; i < totlooptri; i++) {
+		int edges[3];
+		BKE_mesh_looptri_get_real_edges(mesh, &mlooptri[i], edges);
+
+		for (int j = 0; j < 3; j++) {
+			if (edges[j] >= 0 && edge_mode[edges[j]]) {
+				BLI_BITMAP_ENABLE(looptri_has_boundary, i);
+				break;
+			}
+		}
+	}
+
+	data->looptri_has_boundary = looptri_has_boundary;
+
+	/* Find boundary vertices and build a mapping table for compact storage of data. */
+	int *vert_boundary_id = MEM_calloc_arrayN((size_t)mesh->totvert, sizeof(int), "ShrinkwrapBoundaryData::vert_boundary_id");
+
+	for (int i = 0; i < mesh->totedge; i++) {
+		if (edge_mode[i]) {
+			const MEdge *edge = &medge[i];
+
+			vert_boundary_id[edge->v1] = 1;
+			vert_boundary_id[edge->v2] = 1;
+		}
+	}
+
+	unsigned int num_boundary_verts = 0;
+
+	for (int i = 0; i < mesh->totvert; i++) {
+		vert_boundary_id[i] = (vert_boundary_id[i] != 0) ? (int)num_boundary_verts++ : -1;
+	}
+
+	data->vert_boundary_id = vert_boundary_id;
+	data->num_boundary_verts = num_boundary_verts;
+
+	/* Compute average directions. */
+	ShrinkwrapBoundaryVertData *boundary_verts = MEM_calloc_arrayN(num_boundary_verts, sizeof(*boundary_verts), "ShrinkwrapBoundaryData::boundary_verts");
+
+	signed char *vert_status = MEM_calloc_arrayN(num_boundary_verts, sizeof(char), __func__);
+
+	for (int i = 0; i < mesh->totedge; i++) {
+		if (edge_mode[i]) {
+			const MEdge *edge = &medge[i];
+
+			float dir[3];
+			sub_v3_v3v3(dir, mvert[edge->v2].co, mvert[edge->v1].co);
+			normalize_v3(dir);
+
+			merge_vert_dir(boundary_verts, vert_status, vert_boundary_id[edge->v1], dir, 1);
+			merge_vert_dir(boundary_verts, vert_status, vert_boundary_id[edge->v2], dir, 2);
+		}
+	}
+
+	MEM_freeN(vert_status);
+
+	/* Finalize average direction and compute normal. */
+	for (int i = 0; i < mesh->totvert; i++) {
+		int bidx = vert_boundary_id[i];
+
+		if (bidx >= 0) {
+			ShrinkwrapBoundaryVertData *vdata = &boundary_verts[bidx];
+			float no[3], tmp[3];
+
+			normalize_v3(vdata->direction);
+
+			normal_short_to_float_v3(no, mesh->mvert[i].no);
+			cross_v3_v3v3(tmp, no, vdata->direction);
+			cross_v3_v3v3(vdata->normal_plane, tmp, no);
+			normalize_v3(vdata->normal_plane);
+		}
+	}
+
+	data->boundary_verts = boundary_verts;
+
+	MEM_freeN(edge_mode);
+	return data;
+}
+
+void BKE_shrinkwrap_compute_boundary_data(struct Mesh *mesh)
+{
+	BKE_shrinkwrap_discard_boundary_data(mesh);
+
+	mesh->runtime.shrinkwrap_data = shrinkwrap_build_boundary_data(mesh);
+}
+
 /*
  * Shrinkwrap to the nearest vertex
  *
@@ -523,6 +701,350 @@ static void shrinkwrap_calc_normal_projection(ShrinkwrapCalcData *calc)
 }
 
 /*
+ * Shrinkwrap Target Surface Project mode
+ *
+ * It uses Newton's method to find a surface location with its
+ * smooth normal pointing at the original point.
+ *
+ * The equation system on barycentric weights and normal multiplier:
+ *
+ *   (w0*V0 + w1*V1 + w2*V2) + l * (w0*N0 + w1*N1 + w2*N2) - CO = 0
+ *   w0 + w1 + w2 = 1
+ *
+ * The actual solution vector is [ w0, w1, l ], with w2 eliminated.
+ */
+
+//#define TRACE_TARGET_PROJECT
+
+typedef struct TargetProjectTriData {
+	const float **vtri_co;
+	const float (*vtri_no)[3];
+	const float *point_co;
+
+	float n0_minus_n2[3], n1_minus_n2[3];
+	float c0_minus_c2[3], c1_minus_c2[3];
+
+	/* Current interpolated position and normal. */
+	float co_interp[3], no_interp[3];
+} TargetProjectTriData;
+
+/* Computes the deviation of the equation system from goal. */
+static void target_project_tri_deviation(void *userdata, const float x[3], float r_delta[3])
+{
+	TargetProjectTriData *data = userdata;
+
+	float w[3] = { x[0], x[1], 1.0f - x[0] - x[1] };
+	interp_v3_v3v3v3(data->co_interp, data->vtri_co[0], data->vtri_co[1], data->vtri_co[2], w);
+	interp_v3_v3v3v3(data->no_interp, data->vtri_no[0], data->vtri_no[1], data->vtri_no[2], w);
+
+	madd_v3_v3v3fl(r_delta, data->co_interp, data->no_interp, x[2]);
+	sub_v3_v3(r_delta, data->point_co);
+}
+
+/* Computes the Jacobian matrix of the equation system. */
+static void target_project_tri_jacobian(void *userdata, const float x[3], float r_jacobian[3][3])
+{
+	TargetProjectTriData *data = userdata;
+
+	madd_v3_v3v3fl(r_jacobian[0], data->c0_minus_c2, data->n0_minus_n2, x[2]);
+	madd_v3_v3v3fl(r_jacobian[1], data->c1_minus_c2, data->n1_minus_n2, x[2]);
+
+	copy_v3_v3(r_jacobian[2], data->vtri_no[2]);
+	madd_v3_v3fl(r_jacobian[2], data->n0_minus_n2, x[0]);
+	madd_v3_v3fl(r_jacobian[2], data->n1_minus_n2, x[1]);
+}
+
+/* Clamp barycentric weights to the triangle. */
+static void target_project_tri_clamp(float x[3])
+{
+	if (x[0] < 0.0f) {
+		x[0] = 0.0f;
+	}
+	if (x[1] < 0.0f) {
+		x[1] = 0.0f;
+	}
+	if (x[0] + x[1] > 1.0f) {
+		x[0] = x[0] / (x[0] + x[1]);
+		x[1] = 1.0f - x[0];
+	}
+}
+
+/* Correct the Newton's method step to keep the coordinates within the triangle. */
+static bool target_project_tri_correct(void *UNUSED(userdata), const float x[3], float step[3], float x_next[3])
+{
+	/* Insignificant correction threshold */
+	const float epsilon = 1e-6f;
+	const float dir_epsilon = 0.05f;
+	bool fixed = false, locked = false;
+
+	/* Weight 0 and 1 boundary check. */
+	for (int i = 0; i < 2; i++) {
+		if (step[i] > x[i]) {
+			if (step[i] > dir_epsilon * fabsf(step[1 - i])) {
+				/* Abort if the solution is clearly outside the domain. */
+				if (x[i] < epsilon) {
+					return false;
+				}
+
+				/* Scale a significant step down to arrive at the boundary. */
+				mul_v3_fl(step, x[i] / step[i]);
+				fixed = true;
+			}
+			else {
+				/* Reset precision errors to stay at the boundary. */
+				step[i] = x[i];
+				fixed = locked = true;
+			}
+		}
+	}
+
+	/* Weight 2 boundary check. */
+	float sum = x[0] + x[1];
+	float sstep = step[0] + step[1];
+
+	if (sum - sstep > 1.0f) {
+		if (sstep < -dir_epsilon * (fabsf(step[0]) + fabsf(step[1]))) {
+			/* Abort if the solution is clearly outside the domain. */
+			if (sum > 1.0f - epsilon) {
+				return false;
+			}
+
+			/* Scale a significant step down to arrive at the boundary. */
+			mul_v3_fl(step, (1.0f - sum) / -sstep);
+			fixed = true;
+		}
+		else {
+			/* Reset precision errors to stay at the boundary. */
+			if (locked) {
+				step[0] = step[1] = 0.0f;
+			}
+			else {
+				step[0] -= 0.5f * sstep;
+				step[1] = -step[0];
+				fixed = true;
+			}
+		}
+	}
+
+	/* Recompute and clamp the new coordinates after step correction. */
+	if (fixed) {
+		sub_v3_v3v3(x_next, x, step);
+		target_project_tri_clamp(x_next);
+	}
+
+	return true;
+}
+
+static bool target_project_solve_point_tri(
+        const float *vtri_co[3], const float vtri_no[3][3], const float point_co[3],
+        const float hit_co[3], float hit_dist_sq,
+        float r_hit_co[3], float r_hit_no[3])
+{
+	float x[3], tmp[3];
+	float dist = sqrtf(hit_dist_sq);
+	float epsilon = dist * 1.0e-5f;
+
+	CLAMP_MIN(epsilon, 1.0e-5f);
+
+	/* Initial solution vector: barycentric weights plus distance along normal. */
+	interp_weights_tri_v3(x, UNPACK3(vtri_co), hit_co);
+
+	interp_v3_v3v3v3(r_hit_no, UNPACK3(vtri_no), x);
+	sub_v3_v3v3(tmp, point_co, hit_co);
+
+	x[2] = (dot_v3v3(tmp, r_hit_no) < 0) ? -dist : dist;
+
+	/* Solve the equations iteratively. */
+	TargetProjectTriData tri_data = { .vtri_co = vtri_co, .vtri_no = vtri_no, .point_co = point_co };
+
+	sub_v3_v3v3(tri_data.n0_minus_n2, vtri_no[0], vtri_no[2]);
+	sub_v3_v3v3(tri_data.n1_minus_n2, vtri_no[1], vtri_no[2]);
+	sub_v3_v3v3(tri_data.c0_minus_c2, vtri_co[0], vtri_co[2]);
+	sub_v3_v3v3(tri_data.c1_minus_c2, vtri_co[1], vtri_co[2]);
+
+	target_project_tri_clamp(x);
+
+#ifdef TRACE_TARGET_PROJECT
+	const bool trace = true;
+#else
+	const bool trace = false;
+#endif
+
+	bool ok = BLI_newton3d_solve(target_project_tri_deviation, target_project_tri_jacobian, target_project_tri_correct, &tri_data, epsilon, 20, trace, x, x);
+
+	if (ok) {
+		copy_v3_v3(r_hit_co, tri_data.co_interp);
+		copy_v3_v3(r_hit_no, tri_data.no_interp);
+
+		return true;
+	}
+
+	return false;
+}
+
+static bool update_hit(BVHTreeNearest *nearest, int index, const float co[3], const float hit_co[3], const float hit_no[3])
+{
+	float dist_sq = len_squared_v3v3(hit_co, co);
+
+	if (dist_sq < nearest->dist_sq) {
+#ifdef TRACE_TARGET_PROJECT
+		printf("===> %d (%.3f,%.3f,%.3f) %g < %g\n", index, UNPACK3(hit_co), dist_sq, nearest->dist_sq);
+#endif
+		nearest->index = index;
+		nearest->dist_sq = dist_sq;
+		copy_v3_v3(nearest->co, hit_co);
+		normalize_v3_v3(nearest->no, hit_no);
+		return true;
+	}
+
+	return false;
+}
+
+/* Target projection on a non-manifold boundary edge - treats it like an infinitely thin cylinder. */
+static void target_project_edge(const ShrinkwrapTreeData *tree, int index, const float co[3], BVHTreeNearest *nearest, int eidx)
+{
+	const BVHTreeFromMesh *data = &tree->treeData;
+	const MEdge *edge = &tree->mesh->medge[eidx];
+	const float *vedge_co[2] = { data->vert[edge->v1].co, data->vert[edge->v2].co };
+
+#ifdef TRACE_TARGET_PROJECT
+	printf("EDGE %d (%.3f,%.3f,%.3f) (%.3f,%.3f,%.3f)\n", eidx, UNPACK3(vedge_co[0]), UNPACK3(vedge_co[1]));
+#endif
+
+	/* Retrieve boundary vertex IDs */
+	const int *vert_boundary_id = tree->boundary->vert_boundary_id;
+	int bid1 = vert_boundary_id[edge->v1], bid2 = vert_boundary_id[edge->v2];
+
+	if (bid1 < 0 || bid2 < 0) {
+		return;
+	}
+
+	/* Retrieve boundary vertex normals and align them to direction. */
+	const ShrinkwrapBoundaryVertData *boundary_verts = tree->boundary->boundary_verts;
+	float vedge_dir[2][3], dir[3];
+
+	copy_v3_v3(vedge_dir[0], boundary_verts[bid1].normal_plane);
+	copy_v3_v3(vedge_dir[1], boundary_verts[bid2].normal_plane);
+
+	sub_v3_v3v3(dir, vedge_co[1], vedge_co[0]);
+
+	if (dot_v3v3(boundary_verts[bid1].direction, dir) < 0) {
+		negate_v3(vedge_dir[0]);
+	}
+	if (dot_v3v3(boundary_verts[bid2].direction, dir) < 0) {
+		negate_v3(vedge_dir[1]);
+	}
+
+	/* Solve a quadratic equation: lerp(d0,d1,x) * (co - lerp(v0,v1,x)) = 0 */
+	float d0v0 = dot_v3v3(vedge_dir[0], vedge_co[0]), d0v1 = dot_v3v3(vedge_dir[0], vedge_co[1]);
+	float d1v0 = dot_v3v3(vedge_dir[1], vedge_co[0]), d1v1 = dot_v3v3(vedge_dir[1], vedge_co[1]);
+	float d0co = dot_v3v3(vedge_dir[0], co);
+
+	float a = d0v1 - d0v0 + d1v0 - d1v1;
+	float b = 2 * d0v0 - d0v1 - d0co - d1v0 + dot_v3v3(vedge_dir[1], co);
+	float c = d0co - d0v0;
+	float det = b * b - 4 * a * c;
+
+	if (det >= 0) {
+		const float epsilon = 1e-6f;
+		float sdet = sqrtf(det);
+		float hit_co[3], hit_no[3];
+
+		for (int i = (det > 0 ? 2 : 0); i >= 0; i -= 2) {
+			float x = (- b + ((float)i - 1) * sdet) / (2 * a);
+
+			if (x >= -epsilon && x <= 1.0f + epsilon) {
+				CLAMP(x, 0, 1);
+
+				float vedge_no[2][3];
+				normal_short_to_float_v3(vedge_no[0], data->vert[edge->v1].no);
+				normal_short_to_float_v3(vedge_no[1], data->vert[edge->v2].no);
+
+				interp_v3_v3v3(hit_co, vedge_co[0], vedge_co[1], x);
+				interp_v3_v3v3(hit_no, vedge_no[0], vedge_no[1], x);
+
+				update_hit(nearest, index, co, hit_co, hit_no);
+			}
+		}
+	}
+}
+
+/* Target normal projection BVH callback - based on mesh_looptri_nearest_point. */
+static void mesh_looptri_target_project(void *userdata, int index, const float co[3], BVHTreeNearest *nearest)
+{
+	const ShrinkwrapTreeData *tree = (ShrinkwrapTreeData *)userdata;
+	const BVHTreeFromMesh *data = &tree->treeData;
+	const MLoopTri *lt = &data->looptri[index];
+	const MLoop *loop[3] = { &data->loop[lt->tri[0]], &data->loop[lt->tri[1]], &data->loop[lt->tri[2]] };
+	const MVert *vtri[3] = { &data->vert[loop[0]->v], &data->vert[loop[1]->v], &data->vert[loop[2]->v] };
+	const float *vtri_co[3] = { vtri[0]->co, vtri[1]->co, vtri[2]->co };
+	float raw_hit_co[3], hit_co[3], hit_no[3], dist_sq, vtri_no[3][3];
+
+	/* First find the closest point and bail out if it's worse than the current solution. */
+	closest_on_tri_to_point_v3(raw_hit_co, co, UNPACK3(vtri_co));
+	dist_sq = len_squared_v3v3(co, raw_hit_co);
+
+#ifdef TRACE_TARGET_PROJECT
+	printf("TRIANGLE %d (%.3f,%.3f,%.3f) (%.3f,%.3f,%.3f) (%.3f,%.3f,%.3f) %g %g\n",
+	       index, UNPACK3(vtri_co[0]), UNPACK3(vtri_co[1]), UNPACK3(vtri_co[2]), dist_sq, nearest->dist_sq);
+#endif
+
+	if (dist_sq >= nearest->dist_sq)
+		return;
+
+	/* Decode normals */
+	normal_short_to_float_v3(vtri_no[0], vtri[0]->no);
+	normal_short_to_float_v3(vtri_no[1], vtri[1]->no);
+	normal_short_to_float_v3(vtri_no[2], vtri[2]->no);
+
+	/* Solve the equations for the triangle */
+	if (target_project_solve_point_tri(vtri_co, vtri_no, co, raw_hit_co, dist_sq, hit_co, hit_no)) {
+		update_hit(nearest, index, co, hit_co, hit_no);
+	}
+	/* Boundary edges */
+	else if (tree->boundary && BLI_BITMAP_TEST(tree->boundary->looptri_has_boundary, index)) {
+		const BLI_bitmap *is_boundary = tree->boundary->edge_is_boundary;
+		int edges[3];
+
+		BKE_mesh_looptri_get_real_edges(tree->mesh, lt, edges);
+
+		for (int i = 0; i < 3; i++) {
+			if (edges[i] >= 0 && BLI_BITMAP_TEST(is_boundary, edges[i])) {
+				target_project_edge(tree, index, co, nearest, edges[i]);
+			}
+		}
+	}
+}
+
+/*
+ * Maps the point to the nearest surface, either by simple nearest, or by target normal projection.
+ */
+void BKE_shrinkwrap_find_nearest_surface(struct ShrinkwrapTreeData *tree, BVHTreeNearest *nearest, float co[3], int type)
+{
+	BVHTreeFromMesh *treeData = &tree->treeData;
+
+	if (type == MOD_SHRINKWRAP_TARGET_PROJECT) {
+#ifdef TRACE_TARGET_PROJECT
+		printf("====== TARGET PROJECT START ======\n");
+#endif
+
+		BLI_bvhtree_find_nearest_ex(tree->bvh, co, nearest, mesh_looptri_target_project, tree, BVH_NEAREST_OPTIMAL_ORDER);
+
+#ifdef TRACE_TARGET_PROJECT
+		printf("====== TARGET PROJECT END: %d %g ======\n", nearest->index, nearest->dist_sq);
+#endif
+
+		if (nearest->index < 0) {
+			/* fallback to simple nearest */
+			BLI_bvhtree_find_nearest(tree->bvh, co, nearest, treeData->nearest_callback, treeData);
+		}
+	}
+	else {
+		BLI_bvhtree_find_nearest(tree->bvh, co, nearest, treeData->nearest_callback, treeData);
+	}
+}
+
+/*
  * Shrinkwrap moving vertexs to the nearest surface point on the target
  *
  * it builds a BVHTree from the target mesh and then performs a
@@ -536,7 +1058,6 @@ static void shrinkwrap_calc_nearest_surface_point_cb_ex(
 	ShrinkwrapCalcCBData *data = userdata;
 
 	ShrinkwrapCalcData *calc = data->calc;
-	BVHTreeFromMesh *treeData = &data->tree->treeData;
 	BVHTreeNearest *nearest = tls->userdata_chunk;
 
 	float *co = calc->vertexCos[i];
@@ -565,12 +1086,20 @@ static void shrinkwrap_calc_nearest_surface_point_cb_ex(
 	 * If we already had an hit before.. we assume this vertex is going to have a close hit to that other vertex
 	 * so we can initiate the "nearest.dist" with the expected value to that last hit.
 	 * This will lead in pruning of the search tree. */
-	if (nearest->index != -1)
-		nearest->dist_sq = len_squared_v3v3(tmp_co, nearest->co);
+	if (nearest->index != -1) {
+		if (calc->smd->shrinkType == MOD_SHRINKWRAP_TARGET_PROJECT) {
+			/* Heuristic doesn't work because of additional restrictions. */
+			nearest->index = -1;
+			nearest->dist_sq = FLT_MAX;
+		}
+		else {
+			nearest->dist_sq = len_squared_v3v3(tmp_co, nearest->co);
+		}
+	}
 	else
 		nearest->dist_sq = FLT_MAX;
 
-	BLI_bvhtree_find_nearest(treeData->tree, tmp_co, nearest, treeData->nearest_callback, treeData);
+	BKE_shrinkwrap_find_nearest_surface(data->tree, nearest, tmp_co, calc->smd->shrinkType);
 
 	/* Found the nearest vertex */
 	if (nearest->index != -1) {
@@ -838,6 +1367,7 @@ void shrinkwrapModifier_deform(ShrinkwrapModifierData *smd, struct Scene *scene,
 
 		switch (smd->shrinkType) {
 			case MOD_SHRINKWRAP_NEAREST_SURFACE:
+			case MOD_SHRINKWRAP_TARGET_PROJECT:
 				TIMEIT_BENCH(shrinkwrap_calc_nearest_surface_point(&calc), deform_surface);
 				break;