Transform Snap: Optimize edge-snap using BVH tree

changes in BLI_kdopbvh:

- `BLI_bvhtree_find_nearest_to_ray` now takes is_ray_normalized and scale argument.
- `BLI_bvhtree_find_nearest_to_ray_angle` has been added (use for perspective view).

changes in BLI_bvhutils:

- `bvhtree_from_editmesh_edges_ex` was added.

changes in math_geom:

- `dist_squared_ray_to_seg_v3` was added.

other changes:

- `do_ray_start_correction` is no longer necessary to snap to verts.
- the way in which the test of depth was done before is being simulated in callbacks.

7 files changed, 1209 insertions, 256 deletions

diff --git a/source/blender/blenkernel/BKE_bvhutils.h b/source/blender/blenkernel/BKE_bvhutils.h
index 7ec91b006b2..80b248a932a 100644
--- a/source/blender/blenkernel/BKE_bvhutils.h
+++ b/source/blender/blenkernel/BKE_bvhutils.h
@@ -122,6 +122,14 @@ BVHTree *bvhtree_from_mesh_verts_ex(
         const bool vert_allocated, const BLI_bitmap *mask, int verts_num_active,
         float epsilon, int tree_type, int axis);
 
+BVHTree *bvhtree_from_editmesh_edges(
+        BVHTreeFromEditMesh *data, struct BMEditMesh *em,
+        float epsilon, int tree_type, int axis);
+BVHTree *bvhtree_from_editmesh_edges_ex(
+        BVHTreeFromEditMesh *data, struct BMEditMesh *em,
+        const BLI_bitmap *edges_mask, int edges_num_active,
+        float epsilon, int tree_type, int axis);
+
 BVHTree *bvhtree_from_mesh_edges(
         struct BVHTreeFromMesh *data, struct DerivedMesh *mesh,
         float epsilon, int tree_type, int axis);
diff --git a/source/blender/blenkernel/intern/bvhutils.c b/source/blender/blenkernel/intern/bvhutils.c
index a3cfe3f80b4..264d87b86f3 100644
--- a/source/blender/blenkernel/intern/bvhutils.c
+++ b/source/blender/blenkernel/intern/bvhutils.c
@@ -590,6 +590,77 @@ BVHTree *bvhtree_from_mesh_verts_ex(
 /** \name Edge Builder
  * \{ */
 
+static BVHTree *bvhtree_from_editmesh_edges_create_tree(
+        float epsilon, int tree_type, int axis,
+        BMEditMesh *em, const int edges_num,
+        const BLI_bitmap *edges_mask, int edges_num_active)
+{
+	BVHTree *tree = NULL;
+	int i;
+	BM_mesh_elem_table_ensure(em->bm, BM_EDGE);
+	if (edges_mask) {
+		BLI_assert(IN_RANGE_INCL(edges_num_active, 0, edges_num));
+	}
+	else {
+		edges_num_active = edges_num;
+	}
+
+	tree = BLI_bvhtree_new(edges_num_active, epsilon, tree_type, axis);
+
+	if (tree) {
+		BMIter iter;
+		BMEdge *eed;
+		BM_ITER_MESH_INDEX (eed, &iter, em->bm, BM_EDGES_OF_MESH, i) {
+			if (edges_mask && !BLI_BITMAP_TEST_BOOL(edges_mask, i)) {
+				continue;
+			}
+			float co[2][3];
+			copy_v3_v3(co[0], eed->v1->co);
+			copy_v3_v3(co[1], eed->v2->co);
+
+			BLI_bvhtree_insert(tree, i, co[0], 2);
+		}
+		BLI_assert(BLI_bvhtree_get_size(tree) == edges_num_active);
+		BLI_bvhtree_balance(tree);
+	}
+
+	return tree;
+}
+
+/* Builds a bvh tree where nodes are the edges of the given em */
+BVHTree *bvhtree_from_editmesh_edges_ex(
+        BVHTreeFromEditMesh *data, BMEditMesh *em,
+        const BLI_bitmap *edges_mask, int edges_num_active,
+        float epsilon, int tree_type, int axis)
+{
+	int edge_num = em->bm->totedge;
+
+	BVHTree *tree = bvhtree_from_editmesh_edges_create_tree(
+	        epsilon, tree_type, axis,
+	        em, edge_num, edges_mask, edges_num_active);
+
+	if (tree) {
+		memset(data, 0, sizeof(*data));
+		data->tree = tree;
+		data->em = em;
+		data->nearest_callback = NULL;  /* TODO */
+		data->raycast_callback = NULL;  /* TODO */
+		/* TODO: not urgent however since users currently define own callbacks */
+		data->nearest_to_ray_callback = NULL;
+	}
+
+	return tree;
+}
+BVHTree *bvhtree_from_editmesh_edges(
+        BVHTreeFromEditMesh *data, BMEditMesh *em,
+        float epsilon, int tree_type, int axis)
+{
+	return bvhtree_from_editmesh_edges_ex(
+	        data, em,
+	        NULL, -1,
+	        epsilon, tree_type, axis);
+}
+
 /* Builds a bvh tree where nodes are the edges of the given dm */
 BVHTree *bvhtree_from_mesh_edges(
         BVHTreeFromMesh *data, DerivedMesh *dm,
diff --git a/source/blender/blenlib/BLI_kdopbvh.h b/source/blender/blenlib/BLI_kdopbvh.h
index fb8c2520e67..91d39801645 100644
--- a/source/blender/blenlib/BLI_kdopbvh.h
+++ b/source/blender/blenlib/BLI_kdopbvh.h
@@ -96,7 +96,8 @@ typedef void (*BVHTree_NearestPointCallback)(void *userdata, int index, const fl
 typedef void (*BVHTree_RayCastCallback)(void *userdata, int index, const BVHTreeRay *ray, BVHTreeRayHit *hit);
 
 /* callback must update nearest in case it finds a nearest result */
-typedef void (*BVHTree_NearestToRayCallback)(void *userdata, int index, const BVHTreeRay *ray, BVHTreeNearest *nearest);
+typedef void (*BVHTree_NearestToRayCallback)(void *userdata, const float ray_co[3], const float ray_dir[3],
+                                             const float scale[3], int index, BVHTreeNearest *nearest);
 
 /* callback to check if 2 nodes overlap (use thread if intersection results need to be stored) */
 typedef bool (*BVHTree_OverlapCallback)(void *userdata, int index_a, int index_b, int thread);
@@ -142,8 +143,16 @@ int BLI_bvhtree_find_nearest(
         BVHTree *tree, const float co[3], BVHTreeNearest *nearest,
         BVHTree_NearestPointCallback callback, void *userdata);
 
+int BLI_bvhtree_find_nearest_to_ray_angle(
+        BVHTree *tree, const float co[3], const float dir[3],
+        const bool ray_is_normalized, const float scale[3],
+        BVHTreeNearest *nearest,
+        BVHTree_NearestToRayCallback callback, void *userdata);
+
 int BLI_bvhtree_find_nearest_to_ray(
-        BVHTree *tree, const float co[3], const float dir[3], BVHTreeNearest *nearest,
+        BVHTree *tree, const float co[3], const float dir[3],
+        const bool ray_is_normalized, const float scale[3],
+        BVHTreeNearest *nearest,
         BVHTree_NearestToRayCallback callback, void *userdata);
 
 int BLI_bvhtree_ray_cast_ex(
diff --git a/source/blender/blenlib/BLI_math_geom.h b/source/blender/blenlib/BLI_math_geom.h
index 54b824c91ac..84a25f533bf 100644
--- a/source/blender/blenlib/BLI_math_geom.h
+++ b/source/blender/blenlib/BLI_math_geom.h
@@ -115,6 +115,13 @@ float dist_signed_squared_to_corner_v3v3v3(
         const float p[3],
         const float v1[3], const float v2[3], const float v3[3],
         const float axis_ref[3]);
+float dist_squared_to_ray_v3(
+        const float ray_origin[3], const float ray_direction[3],
+        const float co[3], float *r_depth);
+float dist_squared_ray_to_seg_v3(
+        const float ray_origin[3], const float ray_direction[3],
+        const float v0[3], const float v1[3],
+        float r_point[3], float *r_depth);
 float closest_to_line_v2(float r_close[2], const float p[2], const float l1[2], const float l2[2]);
 float closest_to_line_v3(float r_close[3], const float p[3], const float l1[3], const float l2[3]);
 void closest_to_line_segment_v2(float r_close[2], const float p[2], const float l1[2], const float l2[2]);
diff --git a/source/blender/blenlib/intern/BLI_kdopbvh.c b/source/blender/blenlib/intern/BLI_kdopbvh.c
index 6cef1924e33..51b79d1fcd2 100644
--- a/source/blender/blenlib/intern/BLI_kdopbvh.c
+++ b/source/blender/blenlib/intern/BLI_kdopbvh.c
@@ -163,12 +163,23 @@ typedef struct BVHNearestRayData {
 	BVHTree *tree;
 	BVHTree_NearestToRayCallback callback;
 	void    *userdata;
-	BVHTreeRay ray;
 
-	struct NearestRayToAABB_Precalc nearest_precalc;
+	struct {
+		bool sign[3];
+		float origin[3];
+		float direction[3];
+
+		float direction_scaled_square[3];
+		float inv_dir[3];
+
+		float cdot_axis[3];
+	} ray;
 
 	bool pick_smallest[3];
+
 	BVHTreeNearest nearest;
+
+	float scale[3];
 } BVHNearestRayData;
 
 /** \} */
@@ -1889,58 +1900,374 @@ void BLI_bvhtree_ray_cast_all(
 
 /* -------------------------------------------------------------------- */
 
-/** \name BLI_bvhtree_find_nearest_to_ray
+/** \name BLI_bvhtree_find_nearest_to_ray functions
  *
  * \{ */
 
+static void dist_squared_ray_to_aabb_scaled_v3_precalc(
+        BVHNearestRayData *data,
+        const float ray_origin[3], const float ray_direction[3],
+        const bool ray_is_normalized, const float scale[3])
+{
+	if (scale) {
+		copy_v3_v3(data->scale, scale);
+	}
+	else {
+		copy_v3_fl(data->scale, 1.0f);
+	}
+	/* un-normalize ray */
+	if (ray_is_normalized && scale &&
+	   (data->scale[0] != 1.0f || data->scale[1] != 1.0f || data->scale[2] != 1.0f))
+	{
+		data->ray.direction[0] = ray_direction[0] * data->scale[0];
+		data->ray.direction[1] = ray_direction[1] * data->scale[1];
+		data->ray.direction[2] = ray_direction[2] * data->scale[2];
+
+		mul_v3_v3fl(data->ray.direction, ray_direction, 1 / len_v3(data->ray.direction));
+	}
+	else {
+		copy_v3_v3(data->ray.direction, ray_direction);
+	}
+
+	float dir_sq[3];
+
+	for (int i = 0; i < 3; i++) {
+		data->ray.origin[i] = ray_origin[i];
+		data->ray.inv_dir[i] = (data->ray.direction[i] != 0.0f) ?
+		                       (1.0f / data->ray.direction[i]) : FLT_MAX;
+		/* It has to be in function of `ray.inv_dir`,
+		 * since the division of 1 by 0.0f, can be -inf or +inf */
+		data->ray.sign[i] = (data->ray.inv_dir[i] < 0.0f);
+
+		data->ray.direction_scaled_square[i] = data->ray.direction[i] * data->scale[i];
+
+		dir_sq[i] = SQUARE(data->ray.direction_scaled_square[i]);
+
+		data->ray.direction_scaled_square[i] *= data->scale[i];
+	}
+
+	/* `diag_sq` Length square of each face diagonal */
+	float diag_sq[3] = {
+		dir_sq[1] + dir_sq[2],
+		dir_sq[0] + dir_sq[2],
+		dir_sq[0] + dir_sq[1],
+	};
+
+	data->ray.cdot_axis[0] = (diag_sq[0] != 0.0f) ? data->ray.direction[0] / diag_sq[0] : FLT_MAX;
+	data->ray.cdot_axis[1] = (diag_sq[1] != 0.0f) ? data->ray.direction[1] / diag_sq[1] : FLT_MAX;
+	data->ray.cdot_axis[2] = (diag_sq[2] != 0.0f) ? data->ray.direction[2] / diag_sq[2] : FLT_MAX;
+}
+
+/**
+ * Returns the squared distance from a ray to a bound-box `AABB`.
+ * It is based on `fast_ray_nearest_hit` solution to obtain
+ * the coordinates of the nearest edge of Bound Box to the ray
+ */
+MINLINE float dist_squared_ray_to_aabb_scaled_v3__impl(
+        const BVHNearestRayData *data,
+        const float bv[6], float *r_depth_sq, bool r_axis_closest[3])
+{
+
+	/* `tmin` is a vector that has the smaller distances to each of the
+	 * infinite planes of the `AABB` faces (hit in nearest face X plane,
+	 * nearest face Y plane and nearest face Z plane) */
+	float local_bvmin[3], local_bvmax[3];
+
+	if (data->ray.sign[0]) {
+		local_bvmin[0] = bv[1];
+		local_bvmax[0] = bv[0];
+	}
+	else {
+		local_bvmin[0] = bv[0];
+		local_bvmax[0] = bv[1];
+	}
+
+	if (data->ray.sign[1]) {
+		local_bvmin[1] = bv[3];
+		local_bvmax[1] = bv[2];
+	}
+	else {
+		local_bvmin[1] = bv[2];
+		local_bvmax[1] = bv[3];
+	}
+
+	if (data->ray.sign[2]) {
+		local_bvmin[2] = bv[5];
+		local_bvmax[2] = bv[4];
+	}
+	else {
+		local_bvmin[2] = bv[4];
+		local_bvmax[2] = bv[5];
+	}
+
+	sub_v3_v3(local_bvmin, data->ray.origin);
+	sub_v3_v3(local_bvmax, data->ray.origin);
+
+	const float tmin[3] = {
+		local_bvmin[0] * data->ray.inv_dir[0],
+		local_bvmin[1] * data->ray.inv_dir[1],
+		local_bvmin[2] * data->ray.inv_dir[2],
+	};
+
+	/* `tmax` is a vector that has the longer distances to each of the
+	 * infinite planes of the `AABB` faces (hit in farthest face X plane,
+	 * farthest face Y plane and farthest face Z plane) */
+	const float tmax[3] = {
+		local_bvmax[0] * data->ray.inv_dir[0],
+		local_bvmax[1] * data->ray.inv_dir[1],
+		local_bvmax[2] * data->ray.inv_dir[2],
+	};
+	/* `v1` and `v3` is be the coordinates of the nearest `AABB` edge to the ray*/
+	float v1[3], v2[3];
+	/* `rtmin` is the highest value of the smaller distances. == max_axis_v3(tmin)
+	* `rtmax` is the lowest value of longer distances. == min_axis_v3(tmax)*/
+	float rtmin, rtmax, mul;
+	/* `main_axis` is the axis equivalent to edge close to the ray */
+	int main_axis;
+
+	r_axis_closest[0] = false;
+	r_axis_closest[1] = false;
+	r_axis_closest[2] = false;
+
+	/* *** min_axis_v3(tmax) *** */
+	if ((tmax[0] <= tmax[1]) && (tmax[0] <= tmax[2])) {
+		// printf("# Hit in X %s\n", data->sign[0] ? "min", "max");
+		rtmax = tmax[0];
+		v1[0] = v2[0] = local_bvmax[0];
+		mul = local_bvmax[0] * data->ray.direction_scaled_square[0];
+		main_axis = 3;
+		r_axis_closest[0] = data->ray.sign[0];
+	}
+	else if ((tmax[1] <= tmax[0]) && (tmax[1] <= tmax[2])) {
+		// printf("# Hit in Y %s\n", data->sign[1] ? "min", "max");
+		rtmax = tmax[1];
+		v1[1] = v2[1] = local_bvmax[1];
+		mul = local_bvmax[1] * data->ray.direction_scaled_square[1];
+		main_axis = 2;
+		r_axis_closest[1] = data->ray.sign[1];
+	}
+	else {
+		// printf("# Hit in Z %s\n", data->sign[2] ? "min", "max");
+		rtmax = tmax[2];
+		v1[2] = v2[2] = local_bvmax[2];
+		mul = local_bvmax[2] * data->ray.direction_scaled_square[2];
+		main_axis = 1;
+		r_axis_closest[2] = data->ray.sign[2];
+	}
+
+	/* *** max_axis_v3(tmin) *** */
+	if ((tmin[0] >= tmin[1]) && (tmin[0] >= tmin[2])) {
+		// printf("# To X %s\n", data->sign[0] ? "max", "min");
+		rtmin = tmin[0];
+		v1[0] = v2[0] = local_bvmin[0];
+		mul += local_bvmin[0] * data->ray.direction_scaled_square[0];
+		main_axis -= 3;
+		r_axis_closest[0] = !data->ray.sign[0];
+	}
+	else if ((tmin[1] >= tmin[0]) && (tmin[1] >= tmin[2])) {
+		// printf("# To Y %s\n", data->sign[1] ? "max", "min");
+		rtmin = tmin[1];
+		v1[1] = v2[1] = local_bvmin[1];
+		mul += local_bvmin[1] * data->ray.direction_scaled_square[1];
+		main_axis -= 1;
+		r_axis_closest[1] = !data->ray.sign[1];
+	}
+	else {
+		// printf("# To Z %s\n", data->sign[2] ? "max", "min");
+		rtmin = tmin[2];
+		v1[2] = v2[2] = local_bvmin[2];
+		mul += local_bvmin[2] * data->ray.direction_scaled_square[2];
+		main_axis -= 2;
+		r_axis_closest[2] = !data->ray.sign[2];
+	}
+	/* *** end min/max axis *** */
+
+	if (main_axis < 0)
+		main_axis += 3;
+
+	/* if rtmin < rtmax, ray intersect `AABB` */
+	if (rtmin <= rtmax) {
+#ifdef IGNORE_BEHIND_RAY
+		/* `if rtmax < depth_min`, the whole `AABB` is behind us */
+		if (rtmax < min_depth) {
+			return fallback;
+		}
+#endif
+		const float proj = rtmin * data->ray.direction[main_axis];
+
+		if (data->ray.sign[main_axis])
+			r_axis_closest[main_axis] = (proj - local_bvmax[main_axis]) < (local_bvmin[main_axis] - proj);
+		else
+			r_axis_closest[main_axis] = (proj - local_bvmin[main_axis]) < (local_bvmax[main_axis] - proj);
+
+		//if (r_depth_sq)
+		//	*r_depth_sq = SQUARE(rtmin);
+
+		return 0.0f;
+	}
+#ifdef IGNORE_BEHIND_RAY
+	/* `if rtmin < depth_min`, the whole `AABB` is behing us */
+	else if (rtmin < min_depth) {
+		return fallback;
+	}
+#endif
+
+	if (data->ray.sign[main_axis]) {
+		v1[main_axis] = local_bvmax[main_axis];
+		v2[main_axis] = local_bvmin[main_axis];
+	}
+	else {
+		v1[main_axis] = local_bvmin[main_axis];
+		v2[main_axis] = local_bvmax[main_axis];
+	}
+	{
+		/* `proj` equals to nearest point on the ray closest to the edge `v1 v2` of the `AABB`. */
+		const float proj = mul * data->ray.cdot_axis[main_axis];
+		float depth_sq, r_point[3];
+		if (v1[main_axis] > proj) { /* the nearest point to the ray is the point v1 */
+			r_axis_closest[main_axis] = true;
+			/* `depth` is equivalent the distance of the the projection of v1 on the ray */
+			depth_sq = mul + data->ray.direction_scaled_square[main_axis] * v1[main_axis];
+
+			copy_v3_v3(r_point, v1);
+		}
+		else if (v2[main_axis] < proj) { /* the nearest point of the ray is the point v2 */
+			r_axis_closest[main_axis] = false;
+
+			depth_sq = mul + data->ray.direction_scaled_square[main_axis] * v2[main_axis];
+
+			copy_v3_v3(r_point, v2);
+		}
+		else {  /* the nearest point of the ray is on the edge of the `AABB`. */
+			r_axis_closest[main_axis] = (proj - v1[main_axis]) < (v2[main_axis] - proj);
+
+			depth_sq = mul + data->ray.direction_scaled_square[main_axis] * proj;
+#if 0
+			r_point[0] = main_axis == 0 ? proj : v2[0];
+			r_point[1] = main_axis == 1 ? proj : v2[1];
+			r_point[2] = main_axis == 2 ? proj : v2[2];
+#else
+			v2[main_axis] = proj;
+			copy_v3_v3(r_point, v2);
+#endif
+		}
+		depth_sq *= depth_sq;
+
+		if (r_depth_sq)
+			*r_depth_sq = depth_sq;
+
+		/* TODO: scale can be optional */
+		r_point[0] *= data->scale[0];
+		r_point[1] *= data->scale[1];
+		r_point[2] *= data->scale[2];
+
+		return len_squared_v3(r_point) - depth_sq;
+	}
+}
+
+/**
+ * <pre>
+ *  + r_point
+ *  |
+ *  | dist
+ *  |
+ *  +----depth----+orig <-- dir
+ *
+ * tangent = dist/depth
+ * </pre>
+ */
+static float calc_tangent_sq(BVHNearestRayData *data, BVHNode *node)
+{
+	float depth_sq;
+	const float dist_sq = dist_squared_ray_to_aabb_scaled_v3__impl(
+	        data, node->bv, &depth_sq, data->pick_smallest);
+
+	return (dist_sq != 0.0f) ? (dist_sq / depth_sq) : 0.0f;
+}
+
 static float calc_dist_sq_to_ray(BVHNearestRayData *data, BVHNode *node)
 {
-	const float *bv = node->bv;
-	const float bb_min[3] = {bv[0], bv[2], bv[4]};
-	const float bb_max[3] = {bv[1], bv[3], bv[5]};
-	return dist_squared_ray_to_aabb_v3(&data->nearest_precalc, bb_min, bb_max, data->pick_smallest);
+	return dist_squared_ray_to_aabb_scaled_v3__impl(
+	        data, node->bv, NULL,
+	        data->pick_smallest);
 }
 
-static void dfs_find_nearest_to_ray_dfs(BVHNearestRayData *data, BVHNode *node)
+static void dfs_find_lowest_tangent_dfs(BVHNearestRayData *data, BVHNode *node)
 {
 	if (node->totnode == 0) {
 		if (data->callback) {
-			data->callback(data->userdata, node->index, &data->ray, &data->nearest);
+			data->callback(data->userdata, data->ray.origin, data->ray.direction,
+			               data->scale, node->index, &data->nearest);
+		}
+		else {
+			data->nearest.index = node->index;
+			data->nearest.dist_sq = calc_tangent_sq(data, node);
+			/* TODO: return a value to the data->nearest.co
+			 * not urgent however since users currently define own callbacks */
+		}
+	}
+	else {
+		int i;
+		/* First pick the closest node to dive on */
+		if (data->pick_smallest[node->main_axis]) {
+			for (i = 0; i != node->totnode; i++) {
+				if (calc_tangent_sq(data, node->children[i]) < data->nearest.dist_sq) {
+					dfs_find_lowest_tangent_dfs(data, node->children[i]);
+				}
+			}
 		}
 		else {
-			const float dist_sq = calc_dist_sq_to_ray(data, node);
-			if (dist_sq != FLT_MAX) {  /* not an invalid ray */
-				data->nearest.index = node->index;
-				data->nearest.dist_sq = dist_sq;
-				/* TODO: return a value to the data->nearest.co
-				 * not urgent however since users currently define own callbacks */
+			for (i = node->totnode - 1; i >= 0; i--) {
+				if (calc_tangent_sq(data, node->children[i]) < data->nearest.dist_sq) {
+					dfs_find_lowest_tangent_dfs(data, node->children[i]);
+				}
 			}
 		}
 	}
+}
+
+static void dfs_find_nearest_to_ray_dfs(BVHNearestRayData *data, BVHNode *node)
+{
+	if (node->totnode == 0) {
+		if (data->callback) {
+			data->callback(data->userdata, data->ray.origin, data->ray.direction,
+			               data->scale, node->index, &data->nearest);
+		}
+		else {
+			data->nearest.index = node->index;
+			data->nearest.dist_sq = calc_dist_sq_to_ray(data, node);
+			/* TODO: return a value to the data->nearest.co
+			 * not urgent however since users currently define own callbacks */
+		}
+	}
 	else {
 		int i;
 		/* First pick the closest node to dive on */
 		if (data->pick_smallest[node->main_axis]) {
 			for (i = 0; i != node->totnode; i++) {
-				if (calc_dist_sq_to_ray(data, node->children[i]) >= data->nearest.dist_sq) {
-					continue;
+				if (calc_dist_sq_to_ray(data, node->children[i]) < data->nearest.dist_sq) {
+					dfs_find_nearest_to_ray_dfs(data, node->children[i]);
 				}
-				dfs_find_nearest_to_ray_dfs(data, node->children[i]);
 			}
 		}
 		else {
 			for (i = node->totnode - 1; i >= 0; i--) {
-				if (calc_dist_sq_to_ray(data, node->children[i]) >= data->nearest.dist_sq) {
-					continue;
+				if (calc_dist_sq_to_ray(data, node->children[i]) < data->nearest.dist_sq) {
+					dfs_find_nearest_to_ray_dfs(data, node->children[i]);
 				}
-				dfs_find_nearest_to_ray_dfs(data, node->children[i]);
 			}
 		}
 	}
 }
 
-int BLI_bvhtree_find_nearest_to_ray(
-        BVHTree *tree, const float co[3], const float dir[3], BVHTreeNearest *nearest,
+/**
+ * Returns the point whose tangent defined by the angle between the point and ray is the lowest
+ * nearest.dist_sq returns the angle's tangent
+ */
+int BLI_bvhtree_find_nearest_to_ray_angle(
+        BVHTree *tree, const float co[3], const float dir[3],
+        const bool ray_is_normalized, const float scale[3],
+        BVHTreeNearest *nearest,
         BVHTree_NearestToRayCallback callback, void *userdata)
 {
 	BVHNearestRayData data;
@@ -1951,11 +2278,46 @@ int BLI_bvhtree_find_nearest_to_ray(
 	data.callback = callback;
 	data.userdata = userdata;
 
-	copy_v3_v3(data.ray.origin, co);
-	copy_v3_v3(data.ray.direction, dir);
-	data.ray.radius = 0.0f;  /* unused here */
+	dist_squared_ray_to_aabb_scaled_v3_precalc(&data, co, dir, ray_is_normalized, scale);
+
+	if (nearest) {
+		memcpy(&data.nearest, nearest, sizeof(*nearest));
+	}
+	else {
+		data.nearest.index = -1;
+		data.nearest.dist_sq = FLT_MAX;
+	}
+
+	/* dfs search */
+	if (root) {
+		if (calc_tangent_sq(&data, root) < data.nearest.dist_sq)
+			dfs_find_lowest_tangent_dfs(&data, root);
+	}
+
+	/* copy back results */
+	if (nearest) {
+		memcpy(nearest, &data.nearest, sizeof(*nearest));
+	}
+
+	return data.nearest.index;
+}
+
+/* return the nearest point to ray */
+int BLI_bvhtree_find_nearest_to_ray(
+        BVHTree *tree, const float co[3], const float dir[3],
+        const bool ray_is_normalized, const float scale[3],
+        BVHTreeNearest *nearest,
+        BVHTree_NearestToRayCallback callback, void *userdata)
+{
+	BVHNearestRayData data;
+	BVHNode *root = tree->nodes[tree->totleaf];
+
+	data.tree = tree;
+
+	data.callback = callback;
+	data.userdata = userdata;
 
-	dist_squared_ray_to_aabb_v3_precalc(&data.nearest_precalc, co, dir);
+	dist_squared_ray_to_aabb_scaled_v3_precalc(&data, co, dir, ray_is_normalized, scale);
 
 	if (nearest) {
 		memcpy(&data.nearest, nearest, sizeof(*nearest));
diff --git a/source/blender/blenlib/intern/math_geom.c b/source/blender/blenlib/intern/math_geom.c
index 82f527942e8..40454a93ec8 100644
--- a/source/blender/blenlib/intern/math_geom.c
+++ b/source/blender/blenlib/intern/math_geom.c
@@ -572,6 +572,67 @@ float dist_signed_squared_to_corner_v3v3v3(
 	}
 }
 
+/**
+ * return the distance squared of a point to a ray.
+ */
+float dist_squared_to_ray_v3(
+        const float ray_origin[3], const float ray_direction[3],
+        const float co[3], float *r_depth)
+{
+	float dvec[3];
+	sub_v3_v3v3(dvec, co, ray_origin);
+	*r_depth = dot_v3v3(dvec, ray_direction);
+	return len_squared_v3(dvec) - SQUARE(*r_depth);
+}
+/**
+ * Find the closest point in a seg to a ray and return the distance squared.
+ * \param r_point : Is the point on segment closest to ray (or to ray_origin if the ray and the segment are parallel).
+ * \param depth: the distance of r_point projection on ray to the ray_origin.
+ */
+float dist_squared_ray_to_seg_v3(
+        const float ray_origin[3], const float ray_direction[3],
+        const float v0[3], const float v1[3],
+        float r_point[3], float *r_depth)
+{
+	float a[3], t[3], n[3], lambda;
+	sub_v3_v3v3(a, v1, v0);
+	sub_v3_v3v3(t, v0, ray_origin);
+	cross_v3_v3v3(n, a, ray_direction);
+	const float nlen = len_squared_v3(n);
+
+	/* if (nlen == 0.0f) the lines are parallel,
+	 * has no nearest point, only distance squared.*/
+	if (nlen == 0.0f) {
+		/* Calculate the distance to the point v0 then */
+		copy_v3_v3(r_point, v0);
+		*r_depth = dot_v3v3(t, ray_direction);
+	}
+	else {
+		float c[3], cray[3];
+		sub_v3_v3v3(c, n, t);
+		cross_v3_v3v3(cray, c, ray_direction);
+		lambda = dot_v3v3(cray, n) / nlen;
+		if (lambda <= 0) {
+			copy_v3_v3(r_point, v0);
+
+			*r_depth = dot_v3v3(t, ray_direction);
+		}
+		else if (lambda >= 1) {
+			copy_v3_v3(r_point, v1);
+
+			sub_v3_v3v3(t, v1, ray_origin);
+			*r_depth = dot_v3v3(t, ray_direction);
+		}
+		else {
+			madd_v3_v3v3fl(r_point, v0, a, lambda);
+
+			sub_v3_v3v3(t, r_point, ray_origin);
+			*r_depth = dot_v3v3(t, ray_direction);
+		}
+	}
+	return len_squared_v3(t) - SQUARE(*r_depth);
+}
+
 /* Adapted from "Real-Time Collision Detection" by Christer Ericson,
  * published by Morgan Kaufmann Publishers, copyright 2005 Elsevier Inc.
  * 
diff --git a/source/blender/editors/transform/transform_snap_object.c b/source/blender/editors/transform/transform_snap_object.c
index 0e1a65e073f..aa1dafabde7 100644
--- a/source/blender/editors/transform/transform_snap_object.c
+++ b/source/blender/editors/transform/transform_snap_object.c
@@ -68,13 +68,13 @@ typedef struct SnapObjectData {
 
 typedef struct SnapObjectData_Mesh {
 	SnapObjectData sd;
-	BVHTreeFromMesh *bvh_trees[2];
+	BVHTreeFromMesh *bvh_trees[3];
 
 } SnapObjectData_Mesh;
 
 typedef struct SnapObjectData_EditMesh {
 	SnapObjectData sd;
-	BVHTreeFromEditMesh *bvh_trees[2];
+	BVHTreeFromEditMesh *bvh_trees[3];
 
 } SnapObjectData_EditMesh;
 
@@ -223,8 +223,11 @@ static void raycast_all_cb(void *userdata, int index, const BVHTreeRay *ray, BVH
 
 static bool snapEdge(
         const ARegion *ar, const float v1co[3], const short v1no[3], const float v2co[3], const short v2no[3],
-        float obmat[4][4], float timat[3][3], const float mval_fl[2], float *dist_px,
-        const float ray_start[3], const float ray_start_local[3], const float ray_normal_local[3], float *ray_depth,
+        float obmat[4][4], float timat[3][3], const float mval_fl[2],
+        const float ray_start[3], const float ray_start_local[3], const float ray_normal_local[3],
+        /* read/write args */
+        float *ray_depth, float *dist_px,
+        /* return args */
         float r_loc[3], float r_no[3])
 {
 	float intersect[3] = {0, 0, 0}, ray_end[3], dvec[3];
@@ -313,8 +316,11 @@ static bool snapEdge(
 
 static bool snapVertex(
         const ARegion *ar, const float vco[3], const float vno[3],
-        float obmat[4][4], float timat[3][3], const float mval_fl[2], float *dist_px,
-        const float ray_start[3], const float ray_start_local[3], const float ray_normal_local[3], float *ray_depth,
+        float obmat[4][4], float timat[3][3], const float mval_fl[2],
+        const float ray_start[3], const float ray_start_local[3], const float ray_normal_local[3],
+        /* read/write args */
+        float *ray_depth, float *dist_px,
+        /* return args */
         float r_loc[3], float r_no[3])
 {
 	bool retval = false;
@@ -363,8 +369,11 @@ static bool snapVertex(
 
 static bool snapArmature(
         const ARegion *ar, Object *ob, bArmature *arm, float obmat[4][4],
-        const float mval[2], float *dist_px, const short snap_to,
-        const float ray_start[3], const float ray_normal[3], float *ray_depth,
+        const float mval[2], const short snap_to,
+        const float ray_start[3], const float ray_normal[3],
+        /* read/write args */
+        float *ray_depth, float *dist_px,
+        /* return args */
         float r_loc[3], float *UNUSED(r_no))
 {
 	float imat[4][4];
@@ -445,8 +454,11 @@ static bool snapArmature(
 
 static bool snapCurve(
         const ARegion *ar, Object *ob, Curve *cu, float obmat[4][4],
-        const float mval[2], float *dist_px, const short snap_to,
-        const float ray_start[3], const float ray_normal[3], float *ray_depth,
+        const float mval[2], const short snap_to,
+        const float ray_start[3], const float ray_normal[3],
+        /* read/write args */
+        float *ray_depth, float *dist_px,
+        /* return args */
         float r_loc[3], float *UNUSED(r_no))
 {
 	float imat[4][4];
@@ -481,24 +493,27 @@ static bool snapCurve(
 								break;
 							}
 							retval |= snapVertex(
-							        ar, nu->bezt[u].vec[1], NULL, obmat, NULL,  mval, dist_px,
-							        ray_start, ray_start_local, ray_normal_local, ray_depth,
+							        ar, nu->bezt[u].vec[1], NULL, obmat, NULL,  mval,
+							        ray_start, ray_start_local, ray_normal_local,
+							        ray_depth, dist_px,
 							        r_loc, NULL);
 							/* don't snap if handle is selected (moving), or if it is aligning to a moving handle */
 							if (!(nu->bezt[u].f1 & SELECT) &&
 							    !(nu->bezt[u].h1 & HD_ALIGN && nu->bezt[u].f3 & SELECT))
 							{
 								retval |= snapVertex(
-								        ar, nu->bezt[u].vec[0], NULL, obmat, NULL, mval, dist_px,
-								        ray_start, ray_start_local, ray_normal_local, ray_depth,
+								        ar, nu->bezt[u].vec[0], NULL, obmat, NULL, mval,
+								        ray_start, ray_start_local, ray_normal_local,
+								        ray_depth, dist_px,
 								        r_loc, NULL);
 							}
 							if (!(nu->bezt[u].f3 & SELECT) &&
 							    !(nu->bezt[u].h2 & HD_ALIGN && nu->bezt[u].f1 & SELECT))
 							{
 								retval |= snapVertex(
-								        ar, nu->bezt[u].vec[2], NULL, obmat, NULL, mval, dist_px,
-								        ray_start, ray_start_local, ray_normal_local, ray_depth,
+								        ar, nu->bezt[u].vec[2], NULL, obmat, NULL, mval,
+								        ray_start, ray_start_local, ray_normal_local,
+								        ray_depth, dist_px,
 								        r_loc, NULL);
 							}
 						}
@@ -508,8 +523,9 @@ static bool snapCurve(
 								break;
 							}
 							retval |= snapVertex(
-							        ar, nu->bp[u].vec, NULL, obmat, NULL, mval, dist_px,
-							        ray_start, ray_start_local, ray_normal_local, ray_depth,
+							        ar, nu->bp[u].vec, NULL, obmat, NULL, mval,
+							        ray_start, ray_start_local, ray_normal_local,
+							        ray_depth, dist_px,
 							        r_loc, NULL);
 						}
 					}
@@ -518,14 +534,16 @@ static bool snapCurve(
 						if (nu->pntsu > 1) {
 							if (nu->bezt) {
 								retval |= snapVertex(
-								        ar, nu->bezt[u].vec[1], NULL, obmat, NULL, mval, dist_px,
-								        ray_start, ray_start_local, ray_normal_local, ray_depth,
+								        ar, nu->bezt[u].vec[1], NULL, obmat, NULL, mval,
+								        ray_start, ray_start_local, ray_normal_local,
+								        ray_depth, dist_px,
 								        r_loc, NULL);
 							}
 							else {
 								retval |= snapVertex(
-								        ar, nu->bp[u].vec, NULL, obmat, NULL, mval, dist_px,
-								        ray_start, ray_start_local, ray_normal_local, ray_depth,
+								        ar, nu->bp[u].vec, NULL, obmat, NULL, mval,
+								        ray_start, ray_start_local, ray_normal_local,
+								        ray_depth, dist_px,
 								        r_loc, NULL);
 							}
 						}
@@ -543,8 +561,11 @@ static bool snapCurve(
 /* may extend later (for now just snaps to empty center) */
 static bool snapEmpty(
         const ARegion *ar, Object *ob, float obmat[4][4],
-        const float mval[2], float *dist_px, const short snap_to,
-        const float ray_start[3], const float ray_normal[3], float *ray_depth,
+        const float mval[2], const short snap_to,
+        const float ray_start[3], const float ray_normal[3],
+        /* read/write args */
+        float *ray_depth, float *dist_px,
+        /* return args */
         float r_loc[3], float *UNUSED(r_no))
 {
 	float imat[4][4];
@@ -583,8 +604,11 @@ static bool snapEmpty(
 
 static bool snapCamera(
         const ARegion *ar, Scene *scene, Object *object, float obmat[4][4],
-        const float mval[2], float *dist_px, const short snap_to,
-        const float ray_start[3], const float ray_normal[3], float *ray_depth,
+        const float mval[2], const short snap_to,
+        const float ray_start[3], const float ray_normal[3],
+        /* read/write args */
+        float *ray_depth, float *dist_px,
+        /* return args */
         float r_loc[3], float *UNUSED(r_no))
 {
 	float orig_camera_mat[4][4], orig_camera_imat[4][4], imat[4][4];
@@ -674,16 +698,166 @@ static int dm_looptri_to_poly_index(DerivedMesh *dm, const MLoopTri *lt)
 	return index_mp_to_orig ? index_mp_to_orig[lt->poly] : lt->poly;
 }
 
+struct NearestDM_Data {
+	void *bvhdata;
+	bool is_persp;
+	const float *ray_depth_range;
+
+	float *ray_depth;
+};
+
+static bool test_vert(
+        const float vco[3], const float vno[3], const float ray_co[3], const float ray_dir[3],
+        const float ray_depth_range[2], const float scale[3], const bool is_persp,
+        /* read/write args */
+        float *ray_depth, float *dist_to_ray_sq,
+        /* return args */
+        float r_co[3], float r_no[3])
+{
+	const float vco_sc[3] = {
+		vco[0] * scale[0],
+		vco[1] * scale[1],
+		vco[2] * scale[2],
+	};
+	const float co_sc[3] = {
+		ray_co[0] * scale[0],
+		ray_co[1] * scale[1],
+		ray_co[2] * scale[2],
+	};
+	const float dir_sc[3] = {
+		ray_dir[0] * scale[0],
+		ray_dir[1] * scale[1],
+		ray_dir[2] * scale[2],
+	};
+
+	float depth;
+	float dist_sq = dist_squared_to_ray_v3(co_sc, dir_sc, vco_sc, &depth);
+
+	if (depth < ray_depth_range[0])
+		return false;
+
+	if (is_persp)
+		dist_sq /= SQUARE(depth);
+
+	if ((dist_sq < *dist_to_ray_sq) && (depth < *ray_depth)) {
+		*dist_to_ray_sq = dist_sq;
+
+		copy_v3_v3(r_co, vco);
+
+		if (vno) {
+			copy_v3_v3(r_no, vno);
+		}
+
+		*ray_depth = depth;
+		return true;
+	}
+	return false;
+}
+
+static bool test_edge(
+        const float v1[3], const float v2[3], const float ray_co[3], const float ray_dir[3],
+        const float ray_depth_range[2], const float scale[3], const bool is_persp,
+        /* read/write args */
+        float *ray_depth, float *dist_to_ray_sq,
+        /* return args */
+        float r_co[3], float r_no[3])
+{
+	const float v1_sc[3] = {
+		v1[0] * scale[0],
+		v1[1] * scale[1],
+		v1[2] * scale[2],
+	};
+	const float v2_sc[3] = {
+		v2[0] * scale[0],
+		v2[1] * scale[1],
+		v2[2] * scale[2],
+	};
+	const float co_sc[3] = {
+		ray_co[0] * scale[0],
+		ray_co[1] * scale[1],
+		ray_co[2] * scale[2],
+	};
+	const float dir_sc[3] = {
+		ray_dir[0] * scale[0],
+		ray_dir[1] * scale[1],
+		ray_dir[2] * scale[2],
+	};
+
+	float tmp_co[3], depth;
+	float dist_sq = dist_squared_ray_to_seg_v3(co_sc, dir_sc, v1_sc, v2_sc, tmp_co, &depth);
+
+	if (depth < ray_depth_range[0])
+		return false;
+
+	if (is_persp)
+		dist_sq /= SQUARE(depth);
+
+	if ((dist_sq < *dist_to_ray_sq) && (depth < *ray_depth)) {
+		*dist_to_ray_sq = dist_sq;
+
+		tmp_co[0] /= scale[0];
+		tmp_co[1] /= scale[1];
+		tmp_co[2] /= scale[2];
+
+		copy_v3_v3(r_co, tmp_co);
+
+		if (r_no) {
+			sub_v3_v3v3(r_no, v1, v2);
+		}
+
+		*ray_depth = depth;
+		return true;
+	}
+	return false;
+}
+
+static void test_vert_depth_cb(
+        void *userdata, const float origin[3], const float dir[3],
+        const float scale[3], int index, BVHTreeNearest *nearest)
+{
+	struct NearestDM_Data *ndata = userdata;
+	const BVHTreeFromMesh *data = ndata->bvhdata;
+	const MVert *vert = data->vert + index;
+
+	if (test_vert(
+	        vert->co, NULL, origin, dir, ndata->ray_depth_range, scale, ndata->is_persp,
+	        ndata->ray_depth, &nearest->dist_sq,
+	        nearest->co, NULL))
+	{
+		normal_short_to_float_v3(nearest->no, vert->no);
+		nearest->index = index;
+	}
+}
+
+static void test_edge_depth_cb(
+        void *userdata, const float origin[3], const float dir[3],
+        const float scale[3], int index, BVHTreeNearest *nearest)
+{
+	struct NearestDM_Data *ndata = userdata;
+	const BVHTreeFromMesh *data = ndata->bvhdata;
+	const MVert *vert = data->vert;
+	const MEdge *edge = data->edge + index;
+
+	if (test_edge(
+	        vert[edge->v1].co, vert[edge->v2].co, origin, dir, ndata->ray_depth_range, scale, ndata->is_persp,
+	        ndata->ray_depth, &nearest->dist_sq,
+	        nearest->co, nearest->no))
+	{
+		nearest->index = index;
+	}
+}
+
 static bool snapDerivedMesh(
         SnapObjectContext *sctx,
-        Object *ob, DerivedMesh *dm, float obmat[4][4],
-        const float mval[2], float *dist_px, const short snap_to, bool do_bb,
-        const float ray_start[3], const float ray_normal[3], const float ray_origin[3],
-        float *ray_depth, unsigned int ob_index,
+        Object *ob, DerivedMesh *dm, float obmat[4][4], const unsigned int ob_index,
+        const short snap_to, bool do_bb,
+        const float ray_origin[3], const float ray_start[3], const float ray_normal[3], const float ray_depth_range[2],
+        /* read/write args */
+        float *ray_depth, float *dist_to_ray_sq, float *dist_px,
+        /* return args */
         float r_loc[3], float r_no[3], int *r_index,
         ListBase *r_hit_list)
 {
-	const ARegion *ar = sctx->v3d_data.ar;
 	bool retval = false;
 
 	if (snap_to == SCE_SNAP_MODE_FACE) {
@@ -703,10 +877,8 @@ static bool snapDerivedMesh(
 	}
 
 	{
-		const bool do_ray_start_correction = (
-		         ELEM(snap_to, SCE_SNAP_MODE_FACE, SCE_SNAP_MODE_VERTEX) &&
-		         (sctx->use_v3d && !((RegionView3D *)sctx->v3d_data.ar->regiondata)->is_persp));
-		bool need_ray_start_correction_init = do_ray_start_correction;
+		const bool is_persp = sctx->use_v3d && ((RegionView3D *)sctx->v3d_data.ar->regiondata)->is_persp;
+		bool need_ray_start_correction_init = (snap_to == SCE_SNAP_MODE_FACE) && sctx->use_v3d && !is_persp;
 
 		float imat[4][4];
 		float timat[3][3]; /* transpose inverse matrix for normals */
@@ -772,6 +944,9 @@ static bool snapDerivedMesh(
 			int tree_index = -1;
 			switch (snap_to) {
 				case SCE_SNAP_MODE_FACE:
+					tree_index = 2;
+					break;
+				case SCE_SNAP_MODE_EDGE:
 					tree_index = 1;
 					break;
 				case SCE_SNAP_MODE_VERTEX:
@@ -793,10 +968,8 @@ static bool snapDerivedMesh(
 			}
 		}
 		else {
-			if (ELEM(snap_to, SCE_SNAP_MODE_FACE, SCE_SNAP_MODE_VERTEX)) {
-				treedata = &treedata_stack;
-				memset(treedata, 0, sizeof(*treedata));
-			}
+			treedata = &treedata_stack;
+			memset(treedata, 0, sizeof(*treedata));
 		}
 
 		if (treedata && treedata->tree == NULL) {
@@ -804,53 +977,56 @@ static bool snapDerivedMesh(
 				case SCE_SNAP_MODE_FACE:
 					bvhtree_from_mesh_looptri(treedata, dm, 0.0f, 4, 6);
 					break;
+				case SCE_SNAP_MODE_EDGE:
+					bvhtree_from_mesh_edges(treedata, dm, 0.0f, 2, 6);
+					break;
 				case SCE_SNAP_MODE_VERTEX:
 					bvhtree_from_mesh_verts(treedata, dm, 0.0f, 2, 6);
 					break;
 			}
 		}
-
-		if (need_ray_start_correction_init) {
-			/* We *need* a reasonably valid len_diff in this case.
-			 * Use BHVTree to find the closest face from ray_start_local.
-			 */
-			if (treedata && treedata->tree != NULL) {
-				BVHTreeNearest nearest;
-				nearest.index = -1;
-				nearest.dist_sq = FLT_MAX;
-				/* Compute and store result. */
-				BLI_bvhtree_find_nearest(
-				            treedata->tree, ray_start_local, &nearest, treedata->nearest_callback, treedata);
-				if (nearest.index != -1) {
-					len_diff = sqrtf(nearest.dist_sq);
-				}
-			}
-		}
-		/* Only use closer ray_start in case of ortho view! In perspective one, ray_start may already
-		 * been *inside* boundbox, leading to snap failures (see T38409).
-		 * Note also ar might be null (see T38435), in this case we assume ray_start is ok!
-		 */
-		if (do_ray_start_correction) {
-			float ray_org_local[3];
-
-			copy_v3_v3(ray_org_local, ray_origin);
-			mul_m4_v3(imat, ray_org_local);
-
-			/* We pass a temp ray_start, set from object's boundbox, to avoid precision issues with very far
-			 * away ray_start values (as returned in case of ortho view3d), see T38358.
-			 */
-			len_diff -= local_scale;  /* make temp start point a bit away from bbox hit point. */
-			madd_v3_v3v3fl(ray_start_local, ray_org_local, ray_normal_local,
-			               len_diff - len_v3v3(ray_start_local, ray_org_local));
-			local_depth -= len_diff;
-		}
-		else {
-			len_diff = 0.0f;
-		}
-
 		switch (snap_to) {
 			case SCE_SNAP_MODE_FACE:
 			{
+				/* Only use closer ray_start in case of ortho view! In perspective one, ray_start may already
+				 * been *inside* boundbox, leading to snap failures (see T38409).
+				 * Note also ar might be null (see T38435), in this case we assume ray_start is ok!
+				 */
+				if (sctx->use_v3d && !is_persp) { /* do_ray_start_correction */
+					if (need_ray_start_correction_init) {
+						/* We *need* a reasonably valid len_diff in this case.
+						 * Use BHVTree to find the closest face from ray_start_local.
+						 */
+						if (treedata && treedata->tree != NULL) {
+							BVHTreeNearest nearest;
+							nearest.index = -1;
+							nearest.dist_sq = FLT_MAX;
+							/* Compute and store result. */
+							BLI_bvhtree_find_nearest(
+								treedata->tree, ray_start_local, &nearest, treedata->nearest_callback, treedata);
+							if (nearest.index != -1) {
+								float dvec[3];
+								sub_v3_v3v3(dvec, nearest.co, ray_start_local);
+								len_diff = dot_v3v3(dvec, ray_normal_local);
+							}
+						}
+					}
+					float ray_org_local[3];
+
+					copy_v3_v3(ray_org_local, ray_origin);
+					mul_m4_v3(imat, ray_org_local);
+
+					/* We pass a temp ray_start, set from object's boundbox, to avoid precision issues with very far
+					 * away ray_start values (as returned in case of ortho view3d), see T38358.
+					 */
+					len_diff -= local_scale;  /* make temp start point a bit away from bbox hit point. */
+					madd_v3_v3v3fl(ray_start_local, ray_org_local, ray_normal_local,
+					               len_diff + ray_depth_range[0]);
+					local_depth -= len_diff;
+				}
+				else {
+					len_diff = 0.0f;
+				}
 				if (r_hit_list) {
 					struct RayCastAll_Data data;
 
@@ -861,7 +1037,7 @@ static bool snapDerivedMesh(
 					data.len_diff = len_diff;
 					data.local_scale = local_scale;
 					data.ob = ob;
-					data.ob_uuid = ob_index,
+					data.ob_uuid = ob_index;
 					data.dm = dm;
 					data.hit_list = r_hit_list;
 					data.retval = retval;
@@ -907,40 +1083,90 @@ static bool snapDerivedMesh(
 			}
 			case SCE_SNAP_MODE_VERTEX:
 			{
+				float ray_org_local[3];
+
+				copy_v3_v3(ray_org_local, ray_origin);
+				mul_m4_v3(imat, ray_org_local);
+
 				BVHTreeNearest nearest;
 
 				nearest.index = -1;
-				nearest.dist_sq = local_depth * local_depth;
-				if (treedata->tree &&
+				nearest.dist_sq = *dist_to_ray_sq;
+
+				struct NearestDM_Data userdata;
+				userdata.bvhdata = treedata;
+				userdata.is_persp = is_persp;
+				userdata.ray_depth_range = ray_depth_range;
+				userdata.ray_depth = ray_depth;
+
+				float ob_scale[3];
+				mat4_to_size(ob_scale, obmat);
+
+				if (treedata->tree && (
+				    is_persp ?
+				    BLI_bvhtree_find_nearest_to_ray_angle(
+				            treedata->tree, ray_org_local, ray_normal_local,
+				            true, ob_scale, &nearest, test_vert_depth_cb, &userdata) :
 				    BLI_bvhtree_find_nearest_to_ray(
-				        treedata->tree, ray_start_local, ray_normal_local,
-				        &nearest, NULL, NULL) != -1)
+				            treedata->tree, ray_org_local, ray_normal_local,
+				            true, ob_scale, &nearest, test_vert_depth_cb, &userdata)) != -1)
 				{
-					const MVert *v = &treedata->vert[nearest.index];
-					float vno[3];
-					normal_short_to_float_v3(vno, v->no);
-					retval = snapVertex(
-					             ar, v->co, vno, obmat, timat, mval, dist_px,
-					             ray_start, ray_start_local, ray_normal_local, ray_depth,
-					             r_loc, r_no);
+					copy_v3_v3(r_loc, nearest.co);
+					mul_m4_v3(obmat, r_loc);
+					if (r_no) {
+						copy_v3_v3(r_no, nearest.no);
+						mul_m3_v3(timat, r_no);
+						normalize_v3(r_no);
+					}
+					*dist_px *= nearest.dist_sq / (*dist_to_ray_sq);
+					*dist_to_ray_sq = nearest.dist_sq;
+
+					retval = true;
 				}
 				break;
 			}
 			case SCE_SNAP_MODE_EDGE:
 			{
-				MVert *verts = dm->getVertArray(dm);
-				MEdge *edges = dm->getEdgeArray(dm);
-				int totedge = dm->getNumEdges(dm);
-
-				for (int i = 0; i < totedge; i++) {
-					MEdge *e = edges + i;
-					retval |= snapEdge(
-					        ar, verts[e->v1].co, verts[e->v1].no, verts[e->v2].co, verts[e->v2].no,
-					        obmat, timat, mval, dist_px,
-					        ray_start, ray_start_local, ray_normal_local, ray_depth,
-					        r_loc, r_no);
-				}
+				float ray_org_local[3];
+
+				copy_v3_v3(ray_org_local, ray_origin);
+				mul_m4_v3(imat, ray_org_local);
 
+				BVHTreeNearest nearest;
+
+				nearest.index = -1;
+				nearest.dist_sq = *dist_to_ray_sq;
+
+				struct NearestDM_Data userdata;
+				userdata.bvhdata = treedata;
+				userdata.is_persp = is_persp;
+				userdata.ray_depth_range = ray_depth_range;
+				userdata.ray_depth = ray_depth;
+
+				float ob_scale[3];
+				mat4_to_size(ob_scale, obmat);
+
+				if (treedata->tree && (
+				    is_persp ?
+				    BLI_bvhtree_find_nearest_to_ray_angle(
+				            treedata->tree, ray_org_local, ray_normal_local,
+				            true, ob_scale, &nearest, test_edge_depth_cb, &userdata) :
+				    BLI_bvhtree_find_nearest_to_ray(
+				            treedata->tree, ray_org_local, ray_normal_local,
+				            true, ob_scale, &nearest, test_edge_depth_cb, &userdata)) != -1)
+				{
+					copy_v3_v3(r_loc, nearest.co);
+					mul_m4_v3(obmat, r_loc);
+					if (r_no) {
+						copy_v3_v3(r_no, nearest.no);
+						mul_m3_v3(timat, r_no);
+						normalize_v3(r_no);
+					}
+					*dist_px *= nearest.dist_sq / (*dist_to_ray_sq);
+					*dist_to_ray_sq = nearest.dist_sq;
+
+					retval = true;
+				}
 				break;
 			}
 		}
@@ -955,17 +1181,51 @@ static bool snapDerivedMesh(
 	return retval;
 }
 
+static void test_bmvert_depth_cb(
+        void *userdata, const float origin[3], const float dir[3],
+        const float scale[3], int index, BVHTreeNearest *nearest)
+{
+	struct NearestDM_Data *ndata = userdata;
+	const BMEditMesh *em = ndata->bvhdata;
+	BMVert *eve = BM_vert_at_index(em->bm, index);
+
+	if (test_vert(
+	        eve->co, eve->no, origin, dir, ndata->ray_depth_range, scale, ndata->is_persp,
+	        ndata->ray_depth, &nearest->dist_sq,
+	        nearest->co, nearest->no))
+	{
+		nearest->index = index;
+	}
+}
+
+static void test_bmedge_depth_cb(
+        void *userdata, const float origin[3], const float dir[3],
+        const float scale[3], int index, BVHTreeNearest *nearest)
+{
+	struct NearestDM_Data *ndata = userdata;
+	const BMEditMesh *em = ndata->bvhdata;
+	BMEdge *eed = BM_edge_at_index(em->bm, index);
+
+	if (test_edge(
+	        eed->v1->co, eed->v2->co, origin, dir, ndata->ray_depth_range, scale, ndata->is_persp,
+	        ndata->ray_depth, &nearest->dist_sq,
+	        nearest->co, nearest->no))
+	{
+		nearest->index = index;
+	}
+}
 
 static bool snapEditMesh(
         SnapObjectContext *sctx,
-        Object *ob, BMEditMesh *em, float obmat[4][4],
-        const float mval[2], float *dist_px, const short snap_to,
-        const float ray_start[3], const float ray_normal[3], const float ray_origin[3],
-        float *ray_depth, const unsigned int ob_index,
+        Object *ob, BMEditMesh *em, float obmat[4][4], const unsigned int ob_index,
+        float *dist_px, const short snap_to,
+        const float ray_origin[3], const float ray_start[3], const float ray_normal[3], const float ray_depth_range[2],
+        /* read/write args */
+        float *ray_depth, float *dist_to_ray_sq,
+        /* return args */
         float r_loc[3], float r_no[3], int *r_index,
         ListBase *r_hit_list)
 {
-	const ARegion *ar = sctx->v3d_data.ar;
 	bool retval = false;
 
 	if (snap_to == SCE_SNAP_MODE_FACE) {
@@ -985,31 +1245,19 @@ static bool snapEditMesh(
 	}
 
 	{
-		const bool do_ray_start_correction = (
-		        ELEM(snap_to, SCE_SNAP_MODE_FACE, SCE_SNAP_MODE_VERTEX) &&
-		        (sctx->use_v3d && !((RegionView3D *)sctx->v3d_data.ar->regiondata)->is_persp));
+		const bool is_persp = (sctx->use_v3d && ((RegionView3D *)sctx->v3d_data.ar->regiondata)->is_persp);
 
 		float imat[4][4];
 		float timat[3][3]; /* transpose inverse matrix for normals */
-		float ray_start_local[3], ray_normal_local[3];
-		float local_scale, local_depth;
+		float ray_normal_local[3];
 
 		invert_m4_m4(imat, obmat);
 		transpose_m3_m4(timat, imat);
 
-		copy_v3_v3(ray_start_local, ray_start);
 		copy_v3_v3(ray_normal_local, ray_normal);
 
-		mul_m4_v3(imat, ray_start_local);
 		mul_mat3_m4_v3(imat, ray_normal_local);
 
-		/* local scale in normal direction */
-		local_scale = normalize_v3(ray_normal_local);
-		local_depth = *ray_depth;
-		if (local_depth != BVH_RAYCAST_DIST_MAX) {
-			local_depth *= local_scale;
-		}
-
 		SnapObjectData_EditMesh *sod = NULL;
 
 		BVHTreeFromEditMesh *treedata = NULL, treedata_stack;
@@ -1027,6 +1275,9 @@ static bool snapEditMesh(
 			int tree_index = -1;
 			switch (snap_to) {
 				case SCE_SNAP_MODE_FACE:
+					tree_index = 2;
+					break;
+				case SCE_SNAP_MODE_EDGE:
 					tree_index = 1;
 					break;
 				case SCE_SNAP_MODE_VERTEX:
@@ -1041,10 +1292,8 @@ static bool snapEditMesh(
 			}
 		}
 		else {
-			if (ELEM(snap_to, SCE_SNAP_MODE_FACE, SCE_SNAP_MODE_VERTEX)) {
-				treedata = &treedata_stack;
-				memset(treedata, 0, sizeof(*treedata));
-			}
+			treedata = &treedata_stack;
+			memset(treedata, 0, sizeof(*treedata));
 		}
 
 		if (treedata && treedata->tree == NULL) {
@@ -1065,6 +1314,23 @@ static bool snapEditMesh(
 					}
 					break;
 				}
+				case SCE_SNAP_MODE_EDGE:
+				{
+					BLI_bitmap *edges_mask = NULL;
+					int edges_num_active = -1;
+					if (sctx->callbacks.edit_mesh.test_edge_fn) {
+						edges_mask = BLI_BITMAP_NEW(em->bm->totedge, __func__);
+						edges_num_active = BM_iter_mesh_bitmap_from_filter(
+						        BM_EDGES_OF_MESH, em->bm, edges_mask,
+						        (bool (*)(BMElem *, void *))sctx->callbacks.edit_mesh.test_edge_fn,
+						        sctx->callbacks.edit_mesh.user_data);
+					}
+					bvhtree_from_editmesh_edges_ex(treedata, em, edges_mask, edges_num_active, 0.0f, 2, 6);
+					if (edges_mask) {
+						MEM_freeN(edges_mask);
+					}
+					break;
+				}
 				case SCE_SNAP_MODE_VERTEX:
 				{
 					BLI_bitmap *verts_mask = NULL;
@@ -1085,43 +1351,55 @@ static bool snapEditMesh(
 			}
 		}
 
-		/* Only use closer ray_start in case of ortho view! In perspective one, ray_start may already
-		 * been *inside* boundbox, leading to snap failures (see T38409).
-		 * Note also ar might be null (see T38435), in this case we assume ray_start is ok!
-		 */
-		float len_diff = 0.0f;
-		if (do_ray_start_correction) {
-			/* We *need* a reasonably valid len_diff in this case.
-			 * Use BHVTree to find the closest face from ray_start_local.
-			 */
-			if (treedata && treedata->tree != NULL) {
-				BVHTreeNearest nearest;
-				nearest.index = -1;
-				nearest.dist_sq = FLT_MAX;
-				/* Compute and store result. */
-				if (BLI_bvhtree_find_nearest(
-				        treedata->tree, ray_start_local, &nearest, treedata->nearest_callback, treedata) != -1)
-				{
-					len_diff = sqrtf(nearest.dist_sq);
-					float ray_org_local[3];
-
-					copy_v3_v3(ray_org_local, ray_origin);
-					mul_m4_v3(imat, ray_org_local);
-
-					/* We pass a temp ray_start, set from object's boundbox, to avoid precision issues with very far
-					* away ray_start values (as returned in case of ortho view3d), see T38358.
-					*/
-					len_diff -= local_scale;  /* make temp start point a bit away from bbox hit point. */
-					madd_v3_v3v3fl(ray_start_local, ray_org_local, ray_normal_local,
-					               len_diff - len_v3v3(ray_start_local, ray_org_local));
-					local_depth -= len_diff;
-				}
-			}
-		}
-
 		switch (snap_to) {
 			case SCE_SNAP_MODE_FACE:
 			{
+				float ray_start_local[3];
+				copy_v3_v3(ray_start_local, ray_start);
+				mul_m4_v3(imat, ray_start_local);
+
+				/* local scale in normal direction */
+				float local_scale = normalize_v3(ray_normal_local);
+				float local_depth = *ray_depth;
+				if (local_depth != BVH_RAYCAST_DIST_MAX) {
+					local_depth *= local_scale;
+				}
+
+				/* Only use closer ray_start in case of ortho view! In perspective one, ray_start may already
+				 * been *inside* boundbox, leading to snap failures (see T38409).
+				 * Note also ar might be null (see T38435), in this case we assume ray_start is ok!
+				 */
+				float len_diff = 0.0f;
+				if (sctx->use_v3d && !is_persp) { /* do_ray_start_correction */
+					/* We *need* a reasonably valid len_diff in this case.
+					 * Use BHVTree to find the closest face from ray_start_local.
+					 */
+					if (treedata && treedata->tree != NULL) {
+						BVHTreeNearest nearest;
+						nearest.index = -1;
+						nearest.dist_sq = FLT_MAX;
+						/* Compute and store result. */
+						if (BLI_bvhtree_find_nearest(
+						        treedata->tree, ray_start_local, &nearest, NULL, NULL) != -1)
+						{
+							float dvec[3];
+							sub_v3_v3v3(dvec, nearest.co, ray_start_local);
+							len_diff = dot_v3v3(dvec, ray_normal_local);
+							float ray_org_local[3];
+
+							copy_v3_v3(ray_org_local, ray_origin);
+							mul_m4_v3(imat, ray_org_local);
+
+							/* We pass a temp ray_start, set from object's boundbox, to avoid precision issues with very far
+							* away ray_start values (as returned in case of ortho view3d), see T38358.
+							*/
+							len_diff -= local_scale;  /* make temp start point a bit away from bbox hit point. */
+							madd_v3_v3v3fl(ray_start_local, ray_org_local, ray_normal_local,
+							               len_diff + ray_depth_range[0]);
+							local_depth -= len_diff;
+						}
+					}
+				}
 				if (r_hit_list) {
 					struct RayCastAll_Data data;
 
@@ -1176,47 +1454,92 @@ static bool snapEditMesh(
 				}
 				break;
 			}
-			case SCE_SNAP_MODE_VERTEX:
+			case SCE_SNAP_MODE_EDGE:
 			{
+				float ray_org_local[3];
+
+				copy_v3_v3(ray_org_local, ray_origin);
+				mul_m4_v3(imat, ray_org_local);
+
 				BVHTreeNearest nearest;
 
 				nearest.index = -1;
-				nearest.dist_sq = local_depth * local_depth;
-				if (treedata->tree &&
+				nearest.dist_sq = *dist_to_ray_sq;
+
+				struct NearestDM_Data userdata;
+				userdata.bvhdata = em;
+				userdata.is_persp = is_persp;
+				userdata.ray_depth_range = ray_depth_range;
+				userdata.ray_depth = ray_depth;
+
+				float ob_scale[3];
+				mat4_to_size(ob_scale, obmat);
+
+				if (treedata->tree && (
+				    is_persp ?
+				    BLI_bvhtree_find_nearest_to_ray_angle(
+				            treedata->tree, ray_org_local, ray_normal_local,
+				            false, ob_scale, &nearest, test_bmedge_depth_cb, &userdata) :
 				    BLI_bvhtree_find_nearest_to_ray(
-				            treedata->tree, ray_start_local, ray_normal_local,
-				            &nearest, NULL, NULL) != -1)
+				            treedata->tree, ray_org_local, ray_normal_local,
+				            false, ob_scale, &nearest, test_bmedge_depth_cb, &userdata)) != -1)
 				{
-					const BMVert *v = BM_vert_at_index(em->bm, nearest.index);
-					retval = snapVertex(
-					        ar, v->co, v->no, obmat, timat, mval, dist_px,
-					        ray_start, ray_start_local, ray_normal_local, ray_depth,
-					        r_loc, r_no);
+					copy_v3_v3(r_loc, nearest.co);
+					mul_m4_v3(obmat, r_loc);
+					if (r_no) {
+						copy_v3_v3(r_no, nearest.no);
+						mul_m3_v3(timat, r_no);
+						normalize_v3(r_no);
+					}
+					*dist_px *= nearest.dist_sq / (*dist_to_ray_sq);
+					*dist_to_ray_sq = nearest.dist_sq;
+
+					retval = true;
 				}
 				break;
 			}
-			case SCE_SNAP_MODE_EDGE:
+			case SCE_SNAP_MODE_VERTEX:
 			{
-				BM_mesh_elem_table_ensure(em->bm, BM_EDGE);
-				int totedge = em->bm->totedge;
-				for (int i = 0; i < totedge; i++) {
-					BMEdge *eed = BM_edge_at_index(em->bm, i);
-
-					if (!BM_elem_flag_test(eed, BM_ELEM_HIDDEN) &&
-					    !BM_elem_flag_test(eed->v1, BM_ELEM_SELECT) &&
-					    !BM_elem_flag_test(eed->v2, BM_ELEM_SELECT))
-					{
-						short v1no[3], v2no[3];
-						normal_float_to_short_v3(v1no, eed->v1->no);
-						normal_float_to_short_v3(v2no, eed->v2->no);
-						retval |= snapEdge(
-						        ar, eed->v1->co, v1no, eed->v2->co, v2no,
-						        obmat, timat, mval, dist_px,
-						        ray_start, ray_start_local, ray_normal_local, ray_depth,
-						        r_loc, r_no);
+				float ray_org_local[3];
+
+				copy_v3_v3(ray_org_local, ray_origin);
+				mul_m4_v3(imat, ray_org_local);
+
+				BVHTreeNearest nearest;
+
+				nearest.index = -1;
+				nearest.dist_sq = *dist_to_ray_sq;
+
+				struct NearestDM_Data userdata;
+				userdata.bvhdata = em;
+				userdata.is_persp = is_persp;
+				userdata.ray_depth_range = ray_depth_range;
+				userdata.ray_depth = ray_depth;
+
+				float ob_scale[3];
+				mat4_to_size(ob_scale, obmat);
+
+				if (treedata->tree && (
+				    is_persp ?
+				    BLI_bvhtree_find_nearest_to_ray_angle(
+				            treedata->tree, ray_org_local, ray_normal_local,
+				            false, ob_scale, &nearest, test_bmvert_depth_cb, &userdata) :
+				    BLI_bvhtree_find_nearest_to_ray(
+				            treedata->tree, ray_org_local, ray_normal_local,
+				            false, ob_scale, &nearest, test_bmedge_depth_cb, &userdata)) != -1)
+				{
+					copy_v3_v3(r_loc, nearest.co);
+					mul_m4_v3(obmat, r_loc);
+					if (r_no) {
+						copy_v3_v3(r_no, nearest.no);
+						mul_m3_v3(timat, r_no);
+						normalize_v3(r_no);
 					}
-				}
+					*dist_px *= nearest.dist_sq / (*dist_to_ray_sq);
+					*dist_to_ray_sq = nearest.dist_sq;
 
+					retval = true;
+				}
 				break;
 			}
 		}
@@ -1231,12 +1554,22 @@ static bool snapEditMesh(
 	return retval;
 }
 
+/**
+ * \param use_obedit: Uses the coordinates of BMesh (if any) to do the snapping;
+ * \param ray_depth_range:
+ * - 0: distance from the ray_origin to the clipping plane min (can be negative).
+ * - 1: maximum distance, elements outside this are ignored.
+ * \param ray_depth: maximum depth allowed for r_co.
+ *
+ * \note Duplicate args here are documented at #snapObjectsRay
+ */
 static bool snapObject(
         SnapObjectContext *sctx,
-        Object *ob, float obmat[4][4], bool use_obedit, const short snap_to,
-        const float mval[2], float *dist_px, const unsigned int ob_index,
-        const float ray_start[3], const float ray_normal[3], const float ray_origin[3],
-        float *ray_depth,
+        Object *ob, float obmat[4][4], const unsigned int ob_index,
+        bool use_obedit, const short snap_to, const float mval[2],
+        const float ray_origin[3], const float ray_start[3], const float ray_normal[3], const float ray_depth_range[2],
+        /* read/write args */
+        float *ray_depth, float *dist_to_ray_sq, float *dist_px,
         /* return args */
         float r_loc[3], float r_no[3], int *r_index,
         Object **r_ob, float r_obmat[4][4],
@@ -1251,9 +1584,10 @@ static bool snapObject(
 		if (use_obedit) {
 			em = BKE_editmesh_from_object(ob);
 			retval = snapEditMesh(
-			        sctx, ob, em, obmat, mval, dist_px, snap_to,
-			        ray_start, ray_normal, ray_origin,
-			        ray_depth, ob_index,
+			        sctx, ob, em, obmat, ob_index,
+			        dist_px, snap_to,
+			        ray_origin, ray_start, ray_normal, ray_depth_range,
+			        ray_depth, dist_to_ray_sq,
 			        r_loc, r_no, r_index,
 			        r_hit_list);
 		}
@@ -1269,36 +1603,42 @@ static bool snapObject(
 				dm = mesh_get_derived_final(sctx->scene, ob, CD_MASK_BAREMESH);
 			}
 			retval = snapDerivedMesh(
-			        sctx, ob, dm, obmat, mval, dist_px, snap_to, true,
-			        ray_start, ray_normal, ray_origin,
-			        ray_depth, ob_index,
-			        r_loc, r_no, r_index, r_hit_list);
+			        sctx, ob, dm, obmat, ob_index,
+			        snap_to, true,
+			        ray_origin, ray_start, ray_normal, ray_depth_range,
+			        ray_depth, dist_to_ray_sq, dist_px,
+			        r_loc, r_no,
+			        r_index, r_hit_list);
 
 			dm->release(dm);
 		}
 	}
 	else if (ob->type == OB_ARMATURE) {
 		retval = snapArmature(
-		        ar, ob, ob->data, obmat, mval, dist_px, snap_to,
-		        ray_start, ray_normal, ray_depth,
+		        ar, ob, ob->data, obmat, mval, snap_to,
+		        ray_start, ray_normal,
+		        ray_depth, dist_px,
 		        r_loc, r_no);
 	}
 	else if (ob->type == OB_CURVE) {
 		retval = snapCurve(
-		        ar, ob, ob->data, obmat, mval, dist_px, snap_to,
-		        ray_start, ray_normal, ray_depth,
+		        ar, ob, ob->data, obmat, mval, snap_to,
+		        ray_start, ray_normal,
+		        ray_depth, dist_px,
 		        r_loc, r_no);
 	}
 	else if (ob->type == OB_EMPTY) {
 		retval = snapEmpty(
-		        ar, ob, obmat, mval, dist_px, snap_to,
-		        ray_start, ray_normal, ray_depth,
+		        ar, ob, obmat, mval, snap_to,
+		        ray_start, ray_normal,
+		        ray_depth, dist_px,
 		        r_loc, r_no);
 	}
 	else if (ob->type == OB_CAMERA) {
 		retval = snapCamera(
-		        ar, sctx->scene, ob, obmat, mval, dist_px, snap_to,
-		        ray_start, ray_normal, ray_depth,
+		        ar, sctx->scene, ob, obmat, mval, snap_to,
+		        ray_start, ray_normal,
+		        ray_depth, dist_px,
 		        r_loc, r_no);
 	}
 
@@ -1312,18 +1652,68 @@ static bool snapObject(
 	return retval;
 }
 
+/**
+ * Main Snapping Function
+ * ======================
+ *
+ * Walks through all objects in the scene to find the closest snap element ray.
+ *
+ * \param sctx: Snap context to store data.
+ * \param snap_to: Element to snap, Vertice, Edge or Face.
+ * Currently only works one at a time, but can eventually operate as flag.
+ *
+ * \param snap_select: from enum SnapSelect.
+ *
+ * \param use_object_edit_cage: Uses the coordinates of BMesh (if any) to do the snapping.
+ * \param mval: Optional screen-space 2D location we're snapping to (may phase out).
+ * \param ray_origin: ray_start before being moved toward the ray_normal at the distance from vew3d clip_min.
+ * \param ray_start: ray_origin moved for the start clipping plane (clip_min).
+ * \param ray_normal: Unit length direction of the ray.
+ *
+ * Read/Write Args
+ * ---------------
+ *
+ * \param ray_depth: maximum depth allowed for r_co, elements deeper than this value will be ignored.
+ * \param dist_to_ray_sq: Real distance (3D) or Tangent (view cone radius at distance 1.0) squared.
+ * resulting of the function #dist_px_to_dist3d_or_tangent.
+ *
+ * \param dist_px: Pixel distance to element,
+ * note that this will eventually be replaced entirely by \a dist_to_ray_sq.
+ *
+ * Output Args
+ * -----------
+ *
+ * \param r_loc: Hit location.
+ * \param r_no: Hit normal (optional).
+ * \param r_index: Hit index or -1 when no valid index is found.
+ * (currently only set to the polygon index when when using ``snap_to == SCE_SNAP_MODE_FACE``).
+ * \param r_ob: Hit object.
+ * \param r_obmat: Object matrix (may not be #Object.obmat with dupli-instances).
+ * \param r_hit_list: List of #SnapObjectHitDepth (caller must free).
+ *
+ */
 static bool snapObjectsRay(
         SnapObjectContext *sctx,
         const unsigned short snap_to, const SnapSelect snap_select,
-        const bool use_object_edit_cage,
-        const float mval[2], float *dist_px,
-        const float ray_start[3], const float ray_normal[3], const float ray_origin[3], float *ray_depth,
+        const bool use_object_edit_cage, const float mval[2],
+        const float ray_origin[3], const float ray_start[3], const float ray_normal[3],
+        /* read/write args */
+        float *ray_depth, float *dist_to_ray_sq, float *dist_px,
         /* return args */
         float r_loc[3], float r_no[3], int *r_index,
         Object **r_ob, float r_obmat[4][4],
         ListBase *r_hit_list)
 {
 	bool retval = false;
+
+	float dvec[3];
+	sub_v3_v3v3(dvec, ray_start, ray_origin);
+
+	const float ray_depth_range[2] = {
+		dot_v3v3(dvec, ray_normal),
+		*ray_depth,
+	};
+
 	unsigned int ob_index = 0;
 	Object *obedit = use_object_edit_cage ? sctx->scene->obedit : NULL;
 
@@ -1337,9 +1727,10 @@ static bool snapObjectsRay(
 		Object *ob = base_act->object;
 
 		retval |= snapObject(
-		        sctx, ob, ob->obmat, false, snap_to,
-		        mval, dist_px, ob_index++,
-		        ray_start, ray_normal, ray_origin, ray_depth,
+		        sctx, ob, ob->obmat, ob_index++,
+		        false, snap_to, mval,
+		        ray_origin, ray_start, ray_normal, ray_depth_range,
+		        ray_depth, dist_to_ray_sq, dist_px,
 		        r_loc, r_no, r_index, r_ob, r_obmat, r_hit_list);
 	}
 
@@ -1372,9 +1763,10 @@ static bool snapObjectsRay(
 					Object *dupli_snap = (use_obedit_dupli) ? obedit : dupli_ob->ob;
 
 					retval |= snapObject(
-					        sctx, dupli_snap, dupli_ob->mat, use_obedit_dupli, snap_to,
-					        mval, dist_px, ob_index++,
-					        ray_start, ray_normal, ray_origin, ray_depth,
+					        sctx, dupli_snap, dupli_ob->mat, ob_index++,
+					        use_obedit_dupli, snap_to, mval,
+					        ray_origin, ray_start, ray_normal, ray_depth_range,
+					        ray_depth, dist_to_ray_sq, dist_px,
 					        r_loc, r_no, r_index, r_ob, r_obmat, r_hit_list);
 				}
 
@@ -1385,9 +1777,10 @@ static bool snapObjectsRay(
 			Object *ob_snap = use_obedit ? obedit : ob;
 
 			retval |= snapObject(
-			        sctx, ob_snap, ob->obmat, use_obedit, snap_to,
-			        mval, dist_px, ob_index++,
-			        ray_start, ray_normal, ray_origin, ray_depth,
+			        sctx, ob_snap, ob->obmat, ob_index++,
+			        use_obedit, snap_to, mval,
+			        ray_origin, ray_start, ray_normal, ray_depth_range,
+			        ray_depth, dist_to_ray_sq, dist_px,
 			        r_loc, r_no, r_index, r_ob, r_obmat, r_hit_list);
 		}
 	}
@@ -1489,15 +1882,19 @@ bool ED_transform_snap_object_project_ray_ex(
         SnapObjectContext *sctx,
         const unsigned short snap_to,
         const struct SnapObjectParams *params,
-        const float ray_start[3], const float ray_normal[3], float *ray_depth,
+        const float ray_start[3], const float ray_normal[3],
+        float *ray_depth,
         float r_loc[3], float r_no[3], int *r_index,
         Object **r_ob, float r_obmat[4][4])
 {
+	float dist_to_ray_sq = 0.0f;
+
 	return snapObjectsRay(
 	        sctx,
 	        snap_to, params->snap_select, params->use_object_edit_cage,
-	        NULL, NULL,
-	        ray_start, ray_normal, ray_start, ray_depth,
+	        NULL,
+	        ray_start, ray_start, ray_normal,
+	        ray_depth, &dist_to_ray_sq, NULL,
 	        r_loc, r_no, r_index, r_ob, r_obmat, NULL);
 }
 
@@ -1516,6 +1913,8 @@ bool ED_transform_snap_object_project_ray_all(
         float ray_depth, bool sort,
         ListBase *r_hit_list)
 {
+	float dist_to_ray_sq = 0.0f;
+
 	if (ray_depth == -1.0f) {
 		ray_depth = BVH_RAYCAST_DIST_MAX;
 	}
@@ -1526,9 +1925,9 @@ bool ED_transform_snap_object_project_ray_all(
 
 	bool retval = snapObjectsRay(
 	        sctx,
-	        snap_to, params->snap_select, params->use_object_edit_cage,
-	        NULL, NULL,
-	        ray_start, ray_normal, ray_start, &ray_depth,
+	        snap_to, params->snap_select, params->use_object_edit_cage, NULL,
+	        ray_start, ray_start, ray_normal,
+	        &ray_depth, &dist_to_ray_sq, NULL,
 	        NULL, NULL, NULL, NULL, NULL,
 	        r_hit_list);
 
@@ -1637,6 +2036,21 @@ static bool transform_snap_context_project_view3d_mixed_impl(
 }
 
 /**
+ * From a threshold (maximum distance to snap in pixels) returns:
+ *
+ * - The *real* distance (3D) if you are in orthographic-view.
+ * - The *tangent* (view cone radius at distance 1.0) if you are in perspective-view.
+ */
+static float dist_px_to_dist3d_or_tangent(const ARegion *ar, const float dist_px)
+{
+	const RegionView3D *rv3d = ar->regiondata;
+	if (ar->winx >= ar->winy)
+		return 2 * (dist_px / ar->winx) / rv3d->winmat[0][0];
+	else
+		return 2 * (dist_px / ar->winy) / rv3d->winmat[1][1];
+}
+
+/**
  * Convenience function for performing snapping.
  *
  * Given a 2D region value, snap to vert/edge/face.
@@ -1687,11 +2101,32 @@ bool ED_transform_snap_object_project_view3d_ex(
 		return false;
 	}
 
+	float dist_to_ray_sq;
+	{
+		float radius = dist_px_to_dist3d_or_tangent(sctx->v3d_data.ar, *dist_px);
+		/**
+		 * Workaround to use of cone (Instead of project the radius on view plane):
+		 * In perspective view, the radius of the cone may decrease depending on the ray direction.
+		 * This is more evident with small values of the `Viewport lens angle`.
+		 * The threshold becomes distorted that way.
+		 */
+		RegionView3D *rv3d = sctx->v3d_data.ar->regiondata;
+		if (rv3d->is_persp) {
+			float view_dir[3];
+			negate_v3_v3(view_dir, rv3d->viewinv[2]);
+			normalize_v3(view_dir);
+			radius *= dot_v3v3(ray_normal, view_dir);
+		}
+
+		dist_to_ray_sq = SQUARE(radius);
+	}
+
 	return snapObjectsRay(
 	        sctx,
 	        snap_to, params->snap_select, params->use_object_edit_cage,
-	        mval, dist_px,
-	        ray_start, ray_normal, ray_origin, ray_depth,
+	        mval,
+	        ray_origin, ray_start, ray_normal,
+	        ray_depth, &dist_to_ray_sq, dist_px,
 	        r_loc, r_no, r_index, NULL, NULL, NULL);
 }