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

1 files changed, 536 insertions, 6 deletions

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;