BKE: Add 'mesh remap' code.

This is the (big!) core of mesh transfer data, it defines a set of structures
to represent a mapping of mesh elements (verts, edges, polys of loops) between
two arbitrary meshes, and code to compute such mappings.

No similarity is required between source and destination meshes (though results
when using complete different meshes are rather unlikely to be useful!).

This code is not bound to data transfer, it is defined to be as generic as possible,
and easy to reuse or extend as needs arise.

Several methods of mapping generation are defined for each element type,
we probably will have to adjust that in future (remove useless ones, add
new ones...).

For loops, you can also define islands (for UVs e.g.) so that loops of a same
destination polygon do not 'spread' across several source islands.

Heavily reviewed and enhanced by Campbell, thanks a lot!

5 files changed, 2554 insertions, 45 deletions

diff --git a/source/blender/blenkernel/BKE_mesh_mapping.h b/source/blender/blenkernel/BKE_mesh_mapping.h
index ed7e506941c..da44c989146 100644
--- a/source/blender/blenkernel/BKE_mesh_mapping.h
+++ b/source/blender/blenkernel/BKE_mesh_mapping.h
@@ -31,8 +31,9 @@
  *  \ingroup bke
  */
 
-struct MPoly;
+struct MVert;
 struct MEdge;
+struct MPoly;
 struct MLoop;
 struct MLoopUV;
 
@@ -109,6 +110,10 @@ void BKE_mesh_vert_poly_map_create(
         MeshElemMap **r_map, int **r_mem,
         const struct MPoly *mface, const struct MLoop *mloop,
         int totvert, int totface, int totloop);
+void BKE_mesh_vert_loop_map_create(
+        MeshElemMap **r_map, int **r_mem,
+        const struct MPoly *mface, const struct MLoop *mloop,
+        int totvert, int totface, int totloop);
 void BKE_mesh_vert_edge_map_create(
         MeshElemMap **r_map, int **r_mem,
         const struct MEdge *medge, int totvert, int totedge);
@@ -122,7 +127,61 @@ void BKE_mesh_origindex_map_create(
         const int totorig,
         const int *final_origindex, const int totfinal);
 
-/* smoothgroups */
+
+/* islands */
+
+/* Loop islands data helpers. */
+enum {
+	MISLAND_TYPE_NONE = 0,
+	MISLAND_TYPE_VERT = 1,
+	MISLAND_TYPE_EDGE = 2,
+	MISLAND_TYPE_POLY = 3,
+	MISLAND_TYPE_LOOP = 4,
+};
+
+typedef struct MeshIslandStore {
+	short item_type;      /* MISLAND_TYPE_... */
+	short island_type;    /* MISLAND_TYPE_... */
+	short innercut_type;  /* MISLAND_TYPE_... */
+
+	int  items_to_islands_num;
+	int *items_to_islands;  /* map the item to the island index */
+
+	int                  islands_num;
+	size_t               islands_num_alloc;
+	struct MeshElemMap **islands;    /* Array of pointers, one item per island. */
+	struct MeshElemMap **innercuts;  /* Array of pointers, one item per island. */
+
+	struct MemArena *mem;  /* Memory arena, internal use only. */
+} MeshIslandStore;
+
+void BKE_mesh_loop_islands_init(
+        MeshIslandStore *island_store,
+        const short item_type, const int item_num, const short island_type, const short innercut_type);
+void BKE_mesh_loop_islands_clear(MeshIslandStore *island_store);
+void BKE_mesh_loop_islands_free(MeshIslandStore *island_store);
+void BKE_mesh_loop_islands_add(
+        MeshIslandStore *islands, const int item_num, int *item_indices,
+        const int num_island_items, int *island_item_indices,
+        const int num_innercut_items, int *innercut_item_indices);
+
+typedef bool (*MeshRemapIslandsCalc)(
+        struct MVert *verts, const int totvert,
+        struct MEdge *edges, const int totedge,
+        struct MPoly *polys, const int totpoly,
+        struct MLoop *loops, const int totloop,
+        struct MeshIslandStore *r_island_store);
+
+/* Above vert/UV mapping stuff does not do what we need here, but does things we do not need here.
+ * So better keep them separated for now, I think.
+ */
+bool BKE_mesh_calc_islands_loop_poly_uv(
+        struct MVert *verts, const int totvert,
+        struct MEdge *edges, const int totedge,
+        struct MPoly *polys, const int totpoly,
+        struct MLoop *loops, const int totloop,
+        MeshIslandStore *r_island_store);
+
 int *BKE_mesh_calc_smoothgroups(
         const struct MEdge *medge, const int totedge,
         const struct MPoly *mpoly, const int totpoly,
diff --git a/source/blender/blenkernel/BKE_mesh_remap.h b/source/blender/blenkernel/BKE_mesh_remap.h
new file mode 100644
index 00000000000..e1f37a63ff0
--- /dev/null
+++ b/source/blender/blenkernel/BKE_mesh_remap.h
@@ -0,0 +1,171 @@
+/*
+ * ***** BEGIN GPL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ */
+
+#ifndef __BKE_MESH_REMAP_H__
+#define __BKE_MESH_REMAP_H__
+
+/** \file BKE_mesh_remap.h
+ *  \ingroup bke
+ */
+
+struct CustomData;
+struct DerivedMesh;
+struct MVert;
+struct MeshElemMap;
+struct MemArena;
+
+/* Generic ways to map some geometry elements from a source mesh to a dest one. */
+
+typedef struct MeshPairRemapItem {
+	int    sources_num;
+	int   *indices_src;     /* NULL if no source found. */
+	float *weights_src;     /* NULL if no source found, else, always normalized! */
+	/* UNUSED (at the moment)*/
+	// float  hit_dist;     /* FLT_MAX if irrelevant or no source found. */
+	int    island;          /* For loops only. */
+} MeshPairRemapItem;
+
+/* All mapping computing func return this. */
+typedef struct MeshPairRemap {
+	int                items_num;
+	MeshPairRemapItem *items;  /* array, one item per dest element. */
+
+	struct MemArena *mem;  /* memory arena, internal use only. */
+} MeshPairRemap;
+
+/* Helpers! */
+void BKE_mesh_remap_init(MeshPairRemap *map, const int items_num);
+void BKE_mesh_remap_free(MeshPairRemap *map);
+
+void BKE_mesh_remap_item_define_invalid(MeshPairRemap *map, const int index);
+
+/* TODO:
+ * Add other 'from/to' mapping sources, like e.g. using an UVMap, etc.
+ *     http://blenderartists.org/forum/showthread.php?346458-Move-Vertices-to-the-location-of-the-Reference-Mesh-based-on-the-UV-Position
+ * We could also use similar topology mappings inside a same mesh
+ * (cf. Campbell's 'select face islands from similar topology' wip work).
+ * Also, users will have to check, whether we can get rid of some modes here, not sure all will be useful!
+ */
+enum {
+	MREMAP_USE_VERT                      = 1 << 4,
+	MREMAP_USE_EDGE                      = 1 << 5,
+	MREMAP_USE_LOOP                      = 1 << 6,
+	MREMAP_USE_POLY                      = 1 << 7,
+
+	MREMAP_USE_NEAREST                   = 1 << 8,
+	MREMAP_USE_NORPROJ                   = 1 << 9,
+	MREMAP_USE_INTERP                    = 1 << 10,
+	MREMAP_USE_NORMAL                    = 1 << 11,
+
+	/* ***** Target's vertices ***** */
+	MREMAP_MODE_VERT                     = 1 << 24,
+	/* Nearest source vert. */
+	MREMAP_MODE_VERT_NEAREST             = MREMAP_MODE_VERT | MREMAP_USE_VERT | MREMAP_USE_NEAREST,
+
+	/* Nearest vertex of nearest edge. */
+	MREMAP_MODE_VERT_EDGE_NEAREST        = MREMAP_MODE_VERT | MREMAP_USE_EDGE | MREMAP_USE_NEAREST,
+	/* This one uses two verts of selected edge (weighted interpolation). */
+	/* Nearest point on nearest edge. */
+	MREMAP_MODE_VERT_EDGEINTERP_NEAREST  = MREMAP_MODE_VERT | MREMAP_USE_EDGE | MREMAP_USE_NEAREST | MREMAP_USE_INTERP,
+
+	/* Nearest vertex of nearest poly. */
+	MREMAP_MODE_VERT_POLY_NEAREST        = MREMAP_MODE_VERT | MREMAP_USE_POLY | MREMAP_USE_NEAREST,
+	/* Those two use all verts of selected poly (weighted interpolation). */
+	/* Nearest point on nearest poly. */
+	MREMAP_MODE_VERT_POLYINTERP_NEAREST  = MREMAP_MODE_VERT | MREMAP_USE_POLY | MREMAP_USE_NEAREST | MREMAP_USE_INTERP,
+	/* Point on nearest face hit by ray from target vertex's normal. */
+	MREMAP_MODE_VERT_POLYINTERP_VNORPROJ = MREMAP_MODE_VERT | MREMAP_USE_POLY | MREMAP_USE_NORPROJ | MREMAP_USE_INTERP,
+
+	/* ***** Target's edges ***** */
+	MREMAP_MODE_EDGE                     = 1 << 25,
+
+	/* Source edge which both vertices are nearest of dest ones. */
+	MREMAP_MODE_EDGE_VERT_NEAREST        = MREMAP_MODE_EDGE | MREMAP_USE_VERT | MREMAP_USE_NEAREST,
+
+	/* Nearest source edge (using mid-point). */
+	MREMAP_MODE_EDGE_NEAREST             = MREMAP_MODE_EDGE | MREMAP_USE_EDGE | MREMAP_USE_NEAREST,
+
+	/* Nearest edge of nearest poly (using mid-point). */
+	MREMAP_MODE_EDGE_POLY_NEAREST        = MREMAP_MODE_EDGE | MREMAP_USE_POLY | MREMAP_USE_NEAREST,
+
+	/* Cast a set of rays from along dest edge, interpolating its vertices' normals, and use hit source edges. */
+	MREMAP_MODE_EDGE_EDGEINTERP_VNORPROJ = MREMAP_MODE_EDGE | MREMAP_USE_VERT | MREMAP_USE_NORPROJ | MREMAP_USE_INTERP,
+
+	/* ***** Target's loops ***** */
+	/* Note: when islands are given to loop mapping func, all loops from the same destination face will always be mapped
+	 *       to loops of source faces within a same island, regardless of mapping mode. */
+	MREMAP_MODE_LOOP                     = 1 << 26,
+
+	/* Best normal-matching loop from nearest vert. */
+	MREMAP_MODE_LOOP_NEAREST_LOOPNOR     = MREMAP_MODE_LOOP | MREMAP_USE_LOOP | MREMAP_USE_VERT | MREMAP_USE_NEAREST | MREMAP_USE_NORMAL,
+	/* Loop from best normal-matching poly from nearest vert. */
+	MREMAP_MODE_LOOP_NEAREST_POLYNOR     = MREMAP_MODE_LOOP | MREMAP_USE_POLY | MREMAP_USE_VERT | MREMAP_USE_NEAREST | MREMAP_USE_NORMAL,
+
+	/* Loop from nearest vertex of nearest poly. */
+	MREMAP_MODE_LOOP_POLY_NEAREST        = MREMAP_MODE_LOOP | MREMAP_USE_POLY | MREMAP_USE_NEAREST,
+	/* Those two use all verts of selected poly (weighted interpolation). */
+	/* Nearest point on nearest poly. */
+	MREMAP_MODE_LOOP_POLYINTERP_NEAREST  = MREMAP_MODE_LOOP | MREMAP_USE_POLY | MREMAP_USE_NEAREST | MREMAP_USE_INTERP,
+	/* Point on nearest face hit by ray from target loop's normal. */
+	MREMAP_MODE_LOOP_POLYINTERP_LNORPROJ = MREMAP_MODE_LOOP | MREMAP_USE_POLY | MREMAP_USE_NORPROJ | MREMAP_USE_INTERP,
+
+	/* ***** Target's polygons ***** */
+	MREMAP_MODE_POLY                     = 1 << 27,
+
+	/* Nearest source poly. */
+	MREMAP_MODE_POLY_NEAREST             = MREMAP_MODE_POLY | MREMAP_USE_POLY | MREMAP_USE_NEAREST,
+	/* Source poly from best normal-matching dest poly. */
+	MREMAP_MODE_POLY_NOR                 = MREMAP_MODE_POLY | MREMAP_USE_POLY | MREMAP_USE_NORMAL,
+
+	/* Project dest poly onto source mesh using its normal, and use interpolation of all intersecting source polys. */
+	MREMAP_MODE_POLY_POLYINTERP_PNORPROJ = MREMAP_MODE_POLY | MREMAP_USE_POLY | MREMAP_USE_NORPROJ | MREMAP_USE_INTERP,
+
+	/* ***** Same topology, applies to all four elements types. ***** */
+	MREMAP_MODE_TOPOLOGY                 = MREMAP_MODE_VERT | MREMAP_MODE_EDGE | MREMAP_MODE_LOOP | MREMAP_MODE_POLY,
+};
+
+/* TODO add mesh2mesh versions (we'll need mesh versions of bvhtree funcs too, though!). */
+
+void BKE_mesh_remap_calc_verts_from_dm(
+        const int mode, const struct SpaceTransform *space_transform, const float max_dist, const float ray_radius,
+        const struct MVert *verts_dst, const int numverts_dst, const bool dirty_nors_dst,
+        struct DerivedMesh *dm_src, MeshPairRemap *r_map);
+
+void BKE_mesh_remap_calc_edges_from_dm(
+        const int mode, const struct SpaceTransform *space_transform, const float max_dist, const float ray_radius,
+        const struct MVert *verts_dst, const int numverts_dst, const struct MEdge *edges_dst, const int numedges_dst,
+        const bool dirty_nors_dst, struct DerivedMesh *dm_src, MeshPairRemap *r_map);
+
+void BKE_mesh_remap_calc_loops_from_dm(
+        const int mode, const struct SpaceTransform *space_transform, const float max_dist, const float ray_radius,
+        struct MVert *verts_dst, const int numverts_dst, struct MEdge *edges_dst, const int numedges_dst,
+        struct MLoop *loops_dst, const int numloops_dst, struct MPoly *polys_dst, const int numpolys_dst,
+        struct CustomData *ldata_dst, struct CustomData *pdata_dst,
+        const float split_angle_dst, const bool dirty_nors_dst,
+        struct DerivedMesh *dm_src,
+        MeshRemapIslandsCalc gen_islands_src, const float islands_precision_src, struct MeshPairRemap *r_map);
+
+void BKE_mesh_remap_calc_polys_from_dm(
+        const int mode, const struct SpaceTransform *space_transform, const float max_dist, const float ray_radius,
+        struct MVert *verts_dst, const int numverts_dst, struct MLoop *loops_dst, const int numloops_dst,
+        struct MPoly *polys_dst, const int numpolys_dst, struct CustomData *pdata_dst, const bool dirty_nors_dst,
+        struct DerivedMesh *dm_src, struct MeshPairRemap *r_map);
+
+#endif  /* __BKE_MESH_REMAP_H__ */
diff --git a/source/blender/blenkernel/CMakeLists.txt b/source/blender/blenkernel/CMakeLists.txt
index 5850b1614e7..5976e041703 100644
--- a/source/blender/blenkernel/CMakeLists.txt
+++ b/source/blender/blenkernel/CMakeLists.txt
@@ -121,6 +121,7 @@ set(SRC
 	intern/mesh.c
 	intern/mesh_evaluate.c
 	intern/mesh_mapping.c
+	intern/mesh_remap.c
 	intern/mesh_validate.c
 	intern/modifier.c
 	intern/modifiers_bmesh.c
@@ -230,6 +231,7 @@ set(SRC
 	BKE_mball.h
 	BKE_mesh.h
 	BKE_mesh_mapping.h
+	BKE_mesh_remap.h
 	BKE_modifier.h
 	BKE_movieclip.h
 	BKE_multires.h
diff --git a/source/blender/blenkernel/intern/mesh_mapping.c b/source/blender/blenkernel/intern/mesh_mapping.c
index 53d1aae104c..9e490ae6766 100644
--- a/source/blender/blenkernel/intern/mesh_mapping.c
+++ b/source/blender/blenkernel/intern/mesh_mapping.c
@@ -32,10 +32,12 @@
 #include "DNA_meshdata_types.h"
 
 #include "BLI_utildefines.h"
+#include "BLI_bitmap.h"
 #include "BLI_math.h"
 
 #include "BKE_mesh_mapping.h"
 #include "BKE_customdata.h"
+#include "BLI_memarena.h"
 
 #include "BLI_strict_flags.h"
 
@@ -156,18 +158,25 @@ void BKE_mesh_uv_vert_map_free(UvVertMap *vmap)
 	}
 }
 
-/* Generates a map where the key is the vertex and the value is a list
- * of polys that use that vertex as a corner. The lists are allocated
- * from one memory pool. */
-void BKE_mesh_vert_poly_map_create(MeshElemMap **r_map, int **r_mem,
-                                   const MPoly *mpoly, const MLoop *mloop,
-                                   int totvert, int totpoly, int totloop)
+/**
+
+
+ * Generates a map where the key is the vertex and the value is a list
+ * of polys or loops that use that vertex as a corner. The lists are allocated
+ * from one memory pool.
+ *
+ * Wrapped by #BKE_mesh_vert_poly_map_create & BKE_mesh_vert_loop_map_create
+ */
+static void mesh_vert_poly_or_loop_map_create(
+        MeshElemMap **r_map, int **r_mem,
+        const MPoly *mpoly, const MLoop *mloop,
+        int totvert, int totpoly, int totloop, const bool do_loops)
 {
-	MeshElemMap *map = MEM_callocN(sizeof(MeshElemMap) * (size_t)totvert, "vert poly map");
+	MeshElemMap *map = MEM_callocN(sizeof(MeshElemMap) * (size_t)totvert, __func__);
 	int *indices, *index_iter;
 	int i, j;
 
-	indices = index_iter = MEM_mallocN(sizeof(int) * (size_t)totloop, "vert poly map mem");
+	indices = index_iter = MEM_mallocN(sizeof(int) * (size_t)totloop, __func__);
 
 	/* Count number of polys for each vertex */
 	for (i = 0; i < totpoly; i++) {
@@ -193,7 +202,7 @@ void BKE_mesh_vert_poly_map_create(MeshElemMap **r_map, int **r_mem,
 		for (j = 0; j < p->totloop; j++) {
 			unsigned int v = mloop[p->loopstart + j].v;
 
-			map[v].indices[map[v].count] = i;
+			map[v].indices[map[v].count] = do_loops ? p->loopstart + j : i;
 			map[v].count++;
 		}
 	}
@@ -202,6 +211,28 @@ void BKE_mesh_vert_poly_map_create(MeshElemMap **r_map, int **r_mem,
 	*r_mem = indices;
 }
 
+/**
+ * Generates a map where the key is the vertex and the value is a list of polys that use that vertex as a corner.
+ * The lists are allocated from one memory pool.
+ */
+void BKE_mesh_vert_poly_map_create(MeshElemMap **r_map, int **r_mem,
+                                   const MPoly *mpoly, const MLoop *mloop,
+                                   int totvert, int totpoly, int totloop)
+{
+	mesh_vert_poly_or_loop_map_create(r_map, r_mem, mpoly, mloop, totvert, totpoly, totloop, false);
+}
+
+/**
+ * Generates a map where the key is the vertex and the value is a list of loops that use that vertex as a corner.
+ * The lists are allocated from one memory pool.
+ */
+void BKE_mesh_vert_loop_map_create(MeshElemMap **r_map, int **r_mem,
+                                   const MPoly *mpoly, const MLoop *mloop,
+                                   int totvert, int totpoly, int totloop)
+{
+	mesh_vert_poly_or_loop_map_create(r_map, r_mem, mpoly, mloop, totvert, totpoly, totloop, true);
+}
+
 /* Generates a map where the key is the vertex and the value is a list
  * of edges that use that vertex as an endpoint. The lists are allocated
  * from one memory pool. */
@@ -345,26 +376,30 @@ void BKE_mesh_origindex_map_create(MeshElemMap **r_map, int **r_mem,
 /** \} */
 
 
-
 /* -------------------------------------------------------------------- */
 
-/** \name Mesh Smooth Groups
+/** \name Mesh loops/poly islands.
+ * Used currently for UVs and 'smooth groups'.
  * \{ */
 
-/**
- * Calculate smooth groups from sharp edges.
- *
- * \param r_totgroup The total number of groups, 1 or more.
- * \return Polygon aligned array of group index values (bitflags if use_bitflags is true), starting at 1.
+/** Callback deciding whether the given poly/loop/edge define an island boundary or not.
  */
-int *BKE_mesh_calc_smoothgroups(const MEdge *medge, const int totedge,
-                                const MPoly *mpoly, const int totpoly,
-                                const MLoop *mloop, const int totloop,
-                                int *r_totgroup, const bool use_bitflags)
+typedef bool (*MeshRemap_CheckIslandBoundary)(
+        const struct MPoly *mpoly, const struct MLoop *mloop, const struct MEdge *medge,
+        const int nbr_egde_users);
+
+static void poly_edge_loop_islands_calc(
+        const MEdge *medge, const int totedge, const MPoly *mpoly, const int totpoly,
+        const MLoop *mloop, const int totloop, MeshElemMap *edge_poly_map,
+        const bool use_bitflags, MeshRemap_CheckIslandBoundary edge_boundary_check,
+        int **r_poly_groups, int *r_totgroup, BLI_bitmap **r_edge_borders, int *r_totedgeborder)
 {
 	int *poly_groups;
 	int *poly_stack;
 
+	BLI_bitmap *edge_borders = NULL;
+	int num_edgeborders = 0;
+
 	int poly_prev = 0;
 	const int temp_poly_group_id = 3;  /* Placeholder value. */
 	const int poly_group_id_overflowed = 5;  /* Group we could not find any available bit, will be reset to 0 at end */
@@ -372,18 +407,27 @@ int *BKE_mesh_calc_smoothgroups(const MEdge *medge, const int totedge,
 	bool group_id_overflow = false;
 
 	/* map vars */
-	MeshElemMap *edge_poly_map;
-	int *edge_poly_mem;
+	int *edge_poly_mem = NULL;
 
 	if (totpoly == 0) {
 		*r_totgroup = 0;
-		return NULL;
+		*r_poly_groups = NULL;
+		if (r_edge_borders) {
+			*r_edge_borders = NULL;
+			*r_totedgeborder = 0;
+		}
+		return;
 	}
 
-	BKE_mesh_edge_poly_map_create(&edge_poly_map, &edge_poly_mem,
-	                              medge, totedge,
-	                              mpoly, totpoly,
-	                              mloop, totloop);
+	if (r_edge_borders) {
+		edge_borders = BLI_BITMAP_NEW(totedge, __func__);
+		*r_totedgeborder = 0;
+	}
+
+	if (!edge_poly_map) {
+		BKE_mesh_edge_poly_map_create(&edge_poly_map, &edge_poly_mem,
+		                              medge, totedge, mpoly, totpoly, mloop, totloop);
+	}
 
 	poly_groups = MEM_callocN(sizeof(int) * (size_t)totpoly, __func__);
 	poly_stack  = MEM_mallocN(sizeof(int) * (size_t)totpoly, __func__);
@@ -416,22 +460,20 @@ int *BKE_mesh_calc_smoothgroups(const MEdge *medge, const int totedge,
 		while (ps_curr_idx != ps_end_idx) {
 			const MPoly *mp;
 			const MLoop *ml;
-			bool sharp_poly;
 			int j;
 
 			poly = poly_stack[ps_curr_idx++];
 			BLI_assert(poly_groups[poly] == poly_group_id);
 
 			mp = &mpoly[poly];
-			sharp_poly = !(mp->flag & ME_SMOOTH);
 			for (ml = &mloop[mp->loopstart], j = mp->totloop; j--; ml++) {
 				/* loop over poly users */
-				const MeshElemMap *map_ele = &edge_poly_map[ml->e];
+				const int me_idx = (int)ml->e;
+				const MEdge *me = &medge[me_idx];
+				const MeshElemMap *map_ele = &edge_poly_map[me_idx];
 				const int *p = map_ele->indices;
 				int i = map_ele->count;
-				/* Edge is smooth only if its poly is not sharp, edge is not sharp,
-				 * and edge is used by exactly two polygons. */
-				if (!sharp_poly && !(medge[ml->e].flag & ME_SHARP) && i == 2) {
+				if (!edge_boundary_check(mp, ml, me, i)) {
 					for (; i--; p++) {
 						/* if we meet other non initialized its a bug */
 						BLI_assert(ELEM(poly_groups[*p], 0, poly_group_id));
@@ -442,14 +484,20 @@ int *BKE_mesh_calc_smoothgroups(const MEdge *medge, const int totedge,
 						}
 					}
 				}
-				else if (use_bitflags) {
-					/* Find contiguous smooth groups already assigned, these are the values we can't reuse! */
-					for (; i--; p++) {
-						int bit = poly_groups[*p];
-						if (!ELEM(bit, 0, poly_group_id, poly_group_id_overflowed) &&
-						    !(bit_poly_group_mask & bit))
-						{
-							bit_poly_group_mask |= bit;
+				else {
+					if (edge_borders && !BLI_BITMAP_TEST(edge_borders, me_idx)) {
+						BLI_BITMAP_ENABLE(edge_borders, me_idx);
+						num_edgeborders++;
+					}
+					if (use_bitflags) {
+						/* Find contiguous smooth groups already assigned, these are the values we can't reuse! */
+						for (; i--; p++) {
+							int bit = poly_groups[*p];
+							if (!ELEM(bit, 0, poly_group_id, poly_group_id_overflowed) &&
+							    !(bit_poly_group_mask & bit))
+							{
+								bit_poly_group_mask |= bit;
+							}
 						}
 					}
 				}
@@ -502,12 +550,255 @@ int *BKE_mesh_calc_smoothgroups(const MEdge *medge, const int totedge,
 		tot_group++;
 	}
 
-	MEM_freeN(edge_poly_map);
-	MEM_freeN(edge_poly_mem);
+	if (edge_poly_mem) {
+		MEM_freeN(edge_poly_map);
+		MEM_freeN(edge_poly_mem);
+	}
 	MEM_freeN(poly_stack);
 
 	*r_totgroup = tot_group;
+	*r_poly_groups = poly_groups;
+	if (r_edge_borders) {
+		*r_edge_borders = edge_borders;
+		*r_totedgeborder = num_edgeborders;
+	}
+}
+
+static bool poly_is_island_boundary_smooth_cb(
+        const MPoly *mp, const MLoop *UNUSED(ml), const MEdge *me,
+        const int nbr_egde_users)
+{
+	/* Edge is sharp if its poly is sharp, or edge itself is sharp, or edge is not used by exactly two polygons. */
+	return (!(mp->flag & ME_SMOOTH) || (me->flag & ME_SHARP) || (nbr_egde_users != 2));
+}
+
+/**
+ * Calculate smooth groups from sharp edges.
+ *
+ * \param r_totgroup The total number of groups, 1 or more.
+ * \return Polygon aligned array of group index values (bitflags if use_bitflags is true), starting at 1
+ *         (0 being used as 'invalid' flag).
+ *         Note it's callers's responsibility to MEM_freeN returned array.
+ */
+int *BKE_mesh_calc_smoothgroups(const MEdge *medge, const int totedge,
+                                const MPoly *mpoly, const int totpoly,
+                                const MLoop *mloop, const int totloop,
+                                int *r_totgroup, const bool use_bitflags)
+{
+	int *poly_groups = NULL;
+
+	poly_edge_loop_islands_calc(
+	        medge, totedge, mpoly, totpoly, mloop, totloop, NULL, use_bitflags,
+	        poly_is_island_boundary_smooth_cb, &poly_groups, r_totgroup, NULL, NULL);
 
 	return poly_groups;
 }
+
+#define MISLAND_DEFAULT_BUFSIZE 64
+
+void BKE_mesh_loop_islands_init(
+        MeshIslandStore *island_store,
+        const short item_type, const int items_num, const short island_type, const short innercut_type)
+{
+	MemArena *mem = island_store->mem;
+
+	if (mem == NULL) {
+		mem = BLI_memarena_new(BLI_MEMARENA_STD_BUFSIZE, __func__);
+		island_store->mem = mem;
+	}
+	/* else memarena should be cleared */
+
+	BLI_assert(ELEM(item_type, MISLAND_TYPE_VERT, MISLAND_TYPE_EDGE, MISLAND_TYPE_POLY, MISLAND_TYPE_LOOP));
+	BLI_assert(ELEM(island_type, MISLAND_TYPE_VERT, MISLAND_TYPE_EDGE, MISLAND_TYPE_POLY, MISLAND_TYPE_LOOP));
+
+	island_store->item_type = item_type;
+	island_store->items_to_islands_num = items_num;
+	island_store->items_to_islands = BLI_memarena_alloc(mem, sizeof(*island_store->items_to_islands) * (size_t)items_num);
+
+	island_store->island_type = island_type;
+	island_store->islands_num_alloc = MISLAND_DEFAULT_BUFSIZE;
+	island_store->islands = BLI_memarena_alloc(mem, sizeof(*island_store->islands) * island_store->islands_num_alloc);
+
+	island_store->innercut_type = innercut_type;
+	island_store->innercuts = BLI_memarena_alloc(mem, sizeof(*island_store->innercuts) * island_store->islands_num_alloc);
+}
+
+void BKE_mesh_loop_islands_clear(MeshIslandStore *island_store)
+{
+	island_store->item_type = MISLAND_TYPE_NONE;
+	island_store->items_to_islands_num = 0;
+	island_store->items_to_islands = NULL;
+
+	island_store->island_type = MISLAND_TYPE_NONE;
+	island_store->islands_num = 0;
+	island_store->islands = NULL;
+
+	island_store->innercut_type = MISLAND_TYPE_NONE;
+	island_store->innercuts = NULL;
+
+	if (island_store->mem) {
+		BLI_memarena_clear(island_store->mem);
+	}
+
+	island_store->islands_num_alloc = 0;
+}
+
+void BKE_mesh_loop_islands_free(MeshIslandStore *island_store)
+{
+	if (island_store->mem) {
+		BLI_memarena_free(island_store->mem);
+		island_store->mem = NULL;
+	}
+}
+
+void BKE_mesh_loop_islands_add(
+        MeshIslandStore *island_store, const int item_num, int *items_indices,
+        const int num_island_items, int *island_item_indices,
+        const int num_innercut_items, int *innercut_item_indices)
+{
+	MemArena *mem = island_store->mem;
+
+	MeshElemMap *isld, *innrcut;
+	const int curr_island_idx = island_store->islands_num++;
+	const size_t curr_num_islands = (size_t)island_store->islands_num;
+	int i = item_num;
+
+	island_store->items_to_islands_num = item_num;
+	while (i--) {
+		island_store->items_to_islands[items_indices[i]] = curr_island_idx;
+	}
+
+	if (UNLIKELY(curr_num_islands > island_store->islands_num_alloc)) {
+		MeshElemMap **islds, **innrcuts;
+
+		island_store->islands_num_alloc *= 2;
+		islds = BLI_memarena_alloc(mem, sizeof(*islds) * island_store->islands_num_alloc);
+		memcpy(islds, island_store->islands, sizeof(*islds) * (curr_num_islands - 1));
+		island_store->islands = islds;
+
+		innrcuts = BLI_memarena_alloc(mem, sizeof(*innrcuts) * island_store->islands_num_alloc);
+		memcpy(innrcuts, island_store->innercuts, sizeof(*innrcuts) * (curr_num_islands - 1));
+		island_store->innercuts = innrcuts;
+	}
+
+	island_store->islands[curr_island_idx] = isld = BLI_memarena_alloc(mem, sizeof(*isld));
+	isld->count = num_island_items;
+	isld->indices = BLI_memarena_alloc(mem, sizeof(*isld->indices) * (size_t)num_island_items);
+	memcpy(isld->indices, island_item_indices, sizeof(*isld->indices) * (size_t)num_island_items);
+
+	island_store->innercuts[curr_island_idx] = innrcut = BLI_memarena_alloc(mem, sizeof(*innrcut));
+	innrcut->count = num_innercut_items;
+	innrcut->indices = BLI_memarena_alloc(mem, sizeof(*innrcut->indices) * (size_t)num_innercut_items);
+	memcpy(innrcut->indices, innercut_item_indices, sizeof(*innrcut->indices) * (size_t)num_innercut_items);
+}
+
+/* TODO: I'm not sure edge seam flag is enough to define UV islands? Maybe we should also consider UVmaps values
+ *       themselves (i.e. different UV-edges for a same mesh-edge => boundary edge too?).
+ *       Would make things much more complex though, and each UVMap would then need its own mesh mapping,
+ *       not sure we want that at all!
+ */
+static bool mesh_check_island_boundary_uv(
+        const MPoly *UNUSED(mp), const MLoop *UNUSED(ml), const MEdge *me,
+        const int UNUSED(nbr_egde_users))
+{
+	/* Edge is UV boundary if tagged as seam. */
+	return (me->flag & ME_SEAM) != 0;
+}
+
+/**
+ * \note all this could be optimized...
+ * Not sure it would be worth the more complex code, though, those loops
+ * are supposed to be really quick to do...
+ */
+bool BKE_mesh_calc_islands_loop_poly_uv(
+        MVert *UNUSED(verts), const int UNUSED(totvert),
+        MEdge *edges, const int totedge,
+        MPoly *polys, const int totpoly,
+        MLoop *loops, const int totloop,
+        MeshIslandStore *r_island_store)
+{
+	int *poly_groups = NULL;
+	int num_poly_groups;
+
+	/* map vars */
+	MeshElemMap *edge_poly_map;
+	int *edge_poly_mem;
+
+	int *poly_indices = MEM_mallocN(sizeof(*poly_indices) * (size_t)totpoly, __func__);
+	int *loop_indices = MEM_mallocN(sizeof(*loop_indices) * (size_t)totloop, __func__);
+	int num_pidx, num_lidx;
+
+	/* Those are used to detect 'inner cuts', i.e. edges that are borders, and yet have two or more polys of
+	 * a same group using them (typical case: seam used to unwrap properly a cylinder). */
+	BLI_bitmap *edge_borders;
+	int num_edge_borders;
+	char *edge_border_count = NULL;
+	int *edge_innercut_indices = NULL;
+	int num_einnercuts = 0;
+
+	int grp_idx, p_idx, pl_idx, l_idx;
+
+	BKE_mesh_loop_islands_clear(r_island_store);
+	BKE_mesh_loop_islands_init(r_island_store, MISLAND_TYPE_LOOP, totloop, MISLAND_TYPE_POLY, MISLAND_TYPE_EDGE);
+
+	BKE_mesh_edge_poly_map_create(&edge_poly_map, &edge_poly_mem,
+	                              edges, totedge, polys, totpoly, loops, totloop);
+
+	poly_edge_loop_islands_calc(
+	        edges, totedge, polys, totpoly, loops, totloop, edge_poly_map, false,
+	        mesh_check_island_boundary_uv, &poly_groups, &num_poly_groups, &edge_borders, &num_edge_borders);
+
+	if (!num_poly_groups) {
+		/* Should never happen... */
+		return false;
+	}
+
+	if (num_edge_borders) {
+		edge_border_count = MEM_mallocN(sizeof(*edge_border_count) * (size_t)totedge, __func__);
+		edge_innercut_indices = MEM_mallocN(sizeof(*edge_innercut_indices) * (size_t)num_edge_borders, __func__);
+	}
+
+	/* Note: here we ignore '0' invalid group - this should *never* happen in this case anyway? */
+	for (grp_idx = 1; grp_idx <= num_poly_groups; grp_idx++) {
+		num_pidx = num_lidx = 0;
+		if (num_edge_borders) {
+			num_einnercuts = 0;
+			memset(edge_border_count, 0, sizeof(*edge_border_count) * (size_t)totedge);
+		}
+
+		for (p_idx = 0; p_idx < totpoly; p_idx++) {
+			MPoly *mp;
+
+			if (poly_groups[p_idx] != grp_idx) {
+				continue;
+			}
+
+			mp = &polys[p_idx];
+			poly_indices[num_pidx++] = p_idx;
+			for (l_idx = mp->loopstart, pl_idx = 0; pl_idx < mp->totloop; l_idx++, pl_idx++) {
+				MLoop *ml = &loops[l_idx];
+				loop_indices[num_lidx++] = l_idx;
+				if (num_edge_borders && BLI_BITMAP_TEST(edge_borders, ml->e) && (edge_border_count[ml->e] < 2)) {
+					edge_border_count[ml->e]++;
+					if (edge_border_count[ml->e] == 2) {
+						edge_innercut_indices[num_einnercuts++] = (int)ml->e;
+					}
+				}
+			}
+		}
+
+		BKE_mesh_loop_islands_add(r_island_store, num_lidx, loop_indices, num_pidx, poly_indices,
+		                          num_einnercuts, edge_innercut_indices);
+	}
+
+	MEM_freeN(poly_indices);
+	MEM_freeN(loop_indices);
+	MEM_freeN(poly_groups);
+	if (num_edge_borders) {
+		MEM_freeN(edge_border_count);
+		MEM_freeN(edge_innercut_indices);
+	}
+	return true;
+}
+
 /** \} */
diff --git a/source/blender/blenkernel/intern/mesh_remap.c b/source/blender/blenkernel/intern/mesh_remap.c
new file mode 100644
index 00000000000..73ddb18d7b3
--- /dev/null
+++ b/source/blender/blenkernel/intern/mesh_remap.c
@@ -0,0 +1,1986 @@
+/*
+ * ***** BEGIN GPL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ */
+
+/** \file blender/blenkernel/intern/mesh_remap.c
+ *  \ingroup bke
+ *
+ * Functions for mapping data between meshes.
+ */
+
+#include <limits.h>
+
+#include "MEM_guardedalloc.h"
+
+#include "DNA_meshdata_types.h"
+
+#include "BLI_utildefines.h"
+#include "BLI_alloca.h"
+#include "BLI_astar.h"
+#include "BLI_bitmap.h"
+#include "BLI_math.h"
+#include "BLI_memarena.h"
+#include "BLI_polyfill2d.h"
+#include "BLI_rand.h"
+
+#include "BKE_bvhutils.h"
+#include "BKE_customdata.h"
+#include "BKE_DerivedMesh.h"
+#include "BKE_mesh.h"
+#include "BKE_mesh_mapping.h"
+#include "BKE_mesh_remap.h"  /* own include */
+
+#include "BLI_strict_flags.h"
+
+
+/* -------------------------------------------------------------------- */
+
+/** \name Mesh to mesh mapping
+ * \{ */
+
+void BKE_mesh_remap_init(MeshPairRemap *map, const int items_num)
+{
+	MemArena *mem = BLI_memarena_new(BLI_MEMARENA_STD_BUFSIZE, __func__);
+
+	BKE_mesh_remap_free(map);
+
+	map->items = BLI_memarena_alloc(mem, sizeof(*map->items) * (size_t)items_num);
+	map->items_num = items_num;
+
+	map->mem = mem;
+}
+
+void BKE_mesh_remap_free(MeshPairRemap *map)
+{
+	if (map->mem) {
+		BLI_memarena_free((MemArena *)map->mem);
+	}
+
+	map->items_num = 0;
+	map->items = NULL;
+	map->mem = NULL;
+}
+
+static void mesh_remap_item_define(
+        MeshPairRemap *map, const int index, const float UNUSED(hit_dist), const int island,
+        const int sources_num, const int *indices_src, const float *weights_src)
+{
+	MeshPairRemapItem *mapit = &map->items[index];
+	MemArena *mem = map->mem;
+
+	if (sources_num) {
+		mapit->sources_num = sources_num;
+		mapit->indices_src = BLI_memarena_alloc(mem, sizeof(*mapit->indices_src) * (size_t)sources_num);
+		memcpy(mapit->indices_src, indices_src, sizeof(*mapit->indices_src) * (size_t)sources_num);
+		mapit->weights_src = BLI_memarena_alloc(mem, sizeof(*mapit->weights_src) * (size_t)sources_num);
+		memcpy(mapit->weights_src, weights_src, sizeof(*mapit->weights_src) * (size_t)sources_num);
+	}
+	else {
+		mapit->sources_num = 0;
+		mapit->indices_src = NULL;
+		mapit->weights_src = NULL;
+	}
+	/* UNUSED */
+	// mapit->hit_dist = hit_dist;
+	mapit->island = island;
+}
+
+void BKE_mesh_remap_item_define_invalid(MeshPairRemap *map, const int index)
+{
+	mesh_remap_item_define(map, index, FLT_MAX, 0, 0, NULL, NULL);
+}
+
+static int mesh_remap_interp_poly_data_get(
+        const MPoly *mp, MLoop *mloops, const float (*vcos_src)[3], const float point[3],
+        size_t *buff_size, float (**vcos)[3], const bool use_loops, int **indices, float **weights,
+        const bool do_weights, int *r_closest_index)
+{
+	MLoop *ml;
+	float (*vco)[3];
+	float ref_dist_sq = FLT_MAX;
+	int *index;
+	const int sources_num = mp->totloop;
+	int i;
+
+	if ((size_t)sources_num > *buff_size) {
+		*buff_size = (size_t)sources_num;
+		*vcos = MEM_reallocN(*vcos, sizeof(**vcos) * *buff_size);
+		*indices = MEM_reallocN(*indices, sizeof(**indices) * *buff_size);
+		if (do_weights) {
+			*weights = MEM_reallocN(*weights, sizeof(**weights) * *buff_size);
+		}
+	}
+
+	for (i = 0, ml = &mloops[mp->loopstart], vco = *vcos, index = *indices; i < sources_num; i++, ml++, vco++, index++) {
+		*index = use_loops ? (int)mp->loopstart + i : (int)ml->v;
+		copy_v3_v3(*vco, vcos_src[ml->v]);
+		if (r_closest_index) {
+			/* Find closest vert/loop in this case. */
+			const float dist_sq = len_squared_v3v3(point, *vco);
+			if (dist_sq < ref_dist_sq) {
+				ref_dist_sq = dist_sq;
+				*r_closest_index = *index;
+			}
+		}
+	}
+
+	if (do_weights) {
+		interp_weights_poly_v3(*weights, *vcos, sources_num, point);
+	}
+
+	return sources_num;
+}
+
+static bool mesh_remap_bvhtree_query_nearest(
+        BVHTreeFromMesh *treedata, BVHTreeNearest *nearest, const SpaceTransform *space_transform,
+        const float in_co[3], const float max_dist_sq,
+        float *r_hit_dist)
+{
+	float co[3];
+
+	copy_v3_v3(co, in_co);
+
+	/* Convert the vertex to tree coordinates, if needed. */
+	if (space_transform) {
+		BLI_space_transform_apply(space_transform, co);
+	}
+
+	/* Use local proximity heuristics (to reduce the nearest search). */
+	if (nearest->index != -1) {
+		nearest->dist_sq = min_ff(len_squared_v3v3(co, nearest->co), max_dist_sq);
+	}
+	else {
+		nearest->dist_sq = max_dist_sq;
+	}
+	/* Compute and store result. If invalid (-1 index), keep FLT_MAX dist. */
+	BLI_bvhtree_find_nearest(treedata->tree, co, nearest, treedata->nearest_callback, treedata);
+
+	if ((nearest->index != -1) && (nearest->dist_sq <= max_dist_sq)) {
+		*r_hit_dist = sqrtf(nearest->dist_sq);
+		return true;
+	}
+	else {
+		return false;
+	}
+}
+
+static bool mesh_remap_bvhtree_query_raycast(
+        BVHTreeFromMesh *treedata, BVHTreeRayHit *rayhit, const SpaceTransform *space_transform,
+        const float in_co[3], const float in_no[3], const float radius, const float max_dist,
+        float *r_hit_dist)
+{
+	BVHTreeRayHit rayhit_tmp;
+	float co[3], no[3];
+
+	copy_v3_v3(co, in_co);
+	copy_v3_v3(no, in_no);
+
+	/* Convert the vertex to tree coordinates, if needed. */
+	if (space_transform) {
+		BLI_space_transform_apply(space_transform, co);
+		BLI_space_transform_apply_normal(space_transform, no);
+	}
+
+	rayhit->index = -1;
+	rayhit->dist = max_dist;
+	BLI_bvhtree_ray_cast(treedata->tree, co, no, radius, rayhit, treedata->raycast_callback, treedata);
+
+	/* Also cast in the other direction! */
+	rayhit_tmp = *rayhit;
+	negate_v3(no);
+	BLI_bvhtree_ray_cast(treedata->tree, co, no, radius, &rayhit_tmp, treedata->raycast_callback, treedata);
+	if (rayhit_tmp.dist < rayhit->dist) {
+		*rayhit = rayhit_tmp;
+	}
+
+	if ((rayhit->index != -1) && (rayhit->dist <= max_dist)) {
+		*r_hit_dist = rayhit->dist;
+		return true;
+	}
+	else {
+		return false;
+	}
+}
+
+/* Little helper when dealing with source islands */
+typedef struct IslandResult {
+	float factor;           /* A factor, based on which best island for a given set of elements will be selected. */
+	int   index_src;        /* Index of the source. */
+	float hit_dist;         /* The actual hit distance. */
+	float hit_point[3];     /* The hit point, if relevant. */
+} IslandResult;
+
+/* Note about all bvh/raycasting stuff below:
+ *      * We must use our ray radius as BVH epsilon too, else rays not hitting anything but 'passing near' an item
+ *        would be missed (since BVH handling would not detect them, 'refining' callbacks won't be executed,
+ *        even though they would return a valid hit).
+ *      * However, in 'islands' case where each hit gets a weight, 'precise' hits should have a better weight than
+ *        'approximate' hits. To address that, we simplify things with:
+ *        ** A first raycast with default, given rayradius;
+ *        ** If first one fails, we do more raycasting with bigger radius, but if hit is found
+ *           it will get smaller weight.
+ *        This only concerns loops, currently (because of islands), and 'sampled' edges/polys norproj.
+ */
+
+/* At most n raycasts per 'real' ray. */
+#define MREMAP_RAYCAST_APPROXIMATE_NR 3
+/* Each approximated raycasts will have n times bigger radius than previous one. */
+#define MREMAP_RAYCAST_APPROXIMATE_FAC 5.0f
+/* BVH epsilon value we have to give to bvh 'constructor' when doing approximated raycasting. */
+#define MREMAP_RAYCAST_APPROXIMATE_BVHEPSILON(_ray_radius) \
+	((float)MREMAP_RAYCAST_APPROXIMATE_NR * MREMAP_RAYCAST_APPROXIMATE_FAC * (_ray_radius))
+
+/* min 16 rays/face, max 400. */
+#define MREMAP_RAYCAST_TRI_SAMPLES_MIN 4
+#define MREMAP_RAYCAST_TRI_SAMPLES_MAX 20
+
+/* Will be enough in 99% of cases. */
+#define MREMAP_DEFAULT_BUFSIZE 32
+
+void BKE_mesh_remap_calc_verts_from_dm(
+        const int mode, const SpaceTransform *space_transform, const float max_dist, const float ray_radius,
+        const MVert *verts_dst, const int numverts_dst, const bool UNUSED(dirty_nors_dst), DerivedMesh *dm_src,
+        MeshPairRemap *r_map)
+{
+	const float full_weight = 1.0f;
+	const float max_dist_sq = max_dist * max_dist;
+	int i;
+
+	BLI_assert(mode & MREMAP_MODE_VERT);
+
+	BKE_mesh_remap_init(r_map, numverts_dst);
+
+	if (mode == MREMAP_MODE_TOPOLOGY) {
+		BLI_assert(numverts_dst == dm_src->getNumVerts(dm_src));
+		for (i = 0; i < numverts_dst; i++) {
+			mesh_remap_item_define(r_map, i, FLT_MAX, 0, 1, &i, &full_weight);
+		}
+	}
+	else {
+		BVHTreeFromMesh treedata = {NULL};
+		BVHTreeNearest nearest = {0};
+		BVHTreeRayHit rayhit = {0};
+		float hit_dist;
+
+		if (mode == MREMAP_MODE_VERT_NEAREST) {
+			bvhtree_from_mesh_verts(&treedata, dm_src, 0.0f, 2, 6);
+			nearest.index = -1;
+
+			for (i = 0; i < numverts_dst; i++) {
+				float tmp_co[3];
+
+				copy_v3_v3(tmp_co, verts_dst[i].co);
+
+				if (mesh_remap_bvhtree_query_nearest(
+				        &treedata, &nearest, space_transform,
+				        tmp_co, max_dist_sq, &hit_dist))
+				{
+					mesh_remap_item_define(r_map, i, hit_dist, 0, 1, &nearest.index, &full_weight);
+				}
+				else {
+					/* No source for this dest vertex! */
+					BKE_mesh_remap_item_define_invalid(r_map, i);
+				}
+			}
+		}
+		else if (ELEM(mode, MREMAP_MODE_VERT_EDGE_NEAREST, MREMAP_MODE_VERT_EDGEINTERP_NEAREST)) {
+			MEdge *edges_src = dm_src->getEdgeArray(dm_src);
+			float (*vcos_src)[3] = MEM_mallocN(sizeof(*vcos_src) * (size_t)dm_src->getNumVerts(dm_src), __func__);
+			dm_src->getVertCos(dm_src, vcos_src);
+
+			bvhtree_from_mesh_edges(&treedata, dm_src, 0.0f, 2, 6);
+			nearest.index = -1;
+
+			for (i = 0; i < numverts_dst; i++) {
+				float tmp_co[3];
+
+				copy_v3_v3(tmp_co, verts_dst[i].co);
+
+				if (mesh_remap_bvhtree_query_nearest(
+				        &treedata, &nearest, space_transform,
+				        tmp_co, max_dist_sq, &hit_dist))
+				{
+					MEdge *me = &edges_src[nearest.index];
+					const float *v1cos = vcos_src[me->v1];
+					const float *v2cos = vcos_src[me->v2];
+
+					if (mode == MREMAP_MODE_VERT_EDGE_NEAREST) {
+						const float dist_v1 = len_squared_v3v3(tmp_co, v1cos);
+						const float dist_v2 = len_squared_v3v3(tmp_co, v2cos);
+						const int index = (int)((dist_v1 > dist_v2) ? me->v2 : me->v1);
+						mesh_remap_item_define(r_map, i, hit_dist, 0, 1, &index, &full_weight);
+					}
+					else if (mode == MREMAP_MODE_VERT_EDGEINTERP_NEAREST) {
+						int indices[2];
+						float weights[2];
+
+						indices[0] = (int)me->v1;
+						indices[1] = (int)me->v2;
+
+						/* Weight is inverse of point factor here... */
+						weights[0] = line_point_factor_v3(tmp_co, v2cos, v1cos);
+						CLAMP(weights[0], 0.0f, 1.0f);
+						weights[1] = 1.0f - weights[0];
+
+						mesh_remap_item_define(r_map, i, hit_dist, 0, 2, indices, weights);
+					}
+				}
+				else {
+					/* No source for this dest vertex! */
+					BKE_mesh_remap_item_define_invalid(r_map, i);
+				}
+			}
+
+			MEM_freeN(vcos_src);
+		}
+		else if (ELEM(mode, MREMAP_MODE_VERT_POLY_NEAREST, MREMAP_MODE_VERT_POLYINTERP_NEAREST,
+		                    MREMAP_MODE_VERT_POLYINTERP_VNORPROJ))
+		{
+			MPoly *polys_src = dm_src->getPolyArray(dm_src);
+			MLoop *loops_src = dm_src->getLoopArray(dm_src);
+			float (*vcos_src)[3] = MEM_mallocN(sizeof(*vcos_src) * (size_t)dm_src->getNumVerts(dm_src), __func__);
+			int *tessface_to_poly_map_src;
+
+			size_t tmp_buff_size = MREMAP_DEFAULT_BUFSIZE;
+			float (*vcos)[3] = MEM_mallocN(sizeof(*vcos) * tmp_buff_size, __func__);
+			int *indices = MEM_mallocN(sizeof(*indices) * tmp_buff_size, __func__);
+			float *weights = MEM_mallocN(sizeof(*weights) * tmp_buff_size, __func__);
+
+			dm_src->getVertCos(dm_src, vcos_src);
+			bvhtree_from_mesh_faces(&treedata, dm_src, (mode & MREMAP_USE_NORPROJ) ? ray_radius : 0.0f, 2, 6);
+			/* bvhtree here uses tesselated faces... */
+			tessface_to_poly_map_src = dm_src->getTessFaceDataArray(dm_src, CD_ORIGINDEX);
+
+			if (mode == MREMAP_MODE_VERT_POLYINTERP_VNORPROJ) {
+				for (i = 0; i < numverts_dst; i++) {
+					float tmp_co[3], tmp_no[3];
+
+					copy_v3_v3(tmp_co, verts_dst[i].co);
+					normal_short_to_float_v3(tmp_no, verts_dst[i].no);
+
+					if (mesh_remap_bvhtree_query_raycast(
+					        &treedata, &rayhit, space_transform,
+					        tmp_co, tmp_no, ray_radius, max_dist, &hit_dist))
+					{
+						MPoly *mp_src = &polys_src[tessface_to_poly_map_src[rayhit.index]];
+						const int sources_num = mesh_remap_interp_poly_data_get(
+						        mp_src, loops_src, (const float (*)[3])vcos_src, rayhit.co,
+						        &tmp_buff_size, &vcos, false, &indices, &weights, true, NULL);
+
+						mesh_remap_item_define(r_map, i, hit_dist, 0, sources_num, indices, weights);
+					}
+					else {
+						/* No source for this dest vertex! */
+						BKE_mesh_remap_item_define_invalid(r_map, i);
+					}
+				}
+			}
+			else {
+				nearest.index = -1;
+
+				for (i = 0; i < numverts_dst; i++) {
+					float tmp_co[3];
+
+					/* Convert the vertex to tree coordinates. */
+					copy_v3_v3(tmp_co, verts_dst[i].co);
+
+					if (mesh_remap_bvhtree_query_nearest(
+					        &treedata, &nearest, space_transform,
+					        tmp_co, max_dist_sq, &hit_dist))
+					{
+						MPoly *mp = &polys_src[tessface_to_poly_map_src[nearest.index]];
+
+						if (mode == MREMAP_MODE_VERT_POLY_NEAREST) {
+							int index;
+							mesh_remap_interp_poly_data_get(
+							        mp, loops_src, (const float (*)[3])vcos_src, nearest.co,
+							        &tmp_buff_size, &vcos, false, &indices, &weights, false,
+							        &index);
+
+							mesh_remap_item_define(r_map, i, hit_dist, 0, 1, &index, &full_weight);
+						}
+						else if (mode == MREMAP_MODE_VERT_POLYINTERP_NEAREST) {
+							const int sources_num = mesh_remap_interp_poly_data_get(
+							        mp, loops_src, (const float (*)[3])vcos_src, nearest.co,
+							        &tmp_buff_size, &vcos, false, &indices, &weights, true,
+							        NULL);
+
+							mesh_remap_item_define(r_map, i, hit_dist, 0, sources_num, indices, weights);
+						}
+					}
+					else {
+						/* No source for this dest vertex! */
+						BKE_mesh_remap_item_define_invalid(r_map, i);
+					}
+				}
+			}
+
+			MEM_freeN(vcos_src);
+			MEM_freeN(vcos);
+			MEM_freeN(indices);
+			MEM_freeN(weights);
+		}
+		else {
+			printf("WARNING! Unsupported mesh-to-mesh vertex mapping mode (%d)!\n", mode);
+			memset(r_map->items, 0, sizeof(*r_map->items) * (size_t)numverts_dst);
+		}
+
+		free_bvhtree_from_mesh(&treedata);
+	}
+}
+
+void BKE_mesh_remap_calc_edges_from_dm(
+        const int mode, const SpaceTransform *space_transform, const float max_dist, const float ray_radius,
+        const MVert *verts_dst, const int numverts_dst, const MEdge *edges_dst, const int numedges_dst,
+        const bool UNUSED(dirty_nors_dst), DerivedMesh *dm_src, MeshPairRemap *r_map)
+{
+	const float full_weight = 1.0f;
+	const float max_dist_sq = max_dist * max_dist;
+	int i;
+
+	BLI_assert(mode & MREMAP_MODE_EDGE);
+
+	BKE_mesh_remap_init(r_map, numedges_dst);
+
+	if (mode == MREMAP_MODE_TOPOLOGY) {
+		BLI_assert(numedges_dst == dm_src->getNumEdges(dm_src));
+		for (i = 0; i < numedges_dst; i++) {
+			mesh_remap_item_define(r_map, i, FLT_MAX, 0, 1, &i, &full_weight);
+		}
+	}
+	else {
+		BVHTreeFromMesh treedata = {NULL};
+		BVHTreeNearest nearest = {0};
+		BVHTreeRayHit rayhit = {0};
+		float hit_dist;
+
+		if (mode == MREMAP_MODE_EDGE_VERT_NEAREST) {
+			const int num_verts_src = dm_src->getNumVerts(dm_src);
+			const int num_edges_src = dm_src->getNumEdges(dm_src);
+			MEdge *edges_src = dm_src->getEdgeArray(dm_src);
+			float (*vcos_src)[3] = MEM_mallocN(sizeof(*vcos_src) * (size_t)dm_src->getNumVerts(dm_src), __func__);
+
+			MeshElemMap *vert_to_edge_src_map;
+			int         *vert_to_edge_src_map_mem;
+
+			struct {
+				float hit_dist;
+				int   index;
+			} *v_dst_to_src_map = MEM_mallocN(sizeof(*v_dst_to_src_map) * (size_t)numverts_dst, __func__);
+
+			for (i = 0; i < numverts_dst; i++) {
+				v_dst_to_src_map[i].hit_dist = -1.0f;
+			}
+
+			BKE_mesh_vert_edge_map_create(&vert_to_edge_src_map, &vert_to_edge_src_map_mem,
+			                              edges_src, num_verts_src, num_edges_src);
+
+			dm_src->getVertCos(dm_src, vcos_src);
+
+			bvhtree_from_mesh_verts(&treedata, dm_src, 0.0f, 2, 6);
+			nearest.index = -1;
+
+			for (i = 0; i < numedges_dst; i++) {
+				const MEdge *e_dst = &edges_dst[i];
+				float best_totdist = FLT_MAX;
+				int best_eidx_src = -1;
+				int j = 2;
+
+				while (j--) {
+					const unsigned int vidx_dst = j ? e_dst->v1 : e_dst->v2;
+
+					/* Compute closest verts only once! */
+					if (v_dst_to_src_map[vidx_dst].hit_dist == -1.0f) {
+						float tmp_co[3];
+
+						copy_v3_v3(tmp_co, verts_dst[vidx_dst].co);
+
+						if (mesh_remap_bvhtree_query_nearest(
+						        &treedata, &nearest, space_transform,
+						        tmp_co, max_dist_sq, &hit_dist))
+						{
+							v_dst_to_src_map[vidx_dst].hit_dist = hit_dist;
+							v_dst_to_src_map[vidx_dst].index = nearest.index;
+						}
+						else {
+							/* No source for this dest vert! */
+							v_dst_to_src_map[vidx_dst].hit_dist = FLT_MAX;
+						}
+					}
+				}
+
+				/* Now, check all source edges of closest sources vertices, and select the one giving the smallest
+				 * total verts-to-verts distance. */
+				for (j = 2; j--;) {
+					const unsigned int vidx_dst = j ? e_dst->v1 : e_dst->v2;
+					const float first_dist = v_dst_to_src_map[vidx_dst].hit_dist;
+					const int vidx_src = v_dst_to_src_map[vidx_dst].index;
+					int *eidx_src, k;
+
+					if (vidx_src < 0) {
+						continue;
+					}
+
+					eidx_src = vert_to_edge_src_map[vidx_src].indices;
+					k = vert_to_edge_src_map[vidx_src].count;
+
+					for (; k--; eidx_src++) {
+						MEdge *e_src = &edges_src[*eidx_src];
+						const float *other_co_src = vcos_src[BKE_mesh_edge_other_vert(e_src, vidx_src)];
+						const float *other_co_dst = verts_dst[BKE_mesh_edge_other_vert(e_dst, (int)vidx_dst)].co;
+						const float totdist = first_dist + len_v3v3(other_co_src, other_co_dst);
+
+						if (totdist < best_totdist) {
+							best_totdist = totdist;
+							best_eidx_src = *eidx_src;
+						}
+					}
+				}
+
+				if (best_eidx_src >= 0) {
+					const float *co1_src = vcos_src[edges_src[best_eidx_src].v1];
+					const float *co2_src = vcos_src[edges_src[best_eidx_src].v2];
+					const float *co1_dst = verts_dst[e_dst->v1].co;
+					const float *co2_dst = verts_dst[e_dst->v2].co;
+					float co_src[3], co_dst[3];
+
+					/* TODO: would need an isect_seg_seg_v3(), actually! */
+					const int isect_type = isect_line_line_v3(co1_src, co2_src, co1_dst, co2_dst, co_src, co_dst);
+					if (isect_type != 0) {
+						const float fac_src = line_point_factor_v3(co_src, co1_src, co2_src);
+						const float fac_dst = line_point_factor_v3(co_dst, co1_dst, co2_dst);
+						if (fac_src < 0.0f) {
+							copy_v3_v3(co_src, co1_src);
+						}
+						else if (fac_src > 1.0f) {
+							copy_v3_v3(co_src, co2_src);
+						}
+						if (fac_dst < 0.0f) {
+							copy_v3_v3(co_dst, co1_dst);
+						}
+						else if (fac_dst > 1.0f) {
+							copy_v3_v3(co_dst, co2_dst);
+						}
+					}
+					hit_dist = len_v3v3(co_dst, co_src);
+					mesh_remap_item_define(r_map, i, hit_dist, 0, 1, &best_eidx_src, &full_weight);
+				}
+				else {
+					/* No source for this dest edge! */
+					BKE_mesh_remap_item_define_invalid(r_map, i);
+				}
+			}
+
+			MEM_freeN(vcos_src);
+			MEM_freeN(v_dst_to_src_map);
+			MEM_freeN(vert_to_edge_src_map);
+			MEM_freeN(vert_to_edge_src_map_mem);
+		}
+		else if (mode == MREMAP_MODE_EDGE_NEAREST) {
+			bvhtree_from_mesh_edges(&treedata, dm_src, 0.0f, 2, 6);
+			nearest.index = -1;
+
+			for (i = 0; i < numedges_dst; i++) {
+				float tmp_co[3];
+
+				interp_v3_v3v3(tmp_co, verts_dst[edges_dst[i].v1].co, verts_dst[edges_dst[i].v2].co, 0.5f);
+
+				if (mesh_remap_bvhtree_query_nearest(
+				        &treedata, &nearest, space_transform,
+				        tmp_co, max_dist_sq, &hit_dist))
+				{
+					mesh_remap_item_define(r_map, i, hit_dist, 0, 1, &nearest.index, &full_weight);
+				}
+				else {
+					/* No source for this dest edge! */
+					BKE_mesh_remap_item_define_invalid(r_map, i);
+				}
+			}
+		}
+		else if (mode == MREMAP_MODE_EDGE_POLY_NEAREST) {
+			MEdge *edges_src = dm_src->getEdgeArray(dm_src);
+			MPoly *polys_src = dm_src->getPolyArray(dm_src);
+			MLoop *loops_src = dm_src->getLoopArray(dm_src);
+			float (*vcos_src)[3] = MEM_mallocN(sizeof(*vcos_src) * (size_t)dm_src->getNumVerts(dm_src), __func__);
+			int *tessface_to_poly_map_src;
+
+			dm_src->getVertCos(dm_src, vcos_src);
+			bvhtree_from_mesh_faces(&treedata, dm_src, 0.0f, 2, 6);
+			/* bvhtree here uses tesselated faces... */
+			tessface_to_poly_map_src = dm_src->getTessFaceDataArray(dm_src, CD_ORIGINDEX);
+
+			for (i = 0; i < numedges_dst; i++) {
+				float tmp_co[3];
+
+				interp_v3_v3v3(tmp_co, verts_dst[edges_dst[i].v1].co, verts_dst[edges_dst[i].v2].co, 0.5f);
+
+				if (mesh_remap_bvhtree_query_nearest(
+				        &treedata, &nearest, space_transform,
+				        tmp_co, max_dist_sq, &hit_dist))
+				{
+					MPoly *mp_src = &polys_src[tessface_to_poly_map_src[nearest.index]];
+					MLoop *ml_src = &loops_src[mp_src->loopstart];
+					int nloops = mp_src->totloop;
+					float best_dist_sq = FLT_MAX;
+					int best_eidx_src = -1;
+
+					for (; nloops--; ml_src++) {
+						MEdge *me_src = &edges_src[ml_src->e];
+						float *co1_src = vcos_src[me_src->v1];
+						float *co2_src = vcos_src[me_src->v2];
+						float co_src[3];
+						float dist_sq;
+
+						interp_v3_v3v3(co_src, co1_src, co2_src, 0.5f);
+						dist_sq = len_squared_v3v3(tmp_co, co_src);
+						if (dist_sq < best_dist_sq) {
+							best_dist_sq = dist_sq;
+							best_eidx_src = (int)ml_src->e;
+						}
+					}
+					if (best_eidx_src >= 0) {
+						mesh_remap_item_define(r_map, i, hit_dist, 0, 1, &best_eidx_src, &full_weight);
+					}
+				}
+				else {
+					/* No source for this dest edge! */
+					BKE_mesh_remap_item_define_invalid(r_map, i);
+				}
+			}
+
+			MEM_freeN(vcos_src);
+		}
+		else if (mode == MREMAP_MODE_EDGE_EDGEINTERP_VNORPROJ) {
+			const int num_rays_min = 5, num_rays_max = 100;
+			const int numedges_src = dm_src->getNumEdges(dm_src);
+
+			/* Subtleness - this one we can allocate only max number of cast rays per edges! */
+			int *indices = MEM_mallocN(sizeof(*indices) * (size_t)min_ii(numedges_src, num_rays_max), __func__);
+			/* Here it's simpler to just allocate for all edges :/ */
+			float *weights = MEM_mallocN(sizeof(*weights) * (size_t)numedges_src, __func__);
+
+			bvhtree_from_mesh_edges(&treedata, dm_src, MREMAP_RAYCAST_APPROXIMATE_BVHEPSILON(ray_radius), 2, 6);
+
+			for (i = 0; i < numedges_dst; i++) {
+				/* For each dst edge, we sample some rays from it (interpolated from its vertices)
+				 * and use their hits to interpolate from source edges. */
+				const MEdge *me = &edges_dst[i];
+				float tmp_co[3], v1_co[3], v2_co[3];
+				float tmp_no[3], v1_no[3], v2_no[3];
+
+				int grid_size;
+				float edge_dst_len;
+				float grid_step;
+
+				float totweights = 0.0f;
+				float hit_dist_accum = 0.0f;
+				int sources_num = 0;
+				int j;
+
+				copy_v3_v3(v1_co, verts_dst[me->v1].co);
+				copy_v3_v3(v2_co, verts_dst[me->v2].co);
+
+				normal_short_to_float_v3(v1_no, verts_dst[me->v1].no);
+				normal_short_to_float_v3(v2_no, verts_dst[me->v2].no);
+
+				/* We do our transform here, allows to interpolate from normals already in src space. */
+				if (space_transform) {
+					BLI_space_transform_apply(space_transform, v1_co);
+					BLI_space_transform_apply(space_transform, v2_co);
+					BLI_space_transform_apply_normal(space_transform, v1_no);
+					BLI_space_transform_apply_normal(space_transform, v2_no);
+				}
+
+				fill_vn_fl(weights, (int)numedges_src, 0.0f);
+
+				/* We adjust our ray-casting grid to ray_radius (the smaller, the more rays are cast),
+				 * with lower/upper bounds. */
+				edge_dst_len = len_v3v3(v1_co, v2_co);
+
+				grid_size = (int)((edge_dst_len / ray_radius) + 0.5f);
+				CLAMP(grid_size, num_rays_min, num_rays_max);  /* min 5 rays/edge, max 100. */
+
+				grid_step = 1.0f / (float)grid_size;  /* Not actual distance here, rather an interp fac... */
+
+				/* And now we can cast all our rays, and see what we get! */
+				for (j = 0; j < grid_size; j++) {
+					const float fac = grid_step * (float)j;
+
+					int n = (ray_radius > 0.0f) ? MREMAP_RAYCAST_APPROXIMATE_NR : 1;
+					float w = 1.0f;
+
+					interp_v3_v3v3(tmp_co, v1_co, v2_co, fac);
+					interp_v3_v3v3_slerp_safe(tmp_no, v1_no, v2_no, fac);
+
+					while (n--) {
+						/* Note we handle dest to src space conversion ourself, here! */
+						if (mesh_remap_bvhtree_query_raycast(
+						        &treedata, &rayhit, NULL,
+						        tmp_co, tmp_no, ray_radius / w, max_dist, &hit_dist))
+						{
+							weights[rayhit.index] += w;
+							totweights += w;
+							hit_dist_accum += hit_dist;
+							break;
+						}
+						/* Next iteration will get bigger radius but smaller weight! */
+						w /= MREMAP_RAYCAST_APPROXIMATE_FAC;
+					}
+				}
+				/* A sampling is valid (as in, its result can be considered as valid sources) only if at least
+				 * half of the rays found a source! */
+				if (totweights > ((float)grid_size / 2.0f)) {
+					for (j = 0; j < (int)numedges_src; j++) {
+						if (!weights[j]) {
+							continue;
+						}
+						/* Note: sources_num is always <= j! */
+						weights[sources_num] = weights[j] / totweights;
+						indices[sources_num] = j;
+						sources_num++;
+					}
+					mesh_remap_item_define(r_map, i, hit_dist_accum / totweights, 0,
+					                                  sources_num, indices, weights);
+				}
+				else {
+					/* No source for this dest edge! */
+					BKE_mesh_remap_item_define_invalid(r_map, i);
+				}
+			}
+
+			MEM_freeN(indices);
+			MEM_freeN(weights);
+		}
+		else {
+			printf("WARNING! Unsupported mesh-to-mesh edge mapping mode (%d)!\n", mode);
+			memset(r_map->items, 0, sizeof(*r_map->items) * (size_t)numedges_dst);
+		}
+
+		free_bvhtree_from_mesh(&treedata);
+	}
+}
+
+#define POLY_UNSET       0
+#define POLY_CENTER_INIT 1
+#define POLY_COMPLETE    2
+
+static void mesh_island_to_astar_graph_edge_process(
+        MeshIslandStore *islands, const int island_index, BLI_AStarGraph *as_graph,
+        MVert *verts, MPoly *polys, MLoop *loops,
+        const int edge_idx, BLI_bitmap *done_edges, MeshElemMap *edge_to_poly_map, const bool is_edge_innercut,
+        int *poly_island_index_map, float (*poly_centers)[3], unsigned char *poly_status)
+{
+	int *poly_island_indices = BLI_array_alloca(poly_island_indices, (size_t)edge_to_poly_map[edge_idx].count);
+	int i, j;
+
+	for (i = 0; i < edge_to_poly_map[edge_idx].count; i++) {
+		const int pidx = edge_to_poly_map[edge_idx].indices[i];
+		MPoly *mp = &polys[pidx];
+		const int pidx_isld = islands ? poly_island_index_map[pidx] : pidx;
+		void *custom_data = is_edge_innercut ? SET_INT_IN_POINTER(edge_idx) : SET_INT_IN_POINTER(-1);
+
+		if (UNLIKELY(islands && (islands->items_to_islands[mp->loopstart] != island_index))) {
+			/* poly not in current island, happens with border edges... */
+			poly_island_indices[i] = -1;
+			continue;
+		}
+
+		if (poly_status[pidx_isld] == POLY_COMPLETE) {
+			poly_island_indices[i] = pidx_isld;
+			continue;
+		}
+
+		if (poly_status[pidx_isld] == POLY_UNSET) {
+			BKE_mesh_calc_poly_center(mp, &loops[mp->loopstart], verts, poly_centers[pidx_isld]);
+			BLI_astar_node_init(as_graph, pidx_isld, poly_centers[pidx_isld]);
+			poly_status[pidx_isld] = POLY_CENTER_INIT;
+		}
+
+		for (j = i; j--;) {
+			float dist_cost;
+			const int pidx_isld_other = poly_island_indices[j];
+
+			if (pidx_isld_other == -1 || poly_status[pidx_isld_other] == POLY_COMPLETE) {
+				/* If the other poly is complete, that link has already been added! */
+				continue;
+			}
+			dist_cost = len_v3v3(poly_centers[pidx_isld_other], poly_centers[pidx_isld]);
+			BLI_astar_node_link_add(as_graph, pidx_isld_other, pidx_isld, dist_cost, custom_data);
+		}
+
+		poly_island_indices[i] = pidx_isld;
+	}
+
+	BLI_BITMAP_ENABLE(done_edges, edge_idx);
+}
+
+static void mesh_island_to_astar_graph(
+        MeshIslandStore *islands, const int island_index,
+        MVert *verts, MeshElemMap *edge_to_poly_map, const int numedges, MLoop *loops, MPoly *polys, const int numpolys,
+        BLI_AStarGraph *r_as_graph)
+{
+	MeshElemMap *island_poly_map = islands ? islands->islands[island_index] : NULL;
+	MeshElemMap *island_einnercut_map = islands ? islands->innercuts[island_index] : NULL;
+
+	int *poly_island_index_map = NULL;
+	BLI_bitmap *done_edges = BLI_BITMAP_NEW(numedges, __func__);
+
+	const int node_num = islands ? island_poly_map->count : numpolys;
+	unsigned char *poly_status = MEM_callocN(sizeof(*poly_status) * (size_t)node_num, __func__);
+	float (*poly_centers)[3];
+
+	int pidx_isld;
+	int i;
+
+	BLI_astar_graph_init(r_as_graph, node_num, NULL);
+	/* poly_centers is owned by graph memarena. */
+	poly_centers = BLI_memarena_calloc(r_as_graph->mem, sizeof(*poly_centers) * (size_t)node_num);
+
+	if (islands) {
+		/* poly_island_index_map is owned by graph memarena. */
+		poly_island_index_map = BLI_memarena_calloc(r_as_graph->mem, sizeof(*poly_island_index_map) * (size_t)numpolys);
+		for (i = island_poly_map->count; i--;) {
+			poly_island_index_map[island_poly_map->indices[i]] = i;
+		}
+
+		r_as_graph->custom_data = poly_island_index_map;
+
+		for (i = island_einnercut_map->count; i--;) {
+			mesh_island_to_astar_graph_edge_process(
+			        islands, island_index, r_as_graph, verts, polys, loops,
+			        island_einnercut_map->indices[i], done_edges, edge_to_poly_map, true,
+			        poly_island_index_map, poly_centers, poly_status);
+		}
+	}
+
+	for (pidx_isld = node_num; pidx_isld--;) {
+		const int pidx = islands ? island_poly_map->indices[pidx_isld] : pidx_isld;
+		MPoly *mp = &polys[pidx];
+		int pl_idx, l_idx;
+
+		if (poly_status[pidx_isld] == POLY_COMPLETE) {
+			continue;
+		}
+
+		for (pl_idx = 0, l_idx = mp->loopstart; pl_idx < mp->totloop; pl_idx++, l_idx++) {
+			MLoop *ml = &loops[l_idx];
+
+			if (BLI_BITMAP_TEST(done_edges, ml->e)) {
+				continue;
+			}
+
+			mesh_island_to_astar_graph_edge_process(
+			        islands, island_index, r_as_graph, verts, polys, loops,
+			        (int)ml->e, done_edges, edge_to_poly_map, false,
+			        poly_island_index_map, poly_centers, poly_status);
+		}
+		poly_status[pidx_isld] = POLY_COMPLETE;
+	}
+
+	MEM_freeN(done_edges);
+	MEM_freeN(poly_status);
+}
+
+#undef POLY_UNSET
+#undef POLY_CENTER_INIT
+#undef POLY_COMPLETE
+
+/* Our 'f_cost' callback func, to find shortest poly-path between two remapped-loops.
+ * Note we do not want to make innercuts 'walls' here, just detect when the shortest path goes by those. */
+static float mesh_remap_calc_loops_astar_f_cost(
+        BLI_AStarGraph *as_graph, BLI_AStarSolution *as_solution, BLI_AStarGNLink *link,
+        const int node_idx_curr, const int node_idx_next, const int node_idx_dst)
+{
+	float *co_next, *co_dest;
+
+	if (link && (GET_INT_FROM_POINTER(link->custom_data) != -1)) {
+		/* An innercut edge... We tag our solution as potentially crossing innercuts.
+		 * Note it might not be the case in the end (AStar will explore around optimal path), but helps
+		 * trimming off some processing later... */
+		if (!GET_INT_FROM_POINTER(as_solution->custom_data)) {
+			as_solution->custom_data = SET_INT_IN_POINTER(true);
+		}
+	}
+
+	/* Our heuristic part of current f_cost is distance from next node to destination one.
+	 * It is guaranteed to be less than (or equal to) actual shortest poly-path between next node and destination one.
+	 */
+	co_next = (float *)as_graph->nodes[node_idx_next].custom_data;
+	co_dest = (float *)as_graph->nodes[node_idx_dst].custom_data;
+	return (link ? (as_solution->g_costs[node_idx_curr] + link->cost) : 0.0f) + len_v3v3(co_next, co_dest);
+}
+
+#define ASTAR_STEPS_MAX 64
+
+
+void BKE_mesh_remap_calc_loops_from_dm(
+        const int mode, const SpaceTransform *space_transform, const float max_dist, const float ray_radius,
+        MVert *verts_dst, const int numverts_dst, MEdge *edges_dst, const int numedges_dst,
+        MLoop *loops_dst, const int numloops_dst, MPoly *polys_dst, const int numpolys_dst,
+        CustomData *ldata_dst, CustomData *pdata_dst, const float split_angle_dst, const bool dirty_nors_dst,
+        DerivedMesh *dm_src,
+        MeshRemapIslandsCalc gen_islands_src, const float islands_precision_src, MeshPairRemap *r_map)
+{
+	const float full_weight = 1.0f;
+	const float max_dist_sq = max_dist * max_dist;
+
+	int i;
+
+	BLI_assert(mode & MREMAP_MODE_LOOP);
+	BLI_assert((islands_precision_src >= 0.0f) && (islands_precision_src <= 1.0f));
+
+	BKE_mesh_remap_init(r_map, numloops_dst);
+
+	if (mode == MREMAP_MODE_TOPOLOGY) {
+		/* In topology mapping, we assume meshes are identical, islands included! */
+		BLI_assert(numloops_dst == dm_src->getNumLoops(dm_src));
+		for (i = 0; i < numloops_dst; i++) {
+			mesh_remap_item_define(r_map, i, FLT_MAX, 0, 1, &i, &full_weight);
+		}
+	}
+	else {
+		BVHTreeFromMesh *treedata = NULL;
+		BVHTreeNearest nearest = {0};
+		BVHTreeRayHit rayhit = {0};
+		int num_trees = 0;
+		float hit_dist;
+
+		const bool use_from_vert = (mode & MREMAP_USE_VERT);
+
+		MeshIslandStore island_store = {0};
+		bool use_islands = false;
+
+		BLI_AStarGraph *as_graphdata = NULL;
+		BLI_AStarSolution as_solution = {0};
+		const int isld_steps_src = islands_precision_src ?
+		                           max_ii((int)(ASTAR_STEPS_MAX * islands_precision_src + 0.499f), 1) : 0;
+
+		float (*poly_nors_src)[3] = NULL;
+		float (*loop_nors_src)[3] = NULL;
+		float (*poly_nors_dst)[3] = NULL;
+		float (*loop_nors_dst)[3] = NULL;
+
+		MeshElemMap *vert_to_loop_map_src = NULL;
+		int *vert_to_loop_map_src_buff = NULL;
+		MeshElemMap *vert_to_poly_map_src = NULL;
+		int *vert_to_poly_map_src_buff = NULL;
+		MeshElemMap *edge_to_poly_map_src = NULL;
+		int *edge_to_poly_map_src_buff = NULL;
+		MeshElemMap *poly_to_tessface_map_src = NULL;
+		int *poly_to_tessface_map_src_buff = NULL;
+
+		/* Unlike above, those are one-to-one mappings, simpler! */
+		int *loop_to_poly_map_src = NULL;
+		int *tessface_to_poly_map_src = NULL;
+
+		bool verts_allocated_src;
+		MVert *verts_src = DM_get_vert_array(dm_src, &verts_allocated_src);
+		const int num_verts_src = dm_src->getNumVerts(dm_src);
+		float (*vcos_src)[3] = NULL;
+		bool edges_allocated_src;
+		MEdge *edges_src = DM_get_edge_array(dm_src, &edges_allocated_src);
+		const int num_edges_src = dm_src->getNumEdges(dm_src);
+		bool loops_allocated_src;
+		MLoop *loops_src = DM_get_loop_array(dm_src, &loops_allocated_src);
+		const int num_loops_src = dm_src->getNumLoops(dm_src);
+		bool polys_allocated_src;
+		MPoly *polys_src = DM_get_poly_array(dm_src, &polys_allocated_src);
+		const int num_polys_src = dm_src->getNumPolys(dm_src);
+		bool faces_allocated_src = false;
+		MFace *faces_src = NULL;
+		int num_faces_src = 0;
+
+		size_t buff_size_interp = MREMAP_DEFAULT_BUFSIZE;
+		float (*vcos_interp)[3] = NULL;
+		int *indices_interp = NULL;
+		float *weights_interp = NULL;
+
+		int tindex, pidx_dst, lidx_dst, plidx_dst, pidx_src, lidx_src, plidx_src;
+
+		IslandResult **islands_res;
+		size_t islands_res_buff_size = MREMAP_DEFAULT_BUFSIZE;
+
+		const float bvh_epsilon = (mode & MREMAP_USE_NORPROJ) ? MREMAP_RAYCAST_APPROXIMATE_BVHEPSILON(ray_radius) : 0.0f;
+
+		if (!use_from_vert) {
+			vcos_src = MEM_mallocN(sizeof(*vcos_src) * (size_t)num_verts_src, __func__);
+			dm_src->getVertCos(dm_src, vcos_src);
+
+			vcos_interp = MEM_mallocN(sizeof(*vcos_interp) * buff_size_interp, __func__);
+			indices_interp = MEM_mallocN(sizeof(*indices_interp) * buff_size_interp, __func__);
+			weights_interp = MEM_mallocN(sizeof(*weights_interp) * buff_size_interp, __func__);
+		}
+
+		{
+			const bool need_lnors_src = (mode & MREMAP_USE_LOOP) && (mode & MREMAP_USE_NORMAL);
+			const bool need_lnors_dst = need_lnors_src || (mode & MREMAP_USE_NORPROJ);
+			const bool need_pnors_src = need_lnors_src || ((mode & MREMAP_USE_POLY) && (mode & MREMAP_USE_NORMAL));
+			const bool need_pnors_dst = need_lnors_dst || need_pnors_src;
+
+			if (need_pnors_dst) {
+				/* Cache poly nors into a temp CDLayer. */
+				poly_nors_dst = CustomData_get_layer(pdata_dst, CD_NORMAL);
+				if (!poly_nors_dst) {
+					poly_nors_dst = CustomData_add_layer(pdata_dst, CD_NORMAL, CD_CALLOC, NULL, numpolys_dst);
+					CustomData_set_layer_flag(pdata_dst, CD_NORMAL, CD_FLAG_TEMPORARY);
+				}
+				if (dirty_nors_dst) {
+					BKE_mesh_calc_normals_poly(verts_dst, numverts_dst, loops_dst, polys_dst,
+					                           numloops_dst, numpolys_dst, poly_nors_dst, true);
+				}
+			}
+			if (need_lnors_dst) {
+				/* Cache poly nors into a temp CDLayer. */
+				loop_nors_dst = CustomData_get_layer(ldata_dst, CD_NORMAL);
+				if (!loop_nors_dst) {
+					loop_nors_dst = CustomData_add_layer(ldata_dst, CD_NORMAL, CD_CALLOC, NULL, numloops_dst);
+					CustomData_set_layer_flag(ldata_dst, CD_NORMAL, CD_FLAG_TEMPORARY);
+				}
+				if (dirty_nors_dst) {
+					BKE_mesh_normals_loop_split(verts_dst, numverts_dst, edges_dst, numedges_dst,
+					                            loops_dst, loop_nors_dst, numloops_dst,
+					                            polys_dst, poly_nors_dst, numpolys_dst, split_angle_dst);
+				}
+			}
+			if (need_pnors_src || need_lnors_src) {
+				/* Simpler for now, calcNormals never stores pnors :( */
+				dm_src->calcLoopNormals(dm_src, /* TODO */ (float)M_PI);
+
+				if (need_pnors_src) {
+					poly_nors_src = dm_src->getPolyDataArray(dm_src, CD_NORMAL);
+				}
+				if (need_lnors_src) {
+					loop_nors_src = dm_src->getLoopDataArray(dm_src, CD_NORMAL);
+				}
+			}
+		}
+
+		if (use_from_vert) {
+			BKE_mesh_vert_loop_map_create(&vert_to_loop_map_src, &vert_to_loop_map_src_buff,
+			                              polys_src, loops_src, num_verts_src, num_polys_src, num_loops_src);
+			if (mode & MREMAP_USE_POLY) {
+				BKE_mesh_vert_poly_map_create(&vert_to_poly_map_src, &vert_to_poly_map_src_buff,
+				                              polys_src, loops_src, num_verts_src, num_polys_src, num_loops_src);
+			}
+		}
+
+		/* Needed for islands (or plain mesh) to AStar graph conversion. */
+		BKE_mesh_edge_poly_map_create(&edge_to_poly_map_src, &edge_to_poly_map_src_buff,
+		                              edges_src, num_edges_src, polys_src, num_polys_src, loops_src, num_loops_src);
+		if (use_from_vert) {
+			loop_to_poly_map_src = MEM_mallocN(sizeof(*loop_to_poly_map_src) * (size_t)num_loops_src, __func__);
+			for (pidx_src = 0; pidx_src < num_polys_src; pidx_src++) {
+				MPoly *mp = &polys_src[pidx_src];
+				for (plidx_src = 0, lidx_src = mp->loopstart; plidx_src < mp->totloop; plidx_src++, lidx_src++) {
+					loop_to_poly_map_src[lidx_src] = pidx_src;
+				}
+			}
+		}
+
+		/* Island makes things slightly more complex here.
+		 * Basically, we:
+		 *     * Make one treedata for each island's elements.
+		 *     * Check all loops of a same dest poly against all treedata.
+		 *     * Choose the island's elements giving the best results.
+		 */
+
+		/* First, generate the islands, if possible. */
+		if (gen_islands_src) {
+			use_islands = gen_islands_src(
+			        verts_src, num_verts_src,
+			        edges_src, num_edges_src,
+			        polys_src, num_polys_src,
+			        loops_src, num_loops_src,
+			        &island_store);
+
+			num_trees = use_islands ? island_store.islands_num : 1;
+			treedata = MEM_callocN(sizeof(*treedata) * (size_t)num_trees, __func__);
+			if (isld_steps_src) {
+				as_graphdata = MEM_callocN(sizeof(*as_graphdata) * (size_t)num_trees, __func__);
+			}
+
+			if (use_islands) {
+				/* We expect our islands to contain poly indices, with edge indices of 'inner cuts',
+				 * and a mapping loops -> islands indices.
+				 * This implies all loops of a same poly are in the same island. */
+				BLI_assert((island_store.item_type == MISLAND_TYPE_LOOP) &&
+				           (island_store.island_type == MISLAND_TYPE_POLY) &&
+				           (island_store.innercut_type == MISLAND_TYPE_EDGE));
+			}
+		}
+		else {
+			num_trees = 1;
+			treedata = MEM_callocN(sizeof(*treedata), __func__);
+			if (isld_steps_src) {
+				as_graphdata = MEM_callocN(sizeof(*as_graphdata), __func__);
+			}
+		}
+
+		/* Build our AStar graphs. */
+		if (isld_steps_src) {
+			for (tindex = 0; tindex < num_trees; tindex++) {
+				mesh_island_to_astar_graph(
+				        use_islands ? &island_store : NULL, tindex, verts_src, edge_to_poly_map_src, num_edges_src,
+				        loops_src, polys_src, num_polys_src, &as_graphdata[tindex]);
+			}
+		}
+
+		/* Build our BVHtrees, either from verts or tessfaces. */
+		if (use_from_vert) {
+			if (use_islands) {
+				BLI_bitmap *verts_active = BLI_BITMAP_NEW((size_t)num_verts_src, __func__);
+
+				for (tindex = 0; tindex < num_trees; tindex++) {
+					MeshElemMap *isld = island_store.islands[tindex];
+					int num_verts_active = 0;
+					BLI_BITMAP_SET_ALL(verts_active, false, (size_t)num_verts_src);
+					for (i = 0; i < isld->count; i++) {
+						MPoly *mp_src = &polys_src[isld->indices[i]];
+						for (lidx_src = mp_src->loopstart; lidx_src < mp_src->loopstart + mp_src->totloop; lidx_src++) {
+							BLI_BITMAP_ENABLE(verts_active, loops_src[lidx_src].v);
+							num_verts_active++;
+						}
+					}
+					/* verts 'ownership' is transfered to treedata here, which will handle its freeing. */
+					bvhtree_from_mesh_verts_ex(&treedata[tindex], verts_src, num_verts_src, verts_allocated_src,
+					                           verts_active, num_verts_active, bvh_epsilon, 2, 6);
+					if (verts_allocated_src) {
+						verts_allocated_src = false;  /* Only 'give' our verts once, to first tree! */
+					}
+				}
+
+				MEM_freeN(verts_active);
+			}
+			else {
+				BLI_assert(num_trees == 1);
+				bvhtree_from_mesh_verts(&treedata[0], dm_src, bvh_epsilon, 2, 6);
+			}
+		}
+		else {  /* We use polygons. */
+			if (use_islands) {
+				/* bvhtree here uses tesselated faces... */
+				const unsigned int dirty_tess_flag = dm_src->dirty & DM_DIRTY_TESS_CDLAYERS;
+				BLI_bitmap *faces_active;
+
+				/* We do not care about tessellated data here, only geometry itself is important. */
+				if (dirty_tess_flag) {
+					dm_src->dirty &= ~dirty_tess_flag;
+				}
+				DM_ensure_tessface(dm_src);
+				if (dirty_tess_flag) {
+					dm_src->dirty |= dirty_tess_flag;
+				}
+				faces_src = DM_get_tessface_array(dm_src, &faces_allocated_src);
+				num_faces_src = dm_src->getNumTessFaces(dm_src);
+				tessface_to_poly_map_src = dm_src->getTessFaceDataArray(dm_src, CD_ORIGINDEX);
+				faces_active = BLI_BITMAP_NEW((size_t)num_faces_src, __func__);
+
+				for (tindex = 0; tindex < num_trees; tindex++) {
+					int num_faces_active = 0;
+					BLI_BITMAP_SET_ALL(faces_active, false, (size_t)num_faces_src);
+					for (i = 0; i < num_faces_src; i++) {
+						MPoly *mp_src = &polys_src[tessface_to_poly_map_src[i]];
+						if (island_store.items_to_islands[mp_src->loopstart] == tindex) {
+							BLI_BITMAP_ENABLE(faces_active, i);
+							num_faces_active++;
+						}
+					}
+					/* verts 'ownership' is transfered to treedata here, which will handle its freeing. */
+					bvhtree_from_mesh_faces_ex(
+					        &treedata[tindex], verts_src, verts_allocated_src,
+					        faces_src, num_faces_src, faces_allocated_src,
+					        faces_active, num_faces_active, bvh_epsilon, 2, 6);
+					if (verts_allocated_src) {
+						verts_allocated_src = false;  /* Only 'give' our verts once, to first tree! */
+					}
+					if (faces_allocated_src) {
+						faces_allocated_src = false;  /* Only 'give' our faces once, to first tree! */
+					}
+				}
+
+				MEM_freeN(faces_active);
+			}
+			else {
+				BLI_assert(num_trees == 1);
+				bvhtree_from_mesh_faces(&treedata[0], dm_src, bvh_epsilon, 2, 6);
+				tessface_to_poly_map_src = dm_src->getTessFaceDataArray(dm_src, CD_ORIGINDEX);
+			}
+		}
+
+		/* And check each dest poly! */
+		islands_res = MEM_mallocN(sizeof(*islands_res) * (size_t)num_trees, __func__);
+		for (tindex = 0; tindex < num_trees; tindex++) {
+			islands_res[tindex] = MEM_mallocN(sizeof(**islands_res) * islands_res_buff_size, __func__);
+		}
+
+		for (pidx_dst = 0; pidx_dst < numpolys_dst; pidx_dst++) {
+			MPoly *mp_dst = &polys_dst[pidx_dst];
+			float (*pnor_dst)[3] = &poly_nors_dst[pidx_dst];
+
+			if ((size_t)mp_dst->totloop > islands_res_buff_size) {
+				islands_res_buff_size = (size_t)mp_dst->totloop;
+				for (tindex = 0; tindex < num_trees; tindex++) {
+					islands_res[tindex] = MEM_reallocN(islands_res[tindex], sizeof(**islands_res) * islands_res_buff_size);
+				}
+			}
+
+			for (tindex = 0; tindex < num_trees; tindex++) {
+				BVHTreeFromMesh *tdata = &treedata[tindex];
+				MLoop *ml_dst = &loops_dst[mp_dst->loopstart];
+
+				for (plidx_dst = 0; plidx_dst < mp_dst->totloop; plidx_dst++, ml_dst++) {
+					if (use_from_vert) {
+						float tmp_co[3];
+						MeshElemMap *vert_to_refelem_map_src = NULL;
+
+						copy_v3_v3(tmp_co, verts_dst[ml_dst->v].co);
+						nearest.index = -1;
+
+						if (mesh_remap_bvhtree_query_nearest(
+						        tdata, &nearest, space_transform,
+						        tmp_co, max_dist_sq, &hit_dist))
+						{
+							float (*nor_dst)[3];
+							float (*nors_src)[3];
+							float best_nor_dot = -2.0f;
+							int best_index_src = -1;
+
+							if (mode == MREMAP_MODE_LOOP_NEAREST_LOOPNOR) {
+								nor_dst = &loop_nors_dst[plidx_dst + mp_dst->loopstart];
+								nors_src = loop_nors_src;
+								vert_to_refelem_map_src = vert_to_loop_map_src;
+							}
+							else {  /* if (mode == MREMAP_MODE_LOOP_NEAREST_POLYNOR) { */
+								nor_dst = pnor_dst;
+								nors_src = poly_nors_src;
+								vert_to_refelem_map_src = vert_to_poly_map_src;
+							}
+
+							for (i = vert_to_refelem_map_src[nearest.index].count; i--;) {
+								const int index_src = vert_to_refelem_map_src[nearest.index].indices[i];
+								const float dot = dot_v3v3(nors_src[index_src], *nor_dst);
+								if (dot > best_nor_dot) {
+									best_nor_dot = dot;
+									best_index_src = index_src;
+								}
+							}
+							if (mode == MREMAP_MODE_LOOP_NEAREST_POLYNOR) {
+								/* Our best_index_src is a poly one for now!
+								 * Have to find its loop matching our closest vertex. */
+								MPoly *mp_src = &polys_src[best_index_src];
+								MLoop *ml_src = &loops_src[mp_src->loopstart];
+								for (plidx_src = 0; plidx_src < mp_src->totloop; plidx_src++, ml_src++) {
+									if ((int)ml_src->v == nearest.index) {
+										best_index_src = plidx_src + mp_src->loopstart;
+										break;
+									}
+								}
+							}
+							islands_res[tindex][plidx_dst].factor = hit_dist ? (1.0f / (hit_dist * best_nor_dot)) : 1e18f;
+							islands_res[tindex][plidx_dst].hit_dist = hit_dist;
+							islands_res[tindex][plidx_dst].index_src = best_index_src;
+						}
+						else {
+							/* No source for this dest loop! */
+							islands_res[tindex][plidx_dst].factor = 0.0f;
+							islands_res[tindex][plidx_dst].hit_dist = FLT_MAX;
+							islands_res[tindex][plidx_dst].index_src = -1;
+						}
+					}
+					else if (mode & MREMAP_USE_NORPROJ) {
+						float tmp_co[3], tmp_no[3];
+
+						int n = (ray_radius > 0.0f) ? MREMAP_RAYCAST_APPROXIMATE_NR : 1;
+						float w = 1.0f;
+
+						copy_v3_v3(tmp_co, verts_dst[ml_dst->v].co);
+						copy_v3_v3(tmp_no, loop_nors_dst[plidx_dst + mp_dst->loopstart]);
+
+						/* We do our transform here, since we may do several raycast/nearest queries. */
+						if (space_transform) {
+							BLI_space_transform_apply(space_transform, tmp_co);
+							BLI_space_transform_apply_normal(space_transform, tmp_no);
+						}
+
+						while (n--) {
+							/* Note we handle dest to src space conversion ourself, here! */
+							if (mesh_remap_bvhtree_query_raycast(
+							        tdata, &rayhit, NULL,
+							        tmp_co, tmp_no, ray_radius / w, max_dist, &hit_dist))
+							{
+								islands_res[tindex][plidx_dst].factor = (hit_dist ? (1.0f / hit_dist) : 1e18f) * w;
+								islands_res[tindex][plidx_dst].hit_dist = hit_dist;
+								islands_res[tindex][plidx_dst].index_src = tessface_to_poly_map_src[rayhit.index];
+								copy_v3_v3(islands_res[tindex][plidx_dst].hit_point, rayhit.co);
+								break;
+							}
+							/* Next iteration will get bigger radius but smaller weight! */
+							w /= MREMAP_RAYCAST_APPROXIMATE_FAC;
+						}
+						if (n == -1) {
+							/* Fallback to 'nearest' hit here, loops usually comes in 'face group', not good to
+							 * have only part of one dest face's loops to map to source.
+							 * Note that since we give this a null weight, if whole weight for a given face
+							 * is null, it means none of its loop mapped to this source island, hence we can skip it
+							 * later.
+							 */
+							copy_v3_v3(tmp_co, verts_dst[ml_dst->v].co);
+							nearest.index = -1;
+
+							/* In any case, this fallback nearest hit should have no weight at all
+							 * in 'best island' decision! */
+							islands_res[tindex][plidx_dst].factor = 0.0f;
+
+							/* Note we handle dest to src space conversion ourself, here! */
+							if (mesh_remap_bvhtree_query_nearest(
+							        tdata, &nearest, NULL,
+							        tmp_co, max_dist_sq, &hit_dist))
+							{
+								islands_res[tindex][plidx_dst].hit_dist = hit_dist;
+								islands_res[tindex][plidx_dst].index_src = tessface_to_poly_map_src[nearest.index];
+								copy_v3_v3(islands_res[tindex][plidx_dst].hit_point, nearest.co);
+							}
+							else {
+								/* No source for this dest loop! */
+								islands_res[tindex][plidx_dst].hit_dist = FLT_MAX;
+								islands_res[tindex][plidx_dst].index_src = -1;
+							}
+						}
+					}
+					else {  /* Nearest poly either to use all its loops/verts or just closest one. */
+						float tmp_co[3];
+
+						copy_v3_v3(tmp_co, verts_dst[ml_dst->v].co);
+						nearest.index = -1;
+
+						if (mesh_remap_bvhtree_query_nearest(
+						        tdata, &nearest, space_transform,
+						        tmp_co, max_dist_sq, &hit_dist))
+						{
+							islands_res[tindex][plidx_dst].factor = hit_dist ? (1.0f / hit_dist) : 1e18f;
+							islands_res[tindex][plidx_dst].hit_dist = hit_dist;
+							islands_res[tindex][plidx_dst].index_src = tessface_to_poly_map_src[nearest.index];
+							copy_v3_v3(islands_res[tindex][plidx_dst].hit_point, nearest.co);
+						}
+						else {
+							/* No source for this dest loop! */
+							islands_res[tindex][plidx_dst].factor = 0.0f;
+							islands_res[tindex][plidx_dst].hit_dist = FLT_MAX;
+							islands_res[tindex][plidx_dst].index_src = -1;
+						}
+					}
+				}
+			}
+
+			/* And now, find best island to use! */
+			/* We have to first select the 'best source island' for given dst poly and its loops.
+			 * Then, we have to check that poly does not 'spread' across some island's limits
+			 * (like inner seams for UVs, etc.).
+			 * Note we only still partially support that kind of situation here, i.e. polys spreading over actual cracks
+			 * (like a narrow space without faces on src, splitting a 'tube-like' geometry). That kind of situation
+			 * should be relatively rare, though.
+			 */
+			/* XXX This block in itself is big and complex enough to be a separate function but... it uses a bunch
+			 *     of locale vars. Not worth sending all that through parameters (for now at least). */
+			{
+				BLI_AStarGraph *as_graph = NULL;
+				int *poly_island_index_map = NULL;
+				int pidx_src_prev = -1;
+
+				MeshElemMap *best_island = NULL;
+				float best_island_fac = 0.0f;
+				int best_island_index = -1;
+
+				for (tindex = 0; tindex < num_trees; tindex++) {
+					float island_fac = 0.0f;
+
+					for (plidx_dst = 0; plidx_dst < mp_dst->totloop; plidx_dst++) {
+						island_fac += islands_res[tindex][plidx_dst].factor;
+					}
+					island_fac /= (float)mp_dst->totloop;
+
+					if (island_fac > best_island_fac) {
+						best_island_fac = island_fac;
+						best_island_index = tindex;
+					}
+				}
+
+				if (best_island_index != -1 && isld_steps_src) {
+					best_island = island_store.islands[best_island_index];
+					as_graph = &as_graphdata[best_island_index];
+					poly_island_index_map = (int *)as_graph->custom_data;
+					BLI_astar_solution_init(as_graph, &as_solution, false);
+				}
+
+				for (plidx_dst = 0; plidx_dst < mp_dst->totloop; plidx_dst++) {
+					IslandResult *isld_res;
+					lidx_dst = plidx_dst + mp_dst->loopstart;
+
+					if (best_island_index == -1) {
+						/* No source for any loops of our dest poly in any source islands. */
+						BKE_mesh_remap_item_define_invalid(r_map, lidx_dst);
+						continue;
+					}
+
+					as_solution.custom_data = SET_INT_IN_POINTER(false);
+
+					isld_res = &islands_res[best_island_index][plidx_dst];
+					if (use_from_vert) {
+						/* Indices stored in islands_res are those of loops, one per dest loop. */
+						lidx_src = isld_res->index_src;
+						if (lidx_src >= 0) {
+							pidx_src = loop_to_poly_map_src[lidx_src];
+							/* If prev and curr poly are the same, no need to do anything more!!! */
+							if (!ELEM(pidx_src_prev, -1, pidx_src) && isld_steps_src) {
+								BLI_astar_graph_solve(
+								        as_graph, poly_island_index_map[pidx_src_prev], poly_island_index_map[pidx_src],
+								        mesh_remap_calc_loops_astar_f_cost, &as_solution, isld_steps_src);
+								if (GET_INT_FROM_POINTER(as_solution.custom_data) && (as_solution.steps > 0)) {
+									/* Find first 'cutting edge' on path, and bring back lidx_src on poly just
+									 * before that edge.
+									 * Note we could try to be much smarter (like e.g. storing a whole poly's indices,
+									 * and making decision (on which side of cutting edge(s!) to be) on the end,
+									 * but this is one more level of complexity, better to first see if
+									 * simple solution works!
+									 */
+									int pidx_isld_src = poly_island_index_map[pidx_src];
+									int last_valid_pidx_isld_src = -1;
+									/* Note we go backward here, from dest to src poly. */
+									for (i = as_solution.steps - 1; i--;) {
+										BLI_AStarGNLink *as_link = as_solution.prev_links[pidx_isld_src];
+										const int eidx = GET_INT_FROM_POINTER(as_link->custom_data);
+										pidx_isld_src = as_solution.prev_nodes[pidx_isld_src];
+										BLI_assert(pidx_isld_src != -1);
+										if (eidx != -1) {
+											/* we are 'crossing' a cutting edge. */
+											last_valid_pidx_isld_src = pidx_isld_src;
+										}
+									}
+									if (last_valid_pidx_isld_src != -1) {
+										/* Find a new valid loop in that new poly (nearest one for now).
+										 * Note we could be much more subtle here, again that's for later... */
+										MPoly *mp_src;
+										MLoop *ml_src, *ml_dst = &loops_dst[lidx_dst];
+										int j;
+										float best_dist_sq = FLT_MAX;
+										float tmp_co[3];
+
+										copy_v3_v3(tmp_co, verts_dst[ml_dst->v].co);
+
+										/* We do our transform here, since we may do several raycast/nearest queries. */
+										if (space_transform) {
+											BLI_space_transform_apply(space_transform, tmp_co);
+										}
+
+										pidx_src = best_island->indices[last_valid_pidx_isld_src];
+										mp_src = &polys_src[pidx_src];
+										ml_src = &loops_src[mp_src->loopstart];
+										for (j = 0; j < mp_src->totloop; j++, ml_src++) {
+											const float dist_sq = len_squared_v3v3(verts_src[ml_src->v].co, tmp_co);
+											if (dist_sq < best_dist_sq) {
+												best_dist_sq = dist_sq;
+												lidx_src = mp_src->loopstart + j;
+											}
+										}
+									}
+								}
+							}
+							mesh_remap_item_define(
+							        r_map, lidx_dst, isld_res->hit_dist,
+							        best_island_index, 1, &lidx_src, &full_weight);
+							pidx_src_prev = pidx_src;
+						}
+						else {
+							/* No source for this loop in this island. */
+							/* TODO: would probably be better to get a source at all cost in best island anyway? */
+							mesh_remap_item_define(
+							        r_map, lidx_dst, FLT_MAX,
+							        best_island_index, 0, NULL, NULL);
+						}
+					}
+					else {
+						/* Else, we use source poly, indices stored in islands_res are those of polygons. */
+						pidx_src = isld_res->index_src;
+						if (pidx_src >= 0) {
+							MPoly *mp_src = &polys_src[pidx_src];
+							float *hit_co = isld_res->hit_point;
+							int best_loop_index_src;
+
+							/* If prev and curr poly are the same, no need to do anything more!!! */
+							if (!ELEM(pidx_src_prev, -1, pidx_src) && isld_steps_src) {
+								BLI_astar_graph_solve(
+								        as_graph, poly_island_index_map[pidx_src_prev], poly_island_index_map[pidx_src],
+								        mesh_remap_calc_loops_astar_f_cost, &as_solution, isld_steps_src);
+								if (GET_INT_FROM_POINTER(as_solution.custom_data) && (as_solution.steps > 0)) {
+									/* Find first 'cutting edge' on path, and bring back lidx_src on poly just
+									 * before that edge.
+									 * Note we could try to be much smarter (like e.g. storing a whole poly's indices,
+									 * and making decision (one which side of cutting edge(s!) to be on the end,
+									 * but this is one more level of complexity, better to first see if
+									 * simple solution works!
+									 */
+									int pidx_isld_src = poly_island_index_map[pidx_src];
+									int last_valid_pidx_isld_src = -1;
+									/* Note we go backward here, from dest to src poly. */
+									for (i = as_solution.steps - 1; i--;) {
+										BLI_AStarGNLink *as_link = as_solution.prev_links[pidx_isld_src];
+										int eidx = GET_INT_FROM_POINTER(as_link->custom_data);
+
+										pidx_isld_src = as_solution.prev_nodes[pidx_isld_src];
+										BLI_assert(pidx_isld_src != -1);
+										if (eidx != -1) {
+											/* we are 'crossing' a cutting edge. */
+											last_valid_pidx_isld_src = pidx_isld_src;
+										}
+									}
+									if (last_valid_pidx_isld_src != -1) {
+										/* Find a new valid loop in that new poly (nearest point on poly for now).
+										 * Note we could be much more subtle here, again that's for later... */
+										MLoop *ml_dst = &loops_dst[lidx_dst];
+										float best_dist_sq = FLT_MAX;
+										float tmp_co[3];
+										int j;
+
+										copy_v3_v3(tmp_co, verts_dst[ml_dst->v].co);
+
+										/* We do our transform here, since we may do several raycast/nearest queries. */
+										if (space_transform) {
+											BLI_space_transform_apply(space_transform, tmp_co);
+										}
+
+										pidx_src = best_island->indices[last_valid_pidx_isld_src];
+										mp_src = &polys_src[pidx_src];
+
+										/* Create that one on demand. */
+										if (poly_to_tessface_map_src == NULL) {
+											BKE_mesh_origindex_map_create(
+											        &poly_to_tessface_map_src, &poly_to_tessface_map_src_buff,
+											        num_faces_src, tessface_to_poly_map_src, num_polys_src);
+										}
+
+										for (j = poly_to_tessface_map_src[pidx_src].count; j--;) {
+											float h[3];
+											MFace *mf = &faces_src[poly_to_tessface_map_src[pidx_src].indices[j]];
+											float dist_sq;
+
+											closest_on_tri_to_point_v3(
+											        h, tmp_co,
+											        vcos_src[mf->v1], vcos_src[mf->v2], vcos_src[mf->v3]);
+											dist_sq = len_squared_v3v3(tmp_co, h);
+											if (dist_sq < best_dist_sq) {
+												copy_v3_v3(hit_co, h);
+												best_dist_sq = dist_sq;
+											}
+											if (mf->v4) {
+												closest_on_tri_to_point_v3(
+												        h, tmp_co,
+												        vcos_src[mf->v1], vcos_src[mf->v3], vcos_src[mf->v4]);
+												dist_sq = len_squared_v3v3(tmp_co, h);
+												if (dist_sq < best_dist_sq) {
+													copy_v3_v3(hit_co, h);
+													best_dist_sq = dist_sq;
+												}
+											}
+										}
+									}
+								}
+							}
+
+							if (mode == MREMAP_MODE_LOOP_POLY_NEAREST) {
+								mesh_remap_interp_poly_data_get(
+								        mp_src, loops_src, (const float (*)[3])vcos_src, hit_co,
+								        &buff_size_interp, &vcos_interp, true, &indices_interp,
+								        &weights_interp, false, &best_loop_index_src);
+
+								mesh_remap_item_define(
+								        r_map, lidx_dst, isld_res->hit_dist,
+								        best_island_index, 1, &best_loop_index_src, &full_weight);
+							}
+							else {
+								const int sources_num = mesh_remap_interp_poly_data_get(
+								        mp_src, loops_src, (const float (*)[3])vcos_src, hit_co,
+								        &buff_size_interp, &vcos_interp, true, &indices_interp,
+								        &weights_interp, true, NULL);
+
+								mesh_remap_item_define(
+								        r_map, lidx_dst,
+								        isld_res->hit_dist, best_island_index,
+								        sources_num, indices_interp, weights_interp);
+							}
+
+							pidx_src_prev = pidx_src;
+						}
+						else {
+							/* No source for this loop in this island. */
+							/* TODO: would probably be better to get a source at all cost in best island anyway? */
+							mesh_remap_item_define(r_map, lidx_dst, FLT_MAX, best_island_index, 0, NULL, NULL);
+						}
+					}
+				}
+
+				BLI_astar_solution_clear(&as_solution);
+			}
+		}
+
+		for (tindex = 0; tindex < num_trees; tindex++) {
+			MEM_freeN(islands_res[tindex]);
+			free_bvhtree_from_mesh(&treedata[tindex]);
+			if (isld_steps_src) {
+				BLI_astar_graph_free(&as_graphdata[tindex]);
+			}
+		}
+		MEM_freeN(islands_res);
+		BKE_mesh_loop_islands_free(&island_store);
+		MEM_freeN(treedata);
+		if (isld_steps_src) {
+			MEM_freeN(as_graphdata);
+			BLI_astar_solution_free(&as_solution);
+		}
+
+		if (verts_allocated_src) {
+			MEM_freeN(verts_src);
+		}
+		if (vcos_src) {
+			MEM_freeN(vcos_src);
+		}
+		if (edges_allocated_src) {
+			MEM_freeN(edges_src);
+		}
+		if (loops_allocated_src) {
+			MEM_freeN(loops_src);
+		}
+		if (polys_allocated_src) {
+			MEM_freeN(polys_src);
+		}
+		if (faces_allocated_src) {
+			MEM_freeN(faces_src);
+		}
+		if (vert_to_loop_map_src_buff) {
+			MEM_freeN(vert_to_loop_map_src_buff);
+		}
+		if (vert_to_poly_map_src_buff) {
+			MEM_freeN(vert_to_poly_map_src_buff);
+		}
+		if (edge_to_poly_map_src_buff) {
+			MEM_freeN(edge_to_poly_map_src_buff);
+		}
+		if (poly_to_tessface_map_src_buff) {
+			MEM_freeN(poly_to_tessface_map_src_buff);
+		}
+		if (loop_to_poly_map_src) {
+			MEM_freeN(loop_to_poly_map_src);
+		}
+		if (vcos_interp) {
+			MEM_freeN(vcos_interp);
+		}
+		if (indices_interp) {
+			MEM_freeN(indices_interp);
+		}
+		if (weights_interp) {
+			MEM_freeN(weights_interp);
+		}
+	}
+}
+
+void BKE_mesh_remap_calc_polys_from_dm(
+        const int mode, const SpaceTransform *space_transform, const float max_dist, const float ray_radius,
+        MVert *verts_dst, const int numverts_dst, MLoop *loops_dst, const int numloops_dst,
+        MPoly *polys_dst, const int numpolys_dst, CustomData *pdata_dst, const bool dirty_nors_dst,
+        DerivedMesh *dm_src, MeshPairRemap *r_map)
+{
+	const float full_weight = 1.0f;
+	const float max_dist_sq = max_dist * max_dist;
+	float (*poly_nors_dst)[3] = NULL;
+	int i;
+
+	BLI_assert(mode & MREMAP_MODE_POLY);
+
+	if (mode & (MREMAP_USE_NORMAL | MREMAP_USE_NORPROJ)) {
+		/* Cache poly nors into a temp CDLayer. */
+		poly_nors_dst = CustomData_get_layer(pdata_dst, CD_NORMAL);
+		if (!poly_nors_dst) {
+			poly_nors_dst = CustomData_add_layer(pdata_dst, CD_NORMAL, CD_CALLOC, NULL, numpolys_dst);
+			CustomData_set_layer_flag(pdata_dst, CD_NORMAL, CD_FLAG_TEMPORARY);
+		}
+		if (dirty_nors_dst) {
+			BKE_mesh_calc_normals_poly(verts_dst, numverts_dst, loops_dst, polys_dst, numloops_dst, numpolys_dst,
+			                           poly_nors_dst, true);
+		}
+	}
+
+	BKE_mesh_remap_init(r_map, numpolys_dst);
+
+	if (mode == MREMAP_MODE_TOPOLOGY) {
+		BLI_assert(numpolys_dst == dm_src->getNumPolys(dm_src));
+		for (i = 0; i < numpolys_dst; i++) {
+			mesh_remap_item_define(r_map, i, FLT_MAX, 0, 1, &i, &full_weight);
+		}
+	}
+	else {
+		BVHTreeFromMesh treedata = {NULL};
+		BVHTreeNearest nearest = {0};
+		BVHTreeRayHit rayhit = {0};
+		float hit_dist;
+
+		int *tessface_to_poly_map_src;
+
+		bvhtree_from_mesh_faces(
+		        &treedata, dm_src,
+		        (mode & MREMAP_USE_NORPROJ) ? MREMAP_RAYCAST_APPROXIMATE_BVHEPSILON(ray_radius) : 0.0f,
+		        2, 6);
+		/* bvhtree here uses tesselated faces... */
+		tessface_to_poly_map_src = dm_src->getTessFaceDataArray(dm_src, CD_ORIGINDEX);
+
+		if (mode == MREMAP_MODE_POLY_NEAREST) {
+			nearest.index = -1;
+
+			for (i = 0; i < numpolys_dst; i++) {
+				MPoly *mp = &polys_dst[i];
+				float tmp_co[3];
+
+				BKE_mesh_calc_poly_center(mp, &loops_dst[mp->loopstart], verts_dst, tmp_co);
+
+				if (mesh_remap_bvhtree_query_nearest(
+				        &treedata, &nearest, space_transform,
+				        tmp_co, max_dist_sq, &hit_dist))
+				{
+					mesh_remap_item_define(
+					        r_map, i, hit_dist, 0,
+					        1, &tessface_to_poly_map_src[nearest.index], &full_weight);
+				}
+				else {
+					/* No source for this dest poly! */
+					BKE_mesh_remap_item_define_invalid(r_map, i);
+				}
+			}
+		}
+		else if (mode == MREMAP_MODE_POLY_NOR) {
+			BLI_assert(poly_nors_dst);
+
+			for (i = 0; i < numpolys_dst; i++) {
+				MPoly *mp = &polys_dst[i];
+				float tmp_co[3], tmp_no[3];
+
+				BKE_mesh_calc_poly_center(mp, &loops_dst[mp->loopstart], verts_dst, tmp_co);
+				copy_v3_v3(tmp_no, poly_nors_dst[i]);
+
+				if (mesh_remap_bvhtree_query_raycast(
+				        &treedata, &rayhit, space_transform,
+				        tmp_co, tmp_no, ray_radius, max_dist, &hit_dist))
+				{
+					mesh_remap_item_define(
+					        r_map, i, hit_dist, 0,
+					        1, &tessface_to_poly_map_src[rayhit.index], &full_weight);
+				}
+				else {
+					/* No source for this dest poly! */
+					BKE_mesh_remap_item_define_invalid(r_map, i);
+				}
+			}
+		}
+		else if (mode == MREMAP_MODE_POLY_POLYINTERP_PNORPROJ) {
+			/* We cast our rays randomly, with a pseudo-even distribution (since we spread across tessellated tris,
+			 * with additional weighting based on each tri's relative area).
+			 */
+			RNG *rng = BLI_rng_new(0);
+
+			const size_t numpolys_src = (size_t)dm_src->getNumPolys(dm_src);
+
+			/* Here it's simpler to just allocate for all polys :/ */
+			int *indices = MEM_mallocN(sizeof(*indices) * numpolys_src, __func__);
+			float *weights = MEM_mallocN(sizeof(*weights) * numpolys_src, __func__);
+
+			size_t tmp_poly_size = MREMAP_DEFAULT_BUFSIZE;
+			float (*poly_vcos_2d)[2] = MEM_mallocN(sizeof(*poly_vcos_2d) * tmp_poly_size, __func__);
+			/* Tessellated 2D poly, always (num_loops - 2) triangles. */
+			int (*tri_vidx_2d)[3] = MEM_mallocN(sizeof(*tri_vidx_2d) * (tmp_poly_size - 2), __func__);
+
+			for (i = 0; i < numpolys_dst; i++) {
+				/* For each dst poly, we sample some rays from it (2D grid in pnor space)
+				 * and use their hits to interpolate from source polys. */
+				/* Note: dst poly is early-converted into src space! */
+				MPoly *mp = &polys_dst[i];
+				float tmp_co[3], tmp_no[3];
+
+				int tot_rays, done_rays = 0;
+				float poly_area_2d_inv, done_area = 0.0f;
+
+				const float zvec[3] = {0.0f, 0.0f, 1.0f};
+				float pcent_dst[3];
+				float to_pnor_2d_mat[3][3], from_pnor_2d_mat[3][3];
+				float poly_dst_2d_min[2], poly_dst_2d_max[2], poly_dst_2d_z;
+				float poly_dst_2d_size[2];
+
+				float totweights = 0.0f;
+				float hit_dist_accum = 0.0f;
+				int sources_num = 0;
+				const int tris_num = mp->totloop - 2;
+				int j;
+
+				BKE_mesh_calc_poly_center(mp, &loops_dst[mp->loopstart], verts_dst, pcent_dst);
+				copy_v3_v3(tmp_no, poly_nors_dst[i]);
+				/* We do our transform here, else it'd be redone by raycast helper for each ray, ugh! */
+				if (space_transform) {
+					BLI_space_transform_apply(space_transform, pcent_dst);
+					BLI_space_transform_apply_normal(space_transform, tmp_no);
+				}
+
+				fill_vn_fl(weights, (int)numpolys_src, 0.0f);
+
+				if (UNLIKELY((size_t)mp->totloop > tmp_poly_size)) {
+					tmp_poly_size = (size_t)mp->totloop;
+					poly_vcos_2d = MEM_reallocN(poly_vcos_2d, sizeof(*poly_vcos_2d) * tmp_poly_size);
+					tri_vidx_2d = MEM_reallocN(tri_vidx_2d, sizeof(*tri_vidx_2d) * (tmp_poly_size - 2));
+				}
+
+				rotation_between_vecs_to_mat3(to_pnor_2d_mat, tmp_no, zvec);
+				invert_m3_m3(from_pnor_2d_mat, to_pnor_2d_mat);
+
+				mul_m3_v3(to_pnor_2d_mat, pcent_dst);
+				poly_dst_2d_z = pcent_dst[2];
+
+				/* Get (2D) bounding square of our poly. */
+				INIT_MINMAX2(poly_dst_2d_min, poly_dst_2d_max);
+
+				for (j = 0; j < mp->totloop; j++) {
+					MLoop *ml = &loops_dst[j + mp->loopstart];
+					copy_v3_v3(tmp_co, verts_dst[ml->v].co);
+					if (space_transform) {
+						BLI_space_transform_apply(space_transform, tmp_co);
+					}
+					mul_v2_m3v3(poly_vcos_2d[j], to_pnor_2d_mat, tmp_co);
+					minmax_v2v2_v2(poly_dst_2d_min, poly_dst_2d_max, poly_vcos_2d[j]);
+				}
+
+				/* We adjust our ray-casting grid to ray_radius (the smaller, the more rays are cast),
+				 * with lower/upper bounds. */
+				sub_v2_v2v2(poly_dst_2d_size, poly_dst_2d_max, poly_dst_2d_min);
+
+				if (ray_radius) {
+					tot_rays = (int)((max_ff(poly_dst_2d_size[0], poly_dst_2d_size[1]) / ray_radius) + 0.5f);
+					CLAMP(tot_rays, MREMAP_RAYCAST_TRI_SAMPLES_MIN, MREMAP_RAYCAST_TRI_SAMPLES_MAX);
+				}
+				else {
+					/* If no radius (pure rays), give max number of rays! */
+					tot_rays = MREMAP_RAYCAST_TRI_SAMPLES_MIN;
+				}
+				tot_rays *= tot_rays;
+
+				poly_area_2d_inv = 1.0f / area_poly_v2((const float(*)[2])poly_vcos_2d, (unsigned int)mp->totloop);
+
+				/* Tessellate our poly. */
+				if (mp->totloop == 3) {
+					tri_vidx_2d[0][0] = 0;
+					tri_vidx_2d[0][1] = 1;
+					tri_vidx_2d[0][2] = 2;
+				}
+				if (mp->totloop == 4) {
+					tri_vidx_2d[0][0] = 0;
+					tri_vidx_2d[0][1] = 1;
+					tri_vidx_2d[0][2] = 2;
+					tri_vidx_2d[1][0] = 0;
+					tri_vidx_2d[1][1] = 2;
+					tri_vidx_2d[1][2] = 3;
+				}
+				else {
+					BLI_polyfill_calc((const float(*)[2])poly_vcos_2d, (unsigned int)mp->totloop, -1,
+					                  (unsigned int (*)[3])tri_vidx_2d);
+				}
+
+				for (j = 0; j < tris_num; j++) {
+					float *v1 = poly_vcos_2d[tri_vidx_2d[j][0]];
+					float *v2 = poly_vcos_2d[tri_vidx_2d[j][1]];
+					float *v3 = poly_vcos_2d[tri_vidx_2d[j][2]];
+					int rays_num;
+
+					/* All this allows us to get 'absolute' number of rays for each tri, avoiding accumulating
+					 * errors over iterations, and helping better even distribution. */
+					done_area += area_tri_v2(v1, v2, v3);
+					rays_num = (int)((float)tot_rays * done_area * poly_area_2d_inv + 0.5f) - done_rays;
+					done_rays += rays_num;
+
+					while (rays_num--) {
+						int n = (ray_radius > 0.0f) ? MREMAP_RAYCAST_APPROXIMATE_NR : 1;
+						float w = 1.0f;
+
+						BLI_rng_get_tri_sample_float_v2(rng, v1, v2, v3, tmp_co);
+
+						tmp_co[2] = poly_dst_2d_z;
+						mul_m3_v3(from_pnor_2d_mat, tmp_co);
+
+						/* At this point, tmp_co is a point on our poly surface, in mesh_src space! */
+						while (n--) {
+							/* Note we handle dest to src space conversion ourself, here! */
+
+							if (mesh_remap_bvhtree_query_raycast(
+							        &treedata, &rayhit, NULL,
+							        tmp_co, tmp_no, ray_radius / w, max_dist, &hit_dist))
+							{
+								weights[tessface_to_poly_map_src[rayhit.index]] += w;
+								totweights += w;
+								hit_dist_accum += hit_dist;
+								break;
+							}
+							/* Next iteration will get bigger radius but smaller weight! */
+							w /= MREMAP_RAYCAST_APPROXIMATE_FAC;
+						}
+					}
+				}
+
+				if (totweights > 0.0f) {
+					for (j = 0; j < (int)numpolys_src; j++) {
+						if (!weights[j]) {
+							continue;
+						}
+						/* Note: sources_num is always <= j! */
+						weights[sources_num] = weights[j] / totweights;
+						indices[sources_num] = j;
+						sources_num++;
+					}
+					mesh_remap_item_define(r_map, i, hit_dist_accum / totweights, 0, sources_num, indices, weights);
+				}
+				else {
+					/* No source for this dest poly! */
+					BKE_mesh_remap_item_define_invalid(r_map, i);
+				}
+			}
+
+			MEM_freeN(tri_vidx_2d);
+			MEM_freeN(poly_vcos_2d);
+			MEM_freeN(indices);
+			MEM_freeN(weights);
+			BLI_rng_free(rng);
+		}
+		else {
+			printf("WARNING! Unsupported mesh-to-mesh poly mapping mode (%d)!\n", mode);
+			memset(r_map->items, 0, sizeof(*r_map->items) * (size_t)numpolys_dst);
+		}
+
+		free_bvhtree_from_mesh(&treedata);
+	}
+}
+
+#undef MREMAP_RAYCAST_APPROXIMATE_NR
+#undef MREMAP_RAYCAST_APPROXIMATE_FAC
+#undef MREMAP_RAYCAST_APPROXIMATE_BVHEPSILON
+#undef MREMAP_RAYCAST_TRI_SAMPLES_MIN
+#undef MREMAP_RAYCAST_TRI_SAMPLES_MAX
+#undef MREMAP_DEFAULT_BUFSIZE
+
+/** \} */