source/blender/blenkernel/intern/mesh_remap.c

   1 /*
   2  * ***** BEGIN GPL LICENSE BLOCK *****
   3  *
   4  * This program is free software; you can redistribute it and/or
   5  * modify it under the terms of the GNU General Public License
   6  * as published by the Free Software Foundation; either version 2
   7  * of the License, or (at your option) any later version.
   8  *
   9  * This program is distributed in the hope that it will be useful,
  10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12  * GNU General Public License for more details.
  13  *
  14  * You should have received a copy of the GNU General Public License
  15  * along with this program; if not, write to the Free Software Foundation,
  16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  17  *
  18  * ***** END GPL LICENSE BLOCK *****
  19  */
  20
  21 /** \file blender/blenkernel/intern/mesh_remap.c
  22  *  \ingroup bke
  23  *
  24  * Functions for mapping data between meshes.
  25  */
  26
  27 #include <limits.h>
  28
  29 #include "MEM_guardedalloc.h"
  30
  31 #include "DNA_meshdata_types.h"
  32
  33 #include "BLI_utildefines.h"
  34 #include "BLI_alloca.h"
  35 #include "BLI_astar.h"
  36 #include "BLI_bitmap.h"
  37 #include "BLI_math.h"
  38 #include "BLI_memarena.h"
  39 #include "BLI_polyfill_2d.h"
  40 #include "BLI_rand.h"
  41
  42 #include "BKE_bvhutils.h"
  43 #include "BKE_customdata.h"
  44 #include "BKE_DerivedMesh.h"
  45 #include "BKE_mesh.h"
  46 #include "BKE_mesh_mapping.h"
  47 #include "BKE_mesh_remap.h"  /* own include */
  48
  49 #include "BLI_strict_flags.h"
  50
  51
  52 /* -------------------------------------------------------------------- */
  53
  54 /** \name Some generic helpers.
  55  * \{ */
  56
  57 static bool mesh_remap_bvhtree_query_nearest(
  58         BVHTreeFromMesh *treedata, BVHTreeNearest *nearest,
  59         const float co[3], const float max_dist_sq, float *r_hit_dist)
  60 {
  61         /* Use local proximity heuristics (to reduce the nearest search). */
  62         if (nearest->index != -1) {
  63                 nearest->dist_sq = min_ff(len_squared_v3v3(co, nearest->co), max_dist_sq);
  64         }
  65         else {
  66                 nearest->dist_sq = max_dist_sq;
  67         }
  68         /* Compute and store result. If invalid (-1 index), keep FLT_MAX dist. */
  69         BLI_bvhtree_find_nearest(treedata->tree, co, nearest, treedata->nearest_callback, treedata);
  70
  71         if ((nearest->index != -1) && (nearest->dist_sq <= max_dist_sq)) {
  72                 *r_hit_dist = sqrtf(nearest->dist_sq);
  73                 return true;
  74         }
  75         else {
  76                 return false;
  77         }
  78 }
  79
  80 static bool mesh_remap_bvhtree_query_raycast(
  81         BVHTreeFromMesh *treedata, BVHTreeRayHit *rayhit,
  82         const float co[3], const float no[3], const float radius, const float max_dist, float *r_hit_dist)
  83 {
  84         BVHTreeRayHit rayhit_tmp;
  85         float inv_no[3];
  86
  87         rayhit->index = -1;
  88         rayhit->dist = max_dist;
  89         BLI_bvhtree_ray_cast(treedata->tree, co, no, radius, rayhit, treedata->raycast_callback, treedata);
  90
  91         /* Also cast in the other direction! */
  92         rayhit_tmp = *rayhit;
  93         negate_v3_v3(inv_no, no);
  94         BLI_bvhtree_ray_cast(treedata->tree, co, inv_no, radius, &rayhit_tmp, treedata->raycast_callback, treedata);
  95         if (rayhit_tmp.dist < rayhit->dist) {
  96                 *rayhit = rayhit_tmp;
  97         }
  98
  99         if ((rayhit->index != -1) && (rayhit->dist <= max_dist)) {
 100                 *r_hit_dist = rayhit->dist;
 101                 return true;
 102         }
 103         else {
 104                 return false;
 105         }
 106 }
 107
 108 /** \} */
 109
 110 /**
 111  * \name Auto-match.
 112  *
 113  * Find transform of a mesh to get best match with another.
 114  * \{ */
 115
 116 /**
 117  * Compute a value of the difference between both given meshes.
 118  * The smaller the result, the better the match.
 119  *
 120  * We return the inverse of the average of the inversed shortest distance from each dst vertex to src ones.
 121  * In other words, beyond a certain (relatively small) distance, all differences have more or less the same weight
 122  * in final result, which allows to reduce influence of a few high differences, in favor of a global good matching.
 123  */
 124 float BKE_mesh_remap_calc_difference_from_dm(
 125         const SpaceTransform *space_transform, const MVert *verts_dst, const int numverts_dst, DerivedMesh *dm_src)
 126 {
 127         BVHTreeFromMesh treedata = {NULL};
 128         BVHTreeNearest nearest = {0};
 129         float hit_dist;
 130
 131         float result = 0.0f;
 132         int i;
 133
 134         bvhtree_from_mesh_get(&treedata, dm_src, BVHTREE_FROM_VERTS, 2);
 135         nearest.index = -1;
 136
 137         for (i = 0; i < numverts_dst; i++) {
 138                 float tmp_co[3];
 139
 140                 copy_v3_v3(tmp_co, verts_dst[i].co);
 141
 142                 /* Convert the vertex to tree coordinates, if needed. */
 143                 if (space_transform) {
 144                         BLI_space_transform_apply(space_transform, tmp_co);
 145                 }
 146
 147                 if (mesh_remap_bvhtree_query_nearest(&treedata, &nearest, tmp_co, FLT_MAX, &hit_dist)) {
 148                         result += 1.0f / (hit_dist + 1.0f);
 149                 }
 150                 else {
 151                         /* No source for this dest vertex! */
 152                         result += 1e-18f;
 153                 }
 154         }
 155
 156         result = ((float)numverts_dst / result) - 1.0f;
 157
 158 //      printf("%s: Computed difference between meshes (the lower the better): %f\n", __func__, result);
 159
 160         return result;
 161 }
 162
 163 /* This helper computes the eigen values & vectors for covariance matrix of all given vertices coordinates.
 164  *
 165  * Those vectors define the 'average ellipsoid' of the mesh (i.e. the 'best fitting' ellipsoid
 166  * containing 50% of the vertices).
 167  *
 168  * Note that it will not perform fantastic in case two or more eigen values are equal (e.g. a cylinder or
 169  * parallelepiped with a square section give two identical eigenvalues, a sphere or tetrahedron give
 170  * three identical ones, etc.), since you cannot really define all axes in those cases. We default to dummy
 171  * generated orthogonal vectors in this case, instead of using eigen vectors.
 172  */
 173 static void mesh_calc_eigen_matrix(
 174         const MVert *verts, const float (*vcos)[3], const int numverts, float r_mat[4][4])
 175 {
 176         float center[3], covmat[3][3];
 177         float eigen_val[3], eigen_vec[3][3];
 178         float (*cos)[3] = NULL;
 179
 180         bool eigen_success;
 181         int i;
 182
 183         if (verts) {
 184                 const MVert *mv;
 185                 float (*co)[3];
 186
 187                 cos = MEM_mallocN(sizeof(*cos) * (size_t)numverts, __func__);
 188                 for (i = 0, co = cos, mv = verts; i < numverts; i++, co++, mv++) {
 189                         copy_v3_v3(*co, mv->co);
 190                 }
 191                 /* TODO(sergey): For until we officially drop all compilers which
 192                  * doesn't handle casting correct we use workaround to avoid explicit
 193                  * cast here.
 194                  */
 195                 vcos = (void *)cos;
 196         }
 197         unit_m4(r_mat);
 198
 199         /* Note: here we apply sample correction to covariance matrix, since we consider the vertices as a sample
 200          *       of the whole 'surface' population of our mesh... */
 201         BLI_covariance_m3_v3n(vcos, numverts, true, covmat, center);
 202
 203         if (cos) {
 204                 MEM_freeN(cos);
 205         }
 206
 207         eigen_success = BLI_eigen_solve_selfadjoint_m3((const float (*)[3])covmat, eigen_val, eigen_vec);
 208         BLI_assert(eigen_success);
 209         UNUSED_VARS_NDEBUG(eigen_success);
 210
 211         /* Special handling of cases where some eigen values are (nearly) identical. */
 212         if (compare_ff_relative(eigen_val[0], eigen_val[1], FLT_EPSILON, 64)) {
 213                 if (compare_ff_relative(eigen_val[0], eigen_val[2], FLT_EPSILON, 64)) {
 214                         /* No preferred direction, that set of vertices has a spherical average,
 215                          * so we simply returned scaled/translated identity matrix (with no rotation). */
 216                         unit_m3(eigen_vec);
 217                 }
 218                 else {
 219                         /* Ellipsoid defined by eigen values/vectors has a spherical section,
 220                          * we can only define one axis from eigen_vec[2] (two others computed eigen vecs
 221                          * are not so nice for us here, they tend to 'randomly' rotate around valid one).
 222                          * Note that eigen vectors as returned by BLI_eigen_solve_selfadjoint_m3() are normalized. */
 223                         ortho_basis_v3v3_v3(eigen_vec[0], eigen_vec[1], eigen_vec[2]);
 224                 }
 225         }
 226         else if (compare_ff_relative(eigen_val[0], eigen_val[2], FLT_EPSILON, 64)) {
 227                 /* Same as above, but with eigen_vec[1] as valid axis. */
 228                 ortho_basis_v3v3_v3(eigen_vec[2], eigen_vec[0], eigen_vec[1]);
 229         }
 230         else if (compare_ff_relative(eigen_val[1], eigen_val[2], FLT_EPSILON, 64)) {
 231                 /* Same as above, but with eigen_vec[0] as valid axis. */
 232                 ortho_basis_v3v3_v3(eigen_vec[1], eigen_vec[2], eigen_vec[0]);
 233         }
 234
 235         for (i = 0; i < 3; i++) {
 236                 float evi = eigen_val[i];
 237
 238                 /* Protect against 1D/2D degenerated cases! */
 239                 /* Note: not sure why we need square root of eigen values here (which are equivalent to singular values,
 240                  * as far as I have understood), but it seems to heavily reduce (if not completely nullify)
 241                  * the error due to non-uniform scalings... */
 242                 evi = (evi < 1e-6f && evi > -1e-6f) ? ((evi < 0.0f) ? -1e-3f : 1e-3f) : sqrtf_signed(evi);
 243                 mul_v3_fl(eigen_vec[i], evi);
 244         }
 245
 246         copy_m4_m3(r_mat, eigen_vec);
 247         copy_v3_v3(r_mat[3], center);
 248 }
 249
 250 /**
 251  * Set r_space_transform so that best bbox of dst matches best bbox of src.
 252  */
 253 void BKE_mesh_remap_find_best_match_from_dm(
 254         const MVert *verts_dst, const int numverts_dst, DerivedMesh *dm_src, SpaceTransform *r_space_transform)
 255 {
 256         /* Note that those are done so that we successively get actual mirror matrix (by multiplication of columns)... */
 257         const float mirrors[][3] = {
 258             {-1.0f,  1.0f,  1.0f},  /* -> -1,  1,  1 */
 259             { 1.0f, -1.0f,  1.0f},  /* -> -1, -1,  1 */
 260             { 1.0f,  1.0f, -1.0f},  /* -> -1, -1, -1 */
 261             { 1.0f, -1.0f,  1.0f},  /* -> -1,  1, -1 */
 262             {-1.0f,  1.0f,  1.0f},  /* ->  1,  1, -1 */
 263             { 1.0f, -1.0f,  1.0f},  /* ->  1, -1, -1 */
 264             { 1.0f,  1.0f, -1.0f},  /* ->  1, -1,  1 */
 265             {0.0f, 0.0f, 0.0f},
 266         };
 267         const float (*mirr)[3];
 268
 269         float mat_src[4][4], mat_dst[4][4], best_mat_dst[4][4];
 270         float best_match = FLT_MAX, match;
 271
 272         const int numverts_src = dm_src->getNumVerts(dm_src);
 273         float (*vcos_src)[3] = MEM_mallocN(sizeof(*vcos_src) * (size_t)numverts_src, __func__);
 274         dm_src->getVertCos(dm_src, vcos_src);
 275
 276         mesh_calc_eigen_matrix(NULL, (const float (*)[3])vcos_src, numverts_src, mat_src);
 277         mesh_calc_eigen_matrix(verts_dst, NULL, numverts_dst, mat_dst);
 278
 279         BLI_space_transform_global_from_matrices(r_space_transform, mat_dst, mat_src);
 280         match = BKE_mesh_remap_calc_difference_from_dm(r_space_transform, verts_dst, numverts_dst, dm_src);
 281         best_match = match;
 282         copy_m4_m4(best_mat_dst, mat_dst);
 283
 284         /* And now, we have to check the otehr sixth possible mirrored versions... */
 285         for (mirr = mirrors; (*mirr)[0]; mirr++) {
 286                 mul_v3_fl(mat_dst[0], (*mirr)[0]);
 287                 mul_v3_fl(mat_dst[1], (*mirr)[1]);
 288                 mul_v3_fl(mat_dst[2], (*mirr)[2]);
 289
 290                 BLI_space_transform_global_from_matrices(r_space_transform, mat_dst, mat_src);
 291                 match = BKE_mesh_remap_calc_difference_from_dm(r_space_transform, verts_dst, numverts_dst, dm_src);
 292                 if (match < best_match) {
 293                         best_match = match;
 294                         copy_m4_m4(best_mat_dst, mat_dst);
 295                 }
 296         }
 297
 298         BLI_space_transform_global_from_matrices(r_space_transform, best_mat_dst, mat_src);
 299
 300         MEM_freeN(vcos_src);
 301 }
 302
 303 /** \} */
 304
 305 /** \name Mesh to mesh mapping
 306  * \{ */
 307
 308 void BKE_mesh_remap_init(MeshPairRemap *map, const int items_num)
 309 {
 310         MemArena *mem = BLI_memarena_new(BLI_MEMARENA_STD_BUFSIZE, __func__);
 311
 312         BKE_mesh_remap_free(map);
 313
 314         map->items = BLI_memarena_alloc(mem, sizeof(*map->items) * (size_t)items_num);
 315         map->items_num = items_num;
 316
 317         map->mem = mem;
 318 }
 319
 320 void BKE_mesh_remap_free(MeshPairRemap *map)
 321 {
 322         if (map->mem) {
 323                 BLI_memarena_free((MemArena *)map->mem);
 324         }
 325
 326         map->items_num = 0;
 327         map->items = NULL;
 328         map->mem = NULL;
 329 }
 330
 331 static void mesh_remap_item_define(
 332         MeshPairRemap *map, const int index, const float UNUSED(hit_dist), const int island,
 333         const int sources_num, const int *indices_src, const float *weights_src)
 334 {
 335         MeshPairRemapItem *mapit = &map->items[index];
 336         MemArena *mem = map->mem;
 337
 338         if (sources_num) {
 339                 mapit->sources_num = sources_num;
 340                 mapit->indices_src = BLI_memarena_alloc(mem, sizeof(*mapit->indices_src) * (size_t)sources_num);
 341                 memcpy(mapit->indices_src, indices_src, sizeof(*mapit->indices_src) * (size_t)sources_num);
 342                 mapit->weights_src = BLI_memarena_alloc(mem, sizeof(*mapit->weights_src) * (size_t)sources_num);
 343                 memcpy(mapit->weights_src, weights_src, sizeof(*mapit->weights_src) * (size_t)sources_num);
 344         }
 345         else {
 346                 mapit->sources_num = 0;
 347                 mapit->indices_src = NULL;
 348                 mapit->weights_src = NULL;
 349         }
 350         /* UNUSED */
 351         // mapit->hit_dist = hit_dist;
 352         mapit->island = island;
 353 }
 354
 355 void BKE_mesh_remap_item_define_invalid(MeshPairRemap *map, const int index)
 356 {
 357         mesh_remap_item_define(map, index, FLT_MAX, 0, 0, NULL, NULL);
 358 }
 359
 360 static int mesh_remap_interp_poly_data_get(
 361         const MPoly *mp, MLoop *mloops, const float (*vcos_src)[3], const float point[3],
 362         size_t *buff_size, float (**vcos)[3], const bool use_loops, int **indices, float **weights,
 363         const bool do_weights, int *r_closest_index)
 364 {
 365         MLoop *ml;
 366         float (*vco)[3];
 367         float ref_dist_sq = FLT_MAX;
 368         int *index;
 369         const int sources_num = mp->totloop;
 370         int i;
 371
 372         if ((size_t)sources_num > *buff_size) {
 373                 *buff_size = (size_t)sources_num;
 374                 *vcos = MEM_reallocN(*vcos, sizeof(**vcos) * *buff_size);
 375                 *indices = MEM_reallocN(*indices, sizeof(**indices) * *buff_size);
 376                 if (do_weights) {
 377                         *weights = MEM_reallocN(*weights, sizeof(**weights) * *buff_size);
 378                 }
 379         }
 380
 381         for (i = 0, ml = &mloops[mp->loopstart], vco = *vcos, index = *indices; i < sources_num; i++, ml++, vco++, index++) {
 382                 *index = use_loops ? (int)mp->loopstart + i : (int)ml->v;
 383                 copy_v3_v3(*vco, vcos_src[ml->v]);
 384                 if (r_closest_index) {
 385                         /* Find closest vert/loop in this case. */
 386                         const float dist_sq = len_squared_v3v3(point, *vco);
 387                         if (dist_sq < ref_dist_sq) {
 388                                 ref_dist_sq = dist_sq;
 389                                 *r_closest_index = *index;
 390                         }
 391                 }
 392         }
 393
 394         if (do_weights) {
 395                 interp_weights_poly_v3(*weights, *vcos, sources_num, point);
 396         }
 397
 398         return sources_num;
 399 }
 400
 401 /* Little helper when dealing with source islands */
 402 typedef struct IslandResult {
 403         float factor;           /* A factor, based on which best island for a given set of elements will be selected. */
 404         int   index_src;        /* Index of the source. */
 405         float hit_dist;         /* The actual hit distance. */
 406         float hit_point[3];     /* The hit point, if relevant. */
 407 } IslandResult;
 408
 409 /* Note about all bvh/raycasting stuff below:
 410  *      * We must use our ray radius as BVH epsilon too, else rays not hitting anything but 'passing near' an item
 411  *        would be missed (since BVH handling would not detect them, 'refining' callbacks won't be executed,
 412  *        even though they would return a valid hit).
 413  *      * However, in 'islands' case where each hit gets a weight, 'precise' hits should have a better weight than
 414  *        'approximate' hits. To address that, we simplify things with:
 415  *        ** A first raycast with default, given rayradius;
 416  *        ** If first one fails, we do more raycasting with bigger radius, but if hit is found
 417  *           it will get smaller weight.
 418  *        This only concerns loops, currently (because of islands), and 'sampled' edges/polys norproj.
 419  */
 420
 421 /* At most n raycasts per 'real' ray. */
 422 #define MREMAP_RAYCAST_APPROXIMATE_NR 3
 423 /* Each approximated raycasts will have n times bigger radius than previous one. */
 424 #define MREMAP_RAYCAST_APPROXIMATE_FAC 5.0f
 425
 426 /* min 16 rays/face, max 400. */
 427 #define MREMAP_RAYCAST_TRI_SAMPLES_MIN 4
 428 #define MREMAP_RAYCAST_TRI_SAMPLES_MAX 20
 429
 430 /* Will be enough in 99% of cases. */
 431 #define MREMAP_DEFAULT_BUFSIZE 32
 432
 433 void BKE_mesh_remap_calc_verts_from_dm(
 434         const int mode, const SpaceTransform *space_transform, const float max_dist, const float ray_radius,
 435         const MVert *verts_dst, const int numverts_dst, const bool UNUSED(dirty_nors_dst), DerivedMesh *dm_src,
 436         MeshPairRemap *r_map)
 437 {
 438         const float full_weight = 1.0f;
 439         const float max_dist_sq = max_dist * max_dist;
 440         int i;
 441
 442         BLI_assert(mode & MREMAP_MODE_VERT);
 443
 444         BKE_mesh_remap_init(r_map, numverts_dst);
 445
 446         if (mode == MREMAP_MODE_TOPOLOGY) {
 447                 BLI_assert(numverts_dst == dm_src->getNumVerts(dm_src));
 448                 for (i = 0; i < numverts_dst; i++) {
 449                         mesh_remap_item_define(r_map, i, FLT_MAX, 0, 1, &i, &full_weight);
 450                 }
 451         }
 452         else {
 453                 BVHTreeFromMesh treedata = {NULL};
 454                 BVHTreeNearest nearest = {0};
 455                 BVHTreeRayHit rayhit = {0};
 456                 float hit_dist;
 457                 float tmp_co[3], tmp_no[3];
 458
 459                 if (mode == MREMAP_MODE_VERT_NEAREST) {
 460                         bvhtree_from_mesh_get(&treedata, dm_src, BVHTREE_FROM_VERTS, 2);
 461                         nearest.index = -1;
 462
 463                         for (i = 0; i < numverts_dst; i++) {
 464                                 copy_v3_v3(tmp_co, verts_dst[i].co);
 465
 466                                 /* Convert the vertex to tree coordinates, if needed. */
 467                                 if (space_transform) {
 468                                         BLI_space_transform_apply(space_transform, tmp_co);
 469                                 }
 470
 471                                 if (mesh_remap_bvhtree_query_nearest(&treedata, &nearest, tmp_co, max_dist_sq, &hit_dist)) {
 472                                         mesh_remap_item_define(r_map, i, hit_dist, 0, 1, &nearest.index, &full_weight);
 473                                 }
 474                                 else {
 475                                         /* No source for this dest vertex! */
 476                                         BKE_mesh_remap_item_define_invalid(r_map, i);
 477                                 }
 478                         }
 479                 }
 480                 else if (ELEM(mode, MREMAP_MODE_VERT_EDGE_NEAREST, MREMAP_MODE_VERT_EDGEINTERP_NEAREST)) {
 481                         MEdge *edges_src = dm_src->getEdgeArray(dm_src);
 482                         float (*vcos_src)[3] = MEM_mallocN(sizeof(*vcos_src) * (size_t)dm_src->getNumVerts(dm_src), __func__);
 483                         dm_src->getVertCos(dm_src, vcos_src);
 484
 485                         bvhtree_from_mesh_get(&treedata, dm_src, BVHTREE_FROM_EDGES, 2);
 486                         nearest.index = -1;
 487
 488                         for (i = 0; i < numverts_dst; i++) {
 489                                 copy_v3_v3(tmp_co, verts_dst[i].co);
 490
 491                                 /* Convert the vertex to tree coordinates, if needed. */
 492                                 if (space_transform) {
 493                                         BLI_space_transform_apply(space_transform, tmp_co);
 494                                 }
 495
 496                                 if (mesh_remap_bvhtree_query_nearest(&treedata, &nearest, tmp_co, max_dist_sq, &hit_dist)) {
 497                                         MEdge *me = &edges_src[nearest.index];
 498                                         const float *v1cos = vcos_src[me->v1];
 499                                         const float *v2cos = vcos_src[me->v2];
 500
 501                                         if (mode == MREMAP_MODE_VERT_EDGE_NEAREST) {
 502                                                 const float dist_v1 = len_squared_v3v3(tmp_co, v1cos);
 503                                                 const float dist_v2 = len_squared_v3v3(tmp_co, v2cos);
 504                                                 const int index = (int)((dist_v1 > dist_v2) ? me->v2 : me->v1);
 505                                                 mesh_remap_item_define(r_map, i, hit_dist, 0, 1, &index, &full_weight);
 506                                         }
 507                                         else if (mode == MREMAP_MODE_VERT_EDGEINTERP_NEAREST) {
 508                                                 int indices[2];
 509                                                 float weights[2];
 510
 511                                                 indices[0] = (int)me->v1;
 512                                                 indices[1] = (int)me->v2;
 513
 514                                                 /* Weight is inverse of point factor here... */
 515                                                 weights[0] = line_point_factor_v3(tmp_co, v2cos, v1cos);
 516                                                 CLAMP(weights[0], 0.0f, 1.0f);
 517                                                 weights[1] = 1.0f - weights[0];
 518
 519                                                 mesh_remap_item_define(r_map, i, hit_dist, 0, 2, indices, weights);
 520                                         }
 521                                 }
 522                                 else {
 523                                         /* No source for this dest vertex! */
 524                                         BKE_mesh_remap_item_define_invalid(r_map, i);
 525                                 }
 526                         }
 527
 528                         MEM_freeN(vcos_src);
 529                 }
 530                 else if (ELEM(mode, MREMAP_MODE_VERT_POLY_NEAREST, MREMAP_MODE_VERT_POLYINTERP_NEAREST,
 531                                     MREMAP_MODE_VERT_POLYINTERP_VNORPROJ))
 532                 {
 533                         MPoly *polys_src = dm_src->getPolyArray(dm_src);
 534                         MLoop *loops_src = dm_src->getLoopArray(dm_src);
 535                         float (*vcos_src)[3] = MEM_mallocN(sizeof(*vcos_src) * (size_t)dm_src->getNumVerts(dm_src), __func__);
 536
 537                         size_t tmp_buff_size = MREMAP_DEFAULT_BUFSIZE;
 538                         float (*vcos)[3] = MEM_mallocN(sizeof(*vcos) * tmp_buff_size, __func__);
 539                         int *indices = MEM_mallocN(sizeof(*indices) * tmp_buff_size, __func__);
 540                         float *weights = MEM_mallocN(sizeof(*weights) * tmp_buff_size, __func__);
 541
 542                         dm_src->getVertCos(dm_src, vcos_src);
 543
 544                         bvhtree_from_mesh_get(&treedata, dm_src, BVHTREE_FROM_LOOPTRI, 2);
 545
 546                         if (mode == MREMAP_MODE_VERT_POLYINTERP_VNORPROJ) {
 547                                 treedata.sphere_radius = ray_radius;
 548                                 for (i = 0; i < numverts_dst; i++) {
 549                                         copy_v3_v3(tmp_co, verts_dst[i].co);
 550                                         normal_short_to_float_v3(tmp_no, verts_dst[i].no);
 551
 552                                         /* Convert the vertex to tree coordinates, if needed. */
 553                                         if (space_transform) {
 554                                                 BLI_space_transform_apply(space_transform, tmp_co);
 555                                                 BLI_space_transform_apply_normal(space_transform, tmp_no);
 556                                         }
 557
 558                                         if (mesh_remap_bvhtree_query_raycast(
 559                                                 &treedata, &rayhit, tmp_co, tmp_no, ray_radius, max_dist, &hit_dist))
 560                                         {
 561                                                 const MLoopTri *lt = &treedata.looptri[rayhit.index];
 562                                                 MPoly *mp_src = &polys_src[lt->poly];
 563                                                 const int sources_num = mesh_remap_interp_poly_data_get(
 564                                                         mp_src, loops_src, (const float (*)[3])vcos_src, rayhit.co,
 565                                                         &tmp_buff_size, &vcos, false, &indices, &weights, true, NULL);
 566
 567                                                 mesh_remap_item_define(r_map, i, hit_dist, 0, sources_num, indices, weights);
 568                                         }
 569                                         else {
 570                                                 /* No source for this dest vertex! */
 571                                                 BKE_mesh_remap_item_define_invalid(r_map, i);
 572                                         }
 573                                 }
 574                         }
 575                         else {
 576                                 nearest.index = -1;
 577
 578                                 for (i = 0; i < numverts_dst; i++) {
 579                                         copy_v3_v3(tmp_co, verts_dst[i].co);
 580
 581                                         /* Convert the vertex to tree coordinates, if needed. */
 582                                         if (space_transform) {
 583                                                 BLI_space_transform_apply(space_transform, tmp_co);
 584                                         }
 585
 586                                         if (mesh_remap_bvhtree_query_nearest(&treedata, &nearest, tmp_co, max_dist_sq, &hit_dist)) {
 587                                                 const MLoopTri *lt = &treedata.looptri[nearest.index];
 588                                                 MPoly *mp = &polys_src[lt->poly];
 589
 590                                                 if (mode == MREMAP_MODE_VERT_POLY_NEAREST) {
 591                                                         int index;
 592                                                         mesh_remap_interp_poly_data_get(
 593                                                                 mp, loops_src, (const float (*)[3])vcos_src, nearest.co,
 594                                                                 &tmp_buff_size, &vcos, false, &indices, &weights, false,
 595                                                                 &index);
 596
 597                                                         mesh_remap_item_define(r_map, i, hit_dist, 0, 1, &index, &full_weight);
 598                                                 }
 599                                                 else if (mode == MREMAP_MODE_VERT_POLYINTERP_NEAREST) {
 600                                                         const int sources_num = mesh_remap_interp_poly_data_get(
 601                                                                 mp, loops_src, (const float (*)[3])vcos_src, nearest.co,
 602                                                                 &tmp_buff_size, &vcos, false, &indices, &weights, true,
 603                                                                 NULL);
 604
 605                                                         mesh_remap_item_define(r_map, i, hit_dist, 0, sources_num, indices, weights);
 606                                                 }
 607                                         }
 608                                         else {
 609                                                 /* No source for this dest vertex! */
 610                                                 BKE_mesh_remap_item_define_invalid(r_map, i);
 611                                         }
 612                                 }
 613                         }
 614
 615                         MEM_freeN(vcos_src);
 616                         MEM_freeN(vcos);
 617                         MEM_freeN(indices);
 618                         MEM_freeN(weights);
 619                 }
 620                 else {
 621                         printf("WARNING! Unsupported mesh-to-mesh vertex mapping mode (%d)!\n", mode);
 622                         memset(r_map->items, 0, sizeof(*r_map->items) * (size_t)numverts_dst);
 623                 }
 624
 625                 free_bvhtree_from_mesh(&treedata);
 626         }
 627 }
 628
 629 void BKE_mesh_remap_calc_edges_from_dm(
 630         const int mode, const SpaceTransform *space_transform, const float max_dist, const float ray_radius,
 631         const MVert *verts_dst, const int numverts_dst, const MEdge *edges_dst, const int numedges_dst,
 632         const bool UNUSED(dirty_nors_dst), DerivedMesh *dm_src, MeshPairRemap *r_map)
 633 {
 634         const float full_weight = 1.0f;
 635         const float max_dist_sq = max_dist * max_dist;
 636         int i;
 637
 638         BLI_assert(mode & MREMAP_MODE_EDGE);
 639
 640         BKE_mesh_remap_init(r_map, numedges_dst);
 641
 642         if (mode == MREMAP_MODE_TOPOLOGY) {
 643                 BLI_assert(numedges_dst == dm_src->getNumEdges(dm_src));
 644                 for (i = 0; i < numedges_dst; i++) {
 645                         mesh_remap_item_define(r_map, i, FLT_MAX, 0, 1, &i, &full_weight);
 646                 }
 647         }
 648         else {
 649                 BVHTreeFromMesh treedata = {NULL};
 650                 BVHTreeNearest nearest = {0};
 651                 BVHTreeRayHit rayhit = {0};
 652                 float hit_dist;
 653                 float tmp_co[3], tmp_no[3];
 654
 655                 if (mode == MREMAP_MODE_EDGE_VERT_NEAREST) {
 656                         const int num_verts_src = dm_src->getNumVerts(dm_src);
 657                         const int num_edges_src = dm_src->getNumEdges(dm_src);
 658                         MEdge *edges_src = dm_src->getEdgeArray(dm_src);
 659                         float (*vcos_src)[3] = MEM_mallocN(sizeof(*vcos_src) * (size_t)dm_src->getNumVerts(dm_src), __func__);
 660
 661                         MeshElemMap *vert_to_edge_src_map;
 662                         int         *vert_to_edge_src_map_mem;
 663
 664                         struct {
 665                                 float hit_dist;
 666                                 int   index;
 667                         } *v_dst_to_src_map = MEM_mallocN(sizeof(*v_dst_to_src_map) * (size_t)numverts_dst, __func__);
 668
 669                         for (i = 0; i < numverts_dst; i++) {
 670                                 v_dst_to_src_map[i].hit_dist = -1.0f;
 671                         }
 672
 673                         BKE_mesh_vert_edge_map_create(&vert_to_edge_src_map, &vert_to_edge_src_map_mem,
 674                                                       edges_src, num_verts_src, num_edges_src);
 675
 676                         dm_src->getVertCos(dm_src, vcos_src);
 677
 678                         bvhtree_from_mesh_get(&treedata, dm_src, BVHTREE_FROM_VERTS, 2);
 679                         nearest.index = -1;
 680
 681                         for (i = 0; i < numedges_dst; i++) {
 682                                 const MEdge *e_dst = &edges_dst[i];
 683                                 float best_totdist = FLT_MAX;
 684                                 int best_eidx_src = -1;
 685                                 int j = 2;
 686
 687                                 while (j--) {
 688                                         const unsigned int vidx_dst = j ? e_dst->v1 : e_dst->v2;
 689
 690                                         /* Compute closest verts only once! */
 691                                         if (v_dst_to_src_map[vidx_dst].hit_dist == -1.0f) {
 692                                                 copy_v3_v3(tmp_co, verts_dst[vidx_dst].co);
 693
 694                                                 /* Convert the vertex to tree coordinates, if needed. */
 695                                                 if (space_transform) {
 696                                                         BLI_space_transform_apply(space_transform, tmp_co);
 697                                                 }
 698
 699                                                 if (mesh_remap_bvhtree_query_nearest(&treedata, &nearest, tmp_co, max_dist_sq, &hit_dist)) {
 700                                                         v_dst_to_src_map[vidx_dst].hit_dist = hit_dist;
 701                                                         v_dst_to_src_map[vidx_dst].index = nearest.index;
 702                                                 }
 703                                                 else {
 704                                                         /* No source for this dest vert! */
 705                                                         v_dst_to_src_map[vidx_dst].hit_dist = FLT_MAX;
 706                                                 }
 707                                         }
 708                                 }
 709
 710                                 /* Now, check all source edges of closest sources vertices, and select the one giving the smallest
 711                                  * total verts-to-verts distance. */
 712                                 for (j = 2; j--;) {
 713                                         const unsigned int vidx_dst = j ? e_dst->v1 : e_dst->v2;
 714                                         const float first_dist = v_dst_to_src_map[vidx_dst].hit_dist;
 715                                         const int vidx_src = v_dst_to_src_map[vidx_dst].index;
 716                                         int *eidx_src, k;
 717
 718                                         if (vidx_src < 0) {
 719                                                 continue;
 720                                         }
 721
 722                                         eidx_src = vert_to_edge_src_map[vidx_src].indices;
 723                                         k = vert_to_edge_src_map[vidx_src].count;
 724
 725                                         for (; k--; eidx_src++) {
 726                                                 MEdge *e_src = &edges_src[*eidx_src];
 727                                                 const float *other_co_src = vcos_src[BKE_mesh_edge_other_vert(e_src, vidx_src)];
 728                                                 const float *other_co_dst = verts_dst[BKE_mesh_edge_other_vert(e_dst, (int)vidx_dst)].co;
 729                                                 const float totdist = first_dist + len_v3v3(other_co_src, other_co_dst);
 730
 731                                                 if (totdist < best_totdist) {
 732                                                         best_totdist = totdist;
 733                                                         best_eidx_src = *eidx_src;
 734                                                 }
 735                                         }
 736                                 }
 737
 738                                 if (best_eidx_src >= 0) {
 739                                         const float *co1_src = vcos_src[edges_src[best_eidx_src].v1];
 740                                         const float *co2_src = vcos_src[edges_src[best_eidx_src].v2];
 741                                         const float *co1_dst = verts_dst[e_dst->v1].co;
 742                                         const float *co2_dst = verts_dst[e_dst->v2].co;
 743                                         float co_src[3], co_dst[3];
 744
 745                                         /* TODO: would need an isect_seg_seg_v3(), actually! */
 746                                         const int isect_type = isect_line_line_v3(co1_src, co2_src, co1_dst, co2_dst, co_src, co_dst);
 747                                         if (isect_type != 0) {
 748                                                 const float fac_src = line_point_factor_v3(co_src, co1_src, co2_src);
 749                                                 const float fac_dst = line_point_factor_v3(co_dst, co1_dst, co2_dst);
 750                                                 if (fac_src < 0.0f) {
 751                                                         copy_v3_v3(co_src, co1_src);
 752                                                 }
 753                                                 else if (fac_src > 1.0f) {
 754                                                         copy_v3_v3(co_src, co2_src);
 755                                                 }
 756                                                 if (fac_dst < 0.0f) {
 757                                                         copy_v3_v3(co_dst, co1_dst);
 758                                                 }
 759                                                 else if (fac_dst > 1.0f) {
 760                                                         copy_v3_v3(co_dst, co2_dst);
 761                                                 }
 762                                         }
 763                                         hit_dist = len_v3v3(co_dst, co_src);
 764                                         mesh_remap_item_define(r_map, i, hit_dist, 0, 1, &best_eidx_src, &full_weight);
 765                                 }
 766                                 else {
 767                                         /* No source for this dest edge! */
 768                                         BKE_mesh_remap_item_define_invalid(r_map, i);
 769                                 }
 770                         }
 771
 772                         MEM_freeN(vcos_src);
 773                         MEM_freeN(v_dst_to_src_map);
 774                         MEM_freeN(vert_to_edge_src_map);
 775                         MEM_freeN(vert_to_edge_src_map_mem);
 776                 }
 777                 else if (mode == MREMAP_MODE_EDGE_NEAREST) {
 778                         bvhtree_from_mesh_get(&treedata, dm_src, BVHTREE_FROM_EDGES, 2);
 779                         nearest.index = -1;
 780
 781                         for (i = 0; i < numedges_dst; i++) {
 782                                 interp_v3_v3v3(tmp_co, verts_dst[edges_dst[i].v1].co, verts_dst[edges_dst[i].v2].co, 0.5f);
 783
 784                                 /* Convert the vertex to tree coordinates, if needed. */
 785                                 if (space_transform) {
 786                                         BLI_space_transform_apply(space_transform, tmp_co);
 787                                 }
 788
 789                                 if (mesh_remap_bvhtree_query_nearest(&treedata, &nearest, tmp_co, max_dist_sq, &hit_dist)) {
 790                                         mesh_remap_item_define(r_map, i, hit_dist, 0, 1, &nearest.index, &full_weight);
 791                                 }
 792                                 else {
 793                                         /* No source for this dest edge! */
 794                                         BKE_mesh_remap_item_define_invalid(r_map, i);
 795                                 }
 796                         }
 797                 }
 798                 else if (mode == MREMAP_MODE_EDGE_POLY_NEAREST) {
 799                         MEdge *edges_src = dm_src->getEdgeArray(dm_src);
 800                         MPoly *polys_src = dm_src->getPolyArray(dm_src);
 801                         MLoop *loops_src = dm_src->getLoopArray(dm_src);
 802                         float (*vcos_src)[3] = MEM_mallocN(sizeof(*vcos_src) * (size_t)dm_src->getNumVerts(dm_src), __func__);
 803
 804                         dm_src->getVertCos(dm_src, vcos_src);
 805                         bvhtree_from_mesh_get(&treedata, dm_src, BVHTREE_FROM_LOOPTRI, 2);
 806
 807                         for (i = 0; i < numedges_dst; i++) {
 808                                 interp_v3_v3v3(tmp_co, verts_dst[edges_dst[i].v1].co, verts_dst[edges_dst[i].v2].co, 0.5f);
 809
 810                                 /* Convert the vertex to tree coordinates, if needed. */
 811                                 if (space_transform) {
 812                                         BLI_space_transform_apply(space_transform, tmp_co);
 813                                 }
 814
 815                                 if (mesh_remap_bvhtree_query_nearest(&treedata, &nearest, tmp_co, max_dist_sq, &hit_dist)) {
 816                                         const MLoopTri *lt = &treedata.looptri[nearest.index];
 817                                         MPoly *mp_src = &polys_src[lt->poly];
 818                                         MLoop *ml_src = &loops_src[mp_src->loopstart];
 819                                         int nloops = mp_src->totloop;
 820                                         float best_dist_sq = FLT_MAX;
 821                                         int best_eidx_src = -1;
 822
 823                                         for (; nloops--; ml_src++) {
 824                                                 MEdge *me_src = &edges_src[ml_src->e];
 825                                                 float *co1_src = vcos_src[me_src->v1];
 826                                                 float *co2_src = vcos_src[me_src->v2];
 827                                                 float co_src[3];
 828                                                 float dist_sq;
 829
 830                                                 interp_v3_v3v3(co_src, co1_src, co2_src, 0.5f);
 831                                                 dist_sq = len_squared_v3v3(tmp_co, co_src);
 832                                                 if (dist_sq < best_dist_sq) {
 833                                                         best_dist_sq = dist_sq;
 834                                                         best_eidx_src = (int)ml_src->e;
 835                                                 }
 836                                         }
 837                                         if (best_eidx_src >= 0) {
 838                                                 mesh_remap_item_define(r_map, i, hit_dist, 0, 1, &best_eidx_src, &full_weight);
 839                                         }
 840                                 }
 841                                 else {
 842                                         /* No source for this dest edge! */
 843                                         BKE_mesh_remap_item_define_invalid(r_map, i);
 844                                 }
 845                         }
 846
 847                         MEM_freeN(vcos_src);
 848                 }
 849                 else if (mode == MREMAP_MODE_EDGE_EDGEINTERP_VNORPROJ) {
 850                         const int num_rays_min = 5, num_rays_max = 100;
 851                         const int numedges_src = dm_src->getNumEdges(dm_src);
 852
 853                         /* Subtleness - this one we can allocate only max number of cast rays per edges! */
 854                         int *indices = MEM_mallocN(sizeof(*indices) * (size_t)min_ii(numedges_src, num_rays_max), __func__);
 855                         /* Here it's simpler to just allocate for all edges :/ */
 856                         float *weights = MEM_mallocN(sizeof(*weights) * (size_t)numedges_src, __func__);
 857
 858                         bvhtree_from_mesh_get(&treedata, dm_src, BVHTREE_FROM_EDGES, 2);
 859
 860                         for (i = 0; i < numedges_dst; i++) {
 861                                 /* For each dst edge, we sample some rays from it (interpolated from its vertices)
 862                                  * and use their hits to interpolate from source edges. */
 863                                 const MEdge *me = &edges_dst[i];
 864                                 float v1_co[3], v2_co[3];
 865                                 float v1_no[3], v2_no[3];
 866
 867                                 int grid_size;
 868                                 float edge_dst_len;
 869                                 float grid_step;
 870
 871                                 float totweights = 0.0f;
 872                                 float hit_dist_accum = 0.0f;
 873                                 int sources_num = 0;
 874                                 int j;
 875
 876                                 copy_v3_v3(v1_co, verts_dst[me->v1].co);
 877                                 copy_v3_v3(v2_co, verts_dst[me->v2].co);
 878
 879                                 normal_short_to_float_v3(v1_no, verts_dst[me->v1].no);
 880                                 normal_short_to_float_v3(v2_no, verts_dst[me->v2].no);
 881
 882                                 /* We do our transform here, allows to interpolate from normals already in src space. */
 883                                 if (space_transform) {
 884                                         BLI_space_transform_apply(space_transform, v1_co);
 885                                         BLI_space_transform_apply(space_transform, v2_co);
 886                                         BLI_space_transform_apply_normal(space_transform, v1_no);
 887                                         BLI_space_transform_apply_normal(space_transform, v2_no);
 888                                 }
 889
 890                                 copy_vn_fl(weights, (int)numedges_src, 0.0f);
 891
 892                                 /* We adjust our ray-casting grid to ray_radius (the smaller, the more rays are cast),
 893                                  * with lower/upper bounds. */
 894                                 edge_dst_len = len_v3v3(v1_co, v2_co);
 895
 896                                 grid_size = (int)((edge_dst_len / ray_radius) + 0.5f);
 897                                 CLAMP(grid_size, num_rays_min, num_rays_max);  /* min 5 rays/edge, max 100. */
 898
 899                                 grid_step = 1.0f / (float)grid_size;  /* Not actual distance here, rather an interp fac... */
 900
 901                                 /* And now we can cast all our rays, and see what we get! */
 902                                 for (j = 0; j < grid_size; j++) {
 903                                         const float fac = grid_step * (float)j;
 904
 905                                         int n = (ray_radius > 0.0f) ? MREMAP_RAYCAST_APPROXIMATE_NR : 1;
 906                                         float w = 1.0f;
 907
 908                                         interp_v3_v3v3(tmp_co, v1_co, v2_co, fac);
 909                                         interp_v3_v3v3_slerp_safe(tmp_no, v1_no, v2_no, fac);
 910
 911                                         while (n--) {
 912                                                 float radius = (ray_radius / w);
 913                                                 treedata.sphere_radius = radius;
 914                                                 if (mesh_remap_bvhtree_query_raycast(
 915                                                         &treedata, &rayhit, tmp_co, tmp_no, radius, max_dist, &hit_dist))
 916                                                 {
 917                                                         weights[rayhit.index] += w;
 918                                                         totweights += w;
 919                                                         hit_dist_accum += hit_dist;
 920                                                         break;
 921                                                 }
 922                                                 /* Next iteration will get bigger radius but smaller weight! */
 923                                                 w /= MREMAP_RAYCAST_APPROXIMATE_FAC;
 924                                         }
 925                                 }
 926                                 /* A sampling is valid (as in, its result can be considered as valid sources) only if at least
 927                                  * half of the rays found a source! */
 928                                 if (totweights > ((float)grid_size / 2.0f)) {
 929                                         for (j = 0; j < (int)numedges_src; j++) {
 930                                                 if (!weights[j]) {
 931                                                         continue;
 932                                                 }
 933                                                 /* Note: sources_num is always <= j! */
 934                                                 weights[sources_num] = weights[j] / totweights;
 935                                                 indices[sources_num] = j;
 936                                                 sources_num++;
 937                                         }
 938                                         mesh_remap_item_define(r_map, i, hit_dist_accum / totweights, 0,
 939                                                                           sources_num, indices, weights);
 940                                 }
 941                                 else {
 942                                         /* No source for this dest edge! */
 943                                         BKE_mesh_remap_item_define_invalid(r_map, i);
 944                                 }
 945                         }
 946
 947                         MEM_freeN(indices);
 948                         MEM_freeN(weights);
 949                 }
 950                 else {
 951                         printf("WARNING! Unsupported mesh-to-mesh edge mapping mode (%d)!\n", mode);
 952                         memset(r_map->items, 0, sizeof(*r_map->items) * (size_t)numedges_dst);
 953                 }
 954
 955                 free_bvhtree_from_mesh(&treedata);
 956         }
 957 }
 958
 959 #define POLY_UNSET       0
 960 #define POLY_CENTER_INIT 1
 961 #define POLY_COMPLETE    2
 962
 963 static void mesh_island_to_astar_graph_edge_process(
 964         MeshIslandStore *islands, const int island_index, BLI_AStarGraph *as_graph,
 965         MVert *verts, MPoly *polys, MLoop *loops,
 966         const int edge_idx, BLI_bitmap *done_edges, MeshElemMap *edge_to_poly_map, const bool is_edge_innercut,
 967         int *poly_island_index_map, float (*poly_centers)[3], unsigned char *poly_status)
 968 {
 969         int *poly_island_indices = BLI_array_alloca(poly_island_indices, (size_t)edge_to_poly_map[edge_idx].count);
 970         int i, j;
 971
 972         for (i = 0; i < edge_to_poly_map[edge_idx].count; i++) {
 973                 const int pidx = edge_to_poly_map[edge_idx].indices[i];
 974                 MPoly *mp = &polys[pidx];
 975                 const int pidx_isld = islands ? poly_island_index_map[pidx] : pidx;
 976                 void *custom_data = is_edge_innercut ? SET_INT_IN_POINTER(edge_idx) : SET_INT_IN_POINTER(-1);
 977
 978                 if (UNLIKELY(islands && (islands->items_to_islands[mp->loopstart] != island_index))) {
 979                         /* poly not in current island, happens with border edges... */
 980                         poly_island_indices[i] = -1;
 981                         continue;
 982                 }
 983
 984                 if (poly_status[pidx_isld] == POLY_COMPLETE) {
 985                         poly_island_indices[i] = pidx_isld;
 986                         continue;
 987                 }
 988
 989                 if (poly_status[pidx_isld] == POLY_UNSET) {
 990                         BKE_mesh_calc_poly_center(mp, &loops[mp->loopstart], verts, poly_centers[pidx_isld]);
 991                         BLI_astar_node_init(as_graph, pidx_isld, poly_centers[pidx_isld]);
 992                         poly_status[pidx_isld] = POLY_CENTER_INIT;
 993                 }
 994
 995                 for (j = i; j--;) {
 996                         float dist_cost;
 997                         const int pidx_isld_other = poly_island_indices[j];
 998
 999                         if (pidx_isld_other == -1 || poly_status[pidx_isld_other] == POLY_COMPLETE) {
1000                                 /* If the other poly is complete, that link has already been added! */
1001                                 continue;
1002                         }
1003                         dist_cost = len_v3v3(poly_centers[pidx_isld_other], poly_centers[pidx_isld]);
1004                         BLI_astar_node_link_add(as_graph, pidx_isld_other, pidx_isld, dist_cost, custom_data);
1005                 }
1006
1007                 poly_island_indices[i] = pidx_isld;
1008         }
1009
1010         BLI_BITMAP_ENABLE(done_edges, edge_idx);
1011 }
1012
1013 static void mesh_island_to_astar_graph(
1014         MeshIslandStore *islands, const int island_index,
1015         MVert *verts, MeshElemMap *edge_to_poly_map, const int numedges, MLoop *loops, MPoly *polys, const int numpolys,
1016         BLI_AStarGraph *r_as_graph)
1017 {
1018         MeshElemMap *island_poly_map = islands ? islands->islands[island_index] : NULL;
1019         MeshElemMap *island_einnercut_map = islands ? islands->innercuts[island_index] : NULL;
1020
1021         int *poly_island_index_map = NULL;
1022         BLI_bitmap *done_edges = BLI_BITMAP_NEW(numedges, __func__);
1023
1024         const int node_num = islands ? island_poly_map->count : numpolys;
1025         unsigned char *poly_status = MEM_callocN(sizeof(*poly_status) * (size_t)node_num, __func__);
1026         float (*poly_centers)[3];
1027
1028         int pidx_isld;
1029         int i;
1030
1031         BLI_astar_graph_init(r_as_graph, node_num, NULL);
1032         /* poly_centers is owned by graph memarena. */
1033         poly_centers = BLI_memarena_calloc(r_as_graph->mem, sizeof(*poly_centers) * (size_t)node_num);
1034
1035         if (islands) {
1036                 /* poly_island_index_map is owned by graph memarena. */
1037                 poly_island_index_map = BLI_memarena_calloc(r_as_graph->mem, sizeof(*poly_island_index_map) * (size_t)numpolys);
1038                 for (i = island_poly_map->count; i--;) {
1039                         poly_island_index_map[island_poly_map->indices[i]] = i;
1040                 }
1041
1042                 r_as_graph->custom_data = poly_island_index_map;
1043
1044                 for (i = island_einnercut_map->count; i--;) {
1045                         mesh_island_to_astar_graph_edge_process(
1046                                 islands, island_index, r_as_graph, verts, polys, loops,
1047                                 island_einnercut_map->indices[i], done_edges, edge_to_poly_map, true,
1048                                 poly_island_index_map, poly_centers, poly_status);
1049                 }
1050         }
1051
1052         for (pidx_isld = node_num; pidx_isld--;) {
1053                 const int pidx = islands ? island_poly_map->indices[pidx_isld] : pidx_isld;
1054                 MPoly *mp = &polys[pidx];
1055                 int pl_idx, l_idx;
1056
1057                 if (poly_status[pidx_isld] == POLY_COMPLETE) {
1058                         continue;
1059                 }
1060
1061                 for (pl_idx = 0, l_idx = mp->loopstart; pl_idx < mp->totloop; pl_idx++, l_idx++) {
1062                         MLoop *ml = &loops[l_idx];
1063
1064                         if (BLI_BITMAP_TEST(done_edges, ml->e)) {
1065                                 continue;
1066                         }
1067
1068                         mesh_island_to_astar_graph_edge_process(
1069                                 islands, island_index, r_as_graph, verts, polys, loops,
1070                                 (int)ml->e, done_edges, edge_to_poly_map, false,
1071                                 poly_island_index_map, poly_centers, poly_status);
1072                 }
1073                 poly_status[pidx_isld] = POLY_COMPLETE;
1074         }
1075
1076         MEM_freeN(done_edges);
1077         MEM_freeN(poly_status);
1078 }
1079
1080 #undef POLY_UNSET
1081 #undef POLY_CENTER_INIT
1082 #undef POLY_COMPLETE
1083
1084 /* Our 'f_cost' callback func, to find shortest poly-path between two remapped-loops.
1085  * Note we do not want to make innercuts 'walls' here, just detect when the shortest path goes by those. */
1086 static float mesh_remap_calc_loops_astar_f_cost(
1087         BLI_AStarGraph *as_graph, BLI_AStarSolution *as_solution, BLI_AStarGNLink *link,
1088         const int node_idx_curr, const int node_idx_next, const int node_idx_dst)
1089 {
1090         float *co_next, *co_dest;
1091
1092         if (link && (GET_INT_FROM_POINTER(link->custom_data) != -1)) {
1093                 /* An innercut edge... We tag our solution as potentially crossing innercuts.
1094                  * Note it might not be the case in the end (AStar will explore around optimal path), but helps
1095                  * trimming off some processing later... */
1096                 if (!GET_INT_FROM_POINTER(as_solution->custom_data)) {
1097                         as_solution->custom_data = SET_INT_IN_POINTER(true);
1098                 }
1099         }
1100
1101         /* Our heuristic part of current f_cost is distance from next node to destination one.
1102          * It is guaranteed to be less than (or equal to) actual shortest poly-path between next node and destination one.
1103          */
1104         co_next = (float *)as_graph->nodes[node_idx_next].custom_data;
1105         co_dest = (float *)as_graph->nodes[node_idx_dst].custom_data;
1106         return (link ? (as_solution->g_costs[node_idx_curr] + link->cost) : 0.0f) + len_v3v3(co_next, co_dest);
1107 }
1108
1109 #define ASTAR_STEPS_MAX 64
1110
1111
1112 void BKE_mesh_remap_calc_loops_from_dm(
1113         const int mode, const SpaceTransform *space_transform, const float max_dist, const float ray_radius,
1114         MVert *verts_dst, const int numverts_dst, MEdge *edges_dst, const int numedges_dst,
1115         MLoop *loops_dst, const int numloops_dst, MPoly *polys_dst, const int numpolys_dst,
1116         CustomData *ldata_dst, CustomData *pdata_dst,
1117         const bool use_split_nors_dst, const float split_angle_dst, const bool dirty_nors_dst,
1118         DerivedMesh *dm_src, const bool use_split_nors_src, const float split_angle_src,
1119         MeshRemapIslandsCalc gen_islands_src, const float islands_precision_src, MeshPairRemap *r_map)
1120 {
1121         const float full_weight = 1.0f;
1122         const float max_dist_sq = max_dist * max_dist;
1123
1124         int i;
1125
1126         BLI_assert(mode & MREMAP_MODE_LOOP);
1127         BLI_assert((islands_precision_src >= 0.0f) && (islands_precision_src <= 1.0f));
1128
1129         BKE_mesh_remap_init(r_map, numloops_dst);
1130
1131         if (mode == MREMAP_MODE_TOPOLOGY) {
1132                 /* In topology mapping, we assume meshes are identical, islands included! */
1133                 BLI_assert(numloops_dst == dm_src->getNumLoops(dm_src));
1134                 for (i = 0; i < numloops_dst; i++) {
1135                         mesh_remap_item_define(r_map, i, FLT_MAX, 0, 1, &i, &full_weight);
1136                 }
1137         }
1138         else {
1139                 BVHTreeFromMesh *treedata = NULL;
1140                 BVHTreeNearest nearest = {0};
1141                 BVHTreeRayHit rayhit = {0};
1142                 int num_trees = 0;
1143                 float hit_dist;
1144                 float tmp_co[3], tmp_no[3];
1145
1146                 const bool use_from_vert = (mode & MREMAP_USE_VERT);
1147
1148                 MeshIslandStore island_store = {0};
1149                 bool use_islands = false;
1150
1151                 BLI_AStarGraph *as_graphdata = NULL;
1152                 BLI_AStarSolution as_solution = {0};
1153                 const int isld_steps_src = islands_precision_src ?
1154                                            max_ii((int)(ASTAR_STEPS_MAX * islands_precision_src + 0.499f), 1) : 0;
1155
1156                 float (*poly_nors_src)[3] = NULL;
1157                 float (*loop_nors_src)[3] = NULL;
1158                 float (*poly_nors_dst)[3] = NULL;
1159                 float (*loop_nors_dst)[3] = NULL;
1160
1161                 float (*poly_cents_src)[3] = NULL;
1162
1163                 MeshElemMap *vert_to_loop_map_src = NULL;
1164                 int *vert_to_loop_map_src_buff = NULL;
1165                 MeshElemMap *vert_to_poly_map_src = NULL;
1166                 int *vert_to_poly_map_src_buff = NULL;
1167                 MeshElemMap *edge_to_poly_map_src = NULL;
1168                 int *edge_to_poly_map_src_buff = NULL;
1169                 MeshElemMap *poly_to_looptri_map_src = NULL;
1170                 int *poly_to_looptri_map_src_buff = NULL;
1171
1172                 /* Unlike above, those are one-to-one mappings, simpler! */
1173                 int *loop_to_poly_map_src = NULL;
1174
1175                 bool verts_allocated_src;
1176                 MVert *verts_src = DM_get_vert_array(dm_src, &verts_allocated_src);
1177                 const int num_verts_src = dm_src->getNumVerts(dm_src);
1178                 float (*vcos_src)[3] = NULL;
1179                 bool edges_allocated_src;
1180                 MEdge *edges_src = DM_get_edge_array(dm_src, &edges_allocated_src);
1181                 const int num_edges_src = dm_src->getNumEdges(dm_src);
1182                 bool loops_allocated_src;
1183                 MLoop *loops_src = DM_get_loop_array(dm_src, &loops_allocated_src);
1184                 const int num_loops_src = dm_src->getNumLoops(dm_src);
1185                 bool polys_allocated_src;
1186                 MPoly *polys_src = DM_get_poly_array(dm_src, &polys_allocated_src);
1187                 const int num_polys_src = dm_src->getNumPolys(dm_src);
1188                 const MLoopTri *looptri_src = NULL;
1189                 int num_looptri_src = 0;
1190
1191                 size_t buff_size_interp = MREMAP_DEFAULT_BUFSIZE;
1192                 float (*vcos_interp)[3] = NULL;
1193                 int *indices_interp = NULL;
1194                 float *weights_interp = NULL;
1195
1196                 MLoop *ml_src, *ml_dst;
1197                 MPoly *mp_src, *mp_dst;
1198                 int tindex, pidx_dst, lidx_dst, plidx_dst, pidx_src, lidx_src, plidx_src;
1199
1200                 IslandResult **islands_res;
1201                 size_t islands_res_buff_size = MREMAP_DEFAULT_BUFSIZE;
1202
1203                 if (!use_from_vert) {
1204                         vcos_src = MEM_mallocN(sizeof(*vcos_src) * (size_t)num_verts_src, __func__);
1205                         dm_src->getVertCos(dm_src, vcos_src);
1206
1207                         vcos_interp = MEM_mallocN(sizeof(*vcos_interp) * buff_size_interp, __func__);
1208                         indices_interp = MEM_mallocN(sizeof(*indices_interp) * buff_size_interp, __func__);
1209                         weights_interp = MEM_mallocN(sizeof(*weights_interp) * buff_size_interp, __func__);
1210                 }
1211
1212                 {
1213                         const bool need_lnors_src = (mode & MREMAP_USE_LOOP) && (mode & MREMAP_USE_NORMAL);
1214                         const bool need_lnors_dst = need_lnors_src || (mode & MREMAP_USE_NORPROJ);
1215                         const bool need_pnors_src = need_lnors_src || ((mode & MREMAP_USE_POLY) && (mode & MREMAP_USE_NORMAL));
1216                         const bool need_pnors_dst = need_lnors_dst || need_pnors_src;
1217
1218                         if (need_pnors_dst) {
1219                                 /* Cache poly nors into a temp CDLayer. */
1220                                 poly_nors_dst = CustomData_get_layer(pdata_dst, CD_NORMAL);
1221                                 if (!poly_nors_dst) {
1222                                         poly_nors_dst = CustomData_add_layer(pdata_dst, CD_NORMAL, CD_CALLOC, NULL, numpolys_dst);
1223                                         CustomData_set_layer_flag(pdata_dst, CD_NORMAL, CD_FLAG_TEMPORARY);
1224                                 }
1225                                 if (dirty_nors_dst) {
1226                                         BKE_mesh_calc_normals_poly(verts_dst, NULL, numverts_dst, loops_dst, polys_dst,
1227                                                                    numloops_dst, numpolys_dst, poly_nors_dst, true);
1228                                 }
1229                         }
1230                         if (need_lnors_dst) {
1231                                 short (*custom_nors_dst)[2] = CustomData_get_layer(ldata_dst, CD_CUSTOMLOOPNORMAL);
1232
1233                                 /* Cache poly nors into a temp CDLayer. */
1234                                 loop_nors_dst = CustomData_get_layer(ldata_dst, CD_NORMAL);
1235                                 if (dirty_nors_dst || !loop_nors_dst) {
1236                                         if (!loop_nors_dst) {
1237                                                 loop_nors_dst = CustomData_add_layer(ldata_dst, CD_NORMAL, CD_CALLOC, NULL, numloops_dst);
1238                                                 CustomData_set_layer_flag(ldata_dst, CD_NORMAL, CD_FLAG_TEMPORARY);
1239                                         }
1240                                         BKE_mesh_normals_loop_split(verts_dst, numverts_dst, edges_dst, numedges_dst,
1241                                                                     loops_dst, loop_nors_dst, numloops_dst,
1242                                                                     polys_dst, (const float (*)[3])poly_nors_dst, numpolys_dst,
1243                                                                     use_split_nors_dst, split_angle_dst, NULL, custom_nors_dst, NULL);
1244                                 }
1245                         }
1246                         if (need_pnors_src || need_lnors_src) {
1247                                 /* Simpler for now, calcNormals never stores pnors :( */
1248                                 dm_src->calcLoopNormals(dm_src, use_split_nors_src, split_angle_src);
1249
1250                                 if (need_pnors_src) {
1251                                         poly_nors_src = dm_src->getPolyDataArray(dm_src, CD_NORMAL);
1252                                 }
1253                                 if (need_lnors_src) {
1254                                         loop_nors_src = dm_src->getLoopDataArray(dm_src, CD_NORMAL);
1255                                 }
1256                         }
1257                 }
1258
1259                 if (use_from_vert) {
1260                         BKE_mesh_vert_loop_map_create(&vert_to_loop_map_src, &vert_to_loop_map_src_buff,
1261                                                       polys_src, loops_src, num_verts_src, num_polys_src, num_loops_src);
1262                         if (mode & MREMAP_USE_POLY) {
1263                                 BKE_mesh_vert_poly_map_create(&vert_to_poly_map_src, &vert_to_poly_map_src_buff,
1264                                                               polys_src, loops_src, num_verts_src, num_polys_src, num_loops_src);
1265                         }
1266                 }
1267
1268                 /* Needed for islands (or plain mesh) to AStar graph conversion. */
1269                 BKE_mesh_edge_poly_map_create(&edge_to_poly_map_src, &edge_to_poly_map_src_buff,
1270                                               edges_src, num_edges_src, polys_src, num_polys_src, loops_src, num_loops_src);
1271                 if (use_from_vert) {
1272                         loop_to_poly_map_src = MEM_mallocN(sizeof(*loop_to_poly_map_src) * (size_t)num_loops_src, __func__);
1273                         poly_cents_src = MEM_mallocN(sizeof(*poly_cents_src) * (size_t)num_polys_src, __func__);
1274                         for (pidx_src = 0, mp_src = polys_src; pidx_src < num_polys_src; pidx_src++, mp_src++) {
1275                                 ml_src = &loops_src[mp_src->loopstart];
1276                                 for (plidx_src = 0, lidx_src = mp_src->loopstart; plidx_src < mp_src->totloop; plidx_src++, lidx_src++) {
1277                                         loop_to_poly_map_src[lidx_src] = pidx_src;
1278                                 }
1279                                 BKE_mesh_calc_poly_center(mp_src, ml_src, verts_src, poly_cents_src[pidx_src]);
1280                         }
1281                 }
1282
1283                 /* Island makes things slightly more complex here.
1284                  * Basically, we:
1285                  *     * Make one treedata for each island's elements.
1286                  *     * Check all loops of a same dest poly against all treedata.
1287                  *     * Choose the island's elements giving the best results.
1288                  */
1289
1290                 /* First, generate the islands, if possible. */
1291                 if (gen_islands_src) {
1292                         use_islands = gen_islands_src(
1293                                 verts_src, num_verts_src,
1294                                 edges_src, num_edges_src,
1295                                 polys_src, num_polys_src,
1296                                 loops_src, num_loops_src,
1297                                 &island_store);
1298
1299                         num_trees = use_islands ? island_store.islands_num : 1;
1300                         treedata = MEM_callocN(sizeof(*treedata) * (size_t)num_trees, __func__);
1301                         if (isld_steps_src) {
1302                                 as_graphdata = MEM_callocN(sizeof(*as_graphdata) * (size_t)num_trees, __func__);
1303                         }
1304
1305                         if (use_islands) {
1306                                 /* We expect our islands to contain poly indices, with edge indices of 'inner cuts',
1307                                  * and a mapping loops -> islands indices.
1308                                  * This implies all loops of a same poly are in the same island. */
1309                                 BLI_assert((island_store.item_type == MISLAND_TYPE_LOOP) &&
1310                                            (island_store.island_type == MISLAND_TYPE_POLY) &&
1311                                            (island_store.innercut_type == MISLAND_TYPE_EDGE));
1312                         }
1313                 }
1314                 else {
1315                         num_trees = 1;
1316                         treedata = MEM_callocN(sizeof(*treedata), __func__);
1317                         if (isld_steps_src) {
1318                                 as_graphdata = MEM_callocN(sizeof(*as_graphdata), __func__);
1319                         }
1320                 }
1321
1322                 /* Build our AStar graphs. */
1323                 if (isld_steps_src) {
1324                         for (tindex = 0; tindex < num_trees; tindex++) {
1325                                 mesh_island_to_astar_graph(
1326                                         use_islands ? &island_store : NULL, tindex, verts_src, edge_to_poly_map_src, num_edges_src,
1327                                         loops_src, polys_src, num_polys_src, &as_graphdata[tindex]);
1328                         }
1329                 }
1330
1331                 /* Build our BVHtrees, either from verts or tessfaces. */
1332                 if (use_from_vert) {
1333                         if (use_islands) {
1334                                 BLI_bitmap *verts_active = BLI_BITMAP_NEW((size_t)num_verts_src, __func__);
1335
1336                                 for (tindex = 0; tindex < num_trees; tindex++) {
1337                                         MeshElemMap *isld = island_store.islands[tindex];
1338                                         int num_verts_active = 0;
1339                                         BLI_BITMAP_SET_ALL(verts_active, false, (size_t)num_verts_src);
1340                                         for (i = 0; i < isld->count; i++) {
1341                                                 mp_src = &polys_src[isld->indices[i]];
1342                                                 for (lidx_src = mp_src->loopstart; lidx_src < mp_src->loopstart + mp_src->totloop; lidx_src++) {
1343                                                         const unsigned int vidx_src = loops_src[lidx_src].v;
1344                                                         if (!BLI_BITMAP_TEST(verts_active, vidx_src)) {
1345                                                                 BLI_BITMAP_ENABLE(verts_active, loops_src[lidx_src].v);
1346                                                                 num_verts_active++;
1347                                                         }
1348                                                 }
1349                                         }
1350                                         /* verts 'ownership' is transfered to treedata here, which will handle its freeing. */
1351                                         bvhtree_from_mesh_verts_ex(&treedata[tindex], verts_src, num_verts_src, verts_allocated_src,
1352                                                                    verts_active, num_verts_active, 0.0, 2, 6);
1353                                         if (verts_allocated_src) {
1354                                                 verts_allocated_src = false;  /* Only 'give' our verts once, to first tree! */
1355                                         }
1356                                 }
1357
1358                                 MEM_freeN(verts_active);
1359                         }
1360                         else {
1361                                 BLI_assert(num_trees == 1);
1362                                 bvhtree_from_mesh_get(&treedata[0], dm_src, BVHTREE_FROM_VERTS, 2);
1363                         }
1364                 }
1365                 else {  /* We use polygons. */
1366                         if (use_islands) {
1367                                 /* bvhtree here uses looptri faces... */
1368                                 const unsigned int dirty_tess_flag = dm_src->dirty & DM_DIRTY_TESS_CDLAYERS;
1369                                 BLI_bitmap *looptri_active;
1370
1371                                 /* We do not care about tessellated data here, only geometry itself is important. */
1372                                 if (dirty_tess_flag) {
1373                                         dm_src->dirty &= ~dirty_tess_flag;
1374                                 }
1375                                 if (dirty_tess_flag) {
1376                                         dm_src->dirty |= dirty_tess_flag;
1377                                 }
1378
1379                                 looptri_src = dm_src->getLoopTriArray(dm_src);
1380                                 num_looptri_src = dm_src->getNumLoopTri(dm_src);
1381                                 looptri_active = BLI_BITMAP_NEW((size_t)num_looptri_src, __func__);
1382
1383                                 for (tindex = 0; tindex < num_trees; tindex++) {
1384                                         int num_looptri_active = 0;
1385                                         BLI_BITMAP_SET_ALL(looptri_active, false, (size_t)num_looptri_src);
1386                                         for (i = 0; i < num_looptri_src; i++) {
1387                                                 mp_src = &polys_src[looptri_src[i].poly];
1388                                                 if (island_store.items_to_islands[mp_src->loopstart] == tindex) {
1389                                                         BLI_BITMAP_ENABLE(looptri_active, i);
1390                                                         num_looptri_active++;
1391                                                 }
1392                                         }
1393                                         /* verts and faces 'ownership' is transfered to treedata here, which will handle its freeing. */
1394                                         bvhtree_from_mesh_looptri_ex(
1395                                                 &treedata[tindex],
1396                                                 verts_src, verts_allocated_src,
1397                                                 loops_src, loops_allocated_src,
1398                                                 looptri_src, num_looptri_src, false,
1399                                                 looptri_active, num_looptri_active, 0.0, 2, 6);
1400                                         if (verts_allocated_src) {
1401                                                 verts_allocated_src = false;  /* Only 'give' our verts once, to first tree! */
1402                                         }
1403                                         if (loops_allocated_src) {
1404                                                 loops_allocated_src = false;  /* Only 'give' our loops once, to first tree! */
1405                                         }
1406                                 }
1407
1408                                 MEM_freeN(looptri_active);
1409                         }
1410                         else {
1411                                 BLI_assert(num_trees == 1);
1412                                 bvhtree_from_mesh_looptri(&treedata[0], dm_src, 0.0, 2, 6);
1413                         }
1414                 }
1415
1416                 /* And check each dest poly! */
1417                 islands_res = MEM_mallocN(sizeof(*islands_res) * (size_t)num_trees, __func__);
1418                 for (tindex = 0; tindex < num_trees; tindex++) {
1419                         islands_res[tindex] = MEM_mallocN(sizeof(**islands_res) * islands_res_buff_size, __func__);
1420                 }
1421
1422                 for (pidx_dst = 0, mp_dst = polys_dst; pidx_dst < numpolys_dst; pidx_dst++, mp_dst++) {
1423                         float pnor_dst[3];
1424
1425                         /* Only in use_from_vert case, we may need polys' centers as fallback in case we cannot decide which
1426                          * corner to use from normals only. */
1427                         float pcent_dst[3];
1428                         bool pcent_dst_valid = false;
1429
1430                         if (mode == MREMAP_MODE_LOOP_NEAREST_POLYNOR) {
1431                                 copy_v3_v3(pnor_dst, poly_nors_dst[pidx_dst]);
1432                                 if (space_transform) {
1433                                         BLI_space_transform_apply_normal(space_transform, pnor_dst);
1434                                 }
1435                         }
1436
1437                         if ((size_t)mp_dst->totloop > islands_res_buff_size) {
1438                                 islands_res_buff_size = (size_t)mp_dst->totloop + MREMAP_DEFAULT_BUFSIZE;
1439                                 for (tindex = 0; tindex < num_trees; tindex++) {
1440                                         islands_res[tindex] = MEM_reallocN(islands_res[tindex], sizeof(**islands_res) * islands_res_buff_size);
1441                                 }
1442                         }
1443
1444                         for (tindex = 0; tindex < num_trees; tindex++) {
1445                                 BVHTreeFromMesh *tdata = &treedata[tindex];
1446
1447                                 ml_dst = &loops_dst[mp_dst->loopstart];
1448                                 for (plidx_dst = 0; plidx_dst < mp_dst->totloop; plidx_dst++, ml_dst++) {
1449                                         if (use_from_vert) {
1450                                                 MeshElemMap *vert_to_refelem_map_src = NULL;
1451
1452                                                 copy_v3_v3(tmp_co, verts_dst[ml_dst->v].co);
1453                                                 nearest.index = -1;
1454
1455                                                 /* Convert the vertex to tree coordinates, if needed. */
1456                                                 if (space_transform) {
1457                                                         BLI_space_transform_apply(space_transform, tmp_co);
1458                                                 }
1459
1460                                                 if (mesh_remap_bvhtree_query_nearest(tdata, &nearest, tmp_co, max_dist_sq, &hit_dist)) {
1461                                                         float (*nor_dst)[3];
1462                                                         float (*nors_src)[3];
1463                                                         float best_nor_dot = -2.0f;
1464                                                         float best_sqdist_fallback = FLT_MAX;
1465                                                         int best_index_src = -1;
1466
1467                                                         if (mode == MREMAP_MODE_LOOP_NEAREST_LOOPNOR) {
1468                                                                 copy_v3_v3(tmp_no, loop_nors_dst[plidx_dst + mp_dst->loopstart]);
1469                                                                 if (space_transform) {
1470                                                                         BLI_space_transform_apply_normal(space_transform, tmp_no);
1471                                                                 }
1472                                                                 nor_dst = &tmp_no;
1473                                                                 nors_src = loop_nors_src;
1474                                                                 vert_to_refelem_map_src = vert_to_loop_map_src;
1475                                                         }
1476                                                         else {  /* if (mode == MREMAP_MODE_LOOP_NEAREST_POLYNOR) { */
1477                                                                 nor_dst = &pnor_dst;
1478                                                                 nors_src = poly_nors_src;
1479                                                                 vert_to_refelem_map_src = vert_to_poly_map_src;
1480                                                         }
1481
1482                                                         for (i = vert_to_refelem_map_src[nearest.index].count; i--;) {
1483                                                                 const int index_src = vert_to_refelem_map_src[nearest.index].indices[i];
1484                                                                 BLI_assert(index_src != -1);
1485                                                                 const float dot = dot_v3v3(nors_src[index_src], *nor_dst);
1486
1487                                                                 pidx_src = (mode == MREMAP_MODE_LOOP_NEAREST_LOOPNOR) ?
1488                                                                                loop_to_poly_map_src[index_src] : index_src;
1489                                                                 /* WARNING! This is not the *real* lidx_src in case of POLYNOR, we only use it
1490                                                                  *          to check we stay on current island (all loops from a given poly are
1491                                                                  *          on same island!). */
1492                                                                 lidx_src = (mode == MREMAP_MODE_LOOP_NEAREST_LOOPNOR) ?
1493                                                                                index_src : polys_src[pidx_src].loopstart;
1494
1495                                                                 /* A same vert may be at the boundary of several islands! Hence, we have to ensure
1496                                                                  * poly/loop we are currently considering *belongs* to current island! */
1497                                                                 if (use_islands && island_store.items_to_islands[lidx_src] != tindex) {
1498                                                                         continue;
1499                                                                 }
1500
1501                                                                 if (dot > best_nor_dot - 1e-6f) {
1502                                                                         /* We need something as fallback decision in case dest normal matches several
1503                                                                          * source normals (see T44522), using distance between polys' centers here. */
1504                                                                         float *pcent_src;
1505                                                                         float sqdist;
1506
1507                                                                         mp_src = &polys_src[pidx_src];
1508                                                                         ml_src = &loops_src[mp_src->loopstart];
1509
1510                                                                         if (!pcent_dst_valid) {
1511                                                                                 BKE_mesh_calc_poly_center(
1512                                                                                             mp_dst, &loops_dst[mp_dst->loopstart], verts_dst, pcent_dst);
1513                                                                                 pcent_dst_valid = true;
1514                                                                         }
1515                                                                         pcent_src = poly_cents_src[pidx_src];
1516                                                                         sqdist = len_squared_v3v3(pcent_dst, pcent_src);
1517
1518                                                                         if ((dot > best_nor_dot + 1e-6f) || (sqdist < best_sqdist_fallback)) {
1519                                                                                 best_nor_dot = dot;
1520                                                                                 best_sqdist_fallback = sqdist;
1521                                                                                 best_index_src = index_src;
1522                                                                         }
1523                                                                 }
1524                                                         }
1525                                                         if (best_index_src == -1) {
1526                                                                 /* We found no item to map back from closest vertex... */
1527                                                                 best_nor_dot = -1.0f;
1528                                                                 hit_dist = FLT_MAX;
1529                                                         }
1530                                                         else if (mode == MREMAP_MODE_LOOP_NEAREST_POLYNOR) {
1531                                                                 /* Our best_index_src is a poly one for now!
1532                                                                  * Have to find its loop matching our closest vertex. */
1533                                                                 mp_src = &polys_src[best_index_src];
1534                                                                 ml_src = &loops_src[mp_src->loopstart];
1535                                                                 for (plidx_src = 0; plidx_src < mp_src->totloop; plidx_src++, ml_src++) {
1536                                                                         if ((int)ml_src->v == nearest.index) {
1537                                                                                 best_index_src = plidx_src + mp_src->loopstart;
1538                                                                                 break;
1539                                                                         }
1540                                                                 }
1541                                                         }
1542                                                         best_nor_dot = (best_nor_dot + 1.0f) * 0.5f;
1543                                                         islands_res[tindex][plidx_dst].factor = hit_dist ? (best_nor_dot / hit_dist) : 1e18f;
1544                                                         islands_res[tindex][plidx_dst].hit_dist = hit_dist;
1545                                                         islands_res[tindex][plidx_dst].index_src = best_index_src;
1546                                                 }
1547                                                 else {
1548                                                         /* No source for this dest loop! */
1549                                                         islands_res[tindex][plidx_dst].factor = 0.0f;
1550                                                         islands_res[tindex][plidx_dst].hit_dist = FLT_MAX;
1551                                                         islands_res[tindex][plidx_dst].index_src = -1;
1552                                                 }
1553                                         }
1554                                         else if (mode & MREMAP_USE_NORPROJ) {
1555                                                 int n = (ray_radius > 0.0f) ? MREMAP_RAYCAST_APPROXIMATE_NR : 1;
1556                                                 float w = 1.0f;
1557
1558                                                 copy_v3_v3(tmp_co, verts_dst[ml_dst->v].co);
1559                                                 copy_v3_v3(tmp_no, loop_nors_dst[plidx_dst + mp_dst->loopstart]);
1560
1561                                                 /* We do our transform here, since we may do several raycast/nearest queries. */
1562                                                 if (space_transform) {
1563                                                         BLI_space_transform_apply(space_transform, tmp_co);
1564                                                         BLI_space_transform_apply_normal(space_transform, tmp_no);
1565                                                 }
1566
1567                                                 while (n--) {
1568                                                         float radius = ray_radius / w;
1569                                                         tdata->sphere_radius = radius;
1570                                                         if (mesh_remap_bvhtree_query_raycast(
1571                                                                 tdata, &rayhit, tmp_co, tmp_no, radius, max_dist, &hit_dist))
1572                                                         {
1573                                                                 islands_res[tindex][plidx_dst].factor = (hit_dist ? (1.0f / hit_dist) : 1e18f) * w;
1574                                                                 islands_res[tindex][plidx_dst].hit_dist = hit_dist;
1575                                                                 islands_res[tindex][plidx_dst].index_src = (int)tdata->looptri[rayhit.index].poly;
1576                                                                 copy_v3_v3(islands_res[tindex][plidx_dst].hit_point, rayhit.co);
1577                                                                 break;
1578                                                         }
1579                                                         /* Next iteration will get bigger radius but smaller weight! */
1580                                                         w /= MREMAP_RAYCAST_APPROXIMATE_FAC;
1581                                                 }
1582                                                 if (n == -1) {
1583                                                         /* Fallback to 'nearest' hit here, loops usually comes in 'face group', not good to
1584                                                          * have only part of one dest face's loops to map to source.
1585                                                          * Note that since we give this a null weight, if whole weight for a given face
1586                                                          * is null, it means none of its loop mapped to this source island, hence we can skip it
1587                                                          * later.
1588                                                          */
1589                                                         copy_v3_v3(tmp_co, verts_dst[ml_dst->v].co);
1590                                                         nearest.index = -1;
1591
1592                                                         /* Convert the vertex to tree coordinates, if needed. */
1593                                                         if (space_transform) {
1594                                                                 BLI_space_transform_apply(space_transform, tmp_co);
1595                                                         }
1596
1597                                                         /* In any case, this fallback nearest hit should have no weight at all
1598                                                          * in 'best island' decision! */
1599                                                         islands_res[tindex][plidx_dst].factor = 0.0f;
1600
1601                                                         if (mesh_remap_bvhtree_query_nearest(tdata, &nearest, tmp_co, max_dist_sq, &hit_dist)) {
1602                                                                 islands_res[tindex][plidx_dst].hit_dist = hit_dist;
1603                                                                 islands_res[tindex][plidx_dst].index_src = (int)tdata->looptri[nearest.index].poly;
1604                                                                 copy_v3_v3(islands_res[tindex][plidx_dst].hit_point, nearest.co);
1605                                                         }
1606                                                         else {
1607                                                                 /* No source for this dest loop! */
1608                                                                 islands_res[tindex][plidx_dst].hit_dist = FLT_MAX;
1609                                                                 islands_res[tindex][plidx_dst].index_src = -1;
1610                                                         }
1611                                                 }
1612                                         }
1613                                         else {  /* Nearest poly either to use all its loops/verts or just closest one. */
1614                                                 copy_v3_v3(tmp_co, verts_dst[ml_dst->v].co);
1615                                                 nearest.index = -1;
1616
1617                                                 /* Convert the vertex to tree coordinates, if needed. */
1618                                                 if (space_transform) {
1619                                                         BLI_space_transform_apply(space_transform, tmp_co);
1620                                                 }
1621
1622                                                 if (mesh_remap_bvhtree_query_nearest(tdata, &nearest, tmp_co, max_dist_sq, &hit_dist)) {
1623                                                         islands_res[tindex][plidx_dst].factor = hit_dist ? (1.0f / hit_dist) : 1e18f;
1624                                                         islands_res[tindex][plidx_dst].hit_dist = hit_dist;
1625                                                         islands_res[tindex][plidx_dst].index_src = (int)tdata->looptri[nearest.index].poly;
1626                                                         copy_v3_v3(islands_res[tindex][plidx_dst].hit_point, nearest.co);
1627                                                 }
1628                                                 else {
1629                                                         /* No source for this dest loop! */
1630                                                         islands_res[tindex][plidx_dst].factor = 0.0f;
1631                                                         islands_res[tindex][plidx_dst].hit_dist = FLT_MAX;
1632                                                         islands_res[tindex][plidx_dst].index_src = -1;
1633                                                 }
1634                                         }
1635                                 }
1636                         }
1637
1638                         /* And now, find best island to use! */
1639                         /* We have to first select the 'best source island' for given dst poly and its loops.
1640                          * Then, we have to check that poly does not 'spread' across some island's limits
1641                          * (like inner seams for UVs, etc.).
1642                          * Note we only still partially support that kind of situation here, i.e. polys spreading over actual cracks
1643                          * (like a narrow space without faces on src, splitting a 'tube-like' geometry). That kind of situation
1644                          * should be relatively rare, though.
1645                          */
1646                         /* XXX This block in itself is big and complex enough to be a separate function but... it uses a bunch
1647                          *     of locale vars. Not worth sending all that through parameters (for now at least). */
1648                         {
1649                                 BLI_AStarGraph *as_graph = NULL;
1650                                 int *poly_island_index_map = NULL;
1651                                 int pidx_src_prev = -1;
1652
1653                                 MeshElemMap *best_island = NULL;
1654                                 float best_island_fac = 0.0f;
1655                                 int best_island_index = -1;
1656
1657                                 for (tindex = 0; tindex < num_trees; tindex++) {
1658                                         float island_fac = 0.0f;
1659
1660                                         for (plidx_dst = 0; plidx_dst < mp_dst->totloop; plidx_dst++) {
1661                                                 island_fac += islands_res[tindex][plidx_dst].factor;
1662                                         }
1663                                         island_fac /= (float)mp_dst->totloop;
1664
1665                                         if (island_fac > best_island_fac) {
1666                                                 best_island_fac = island_fac;
1667                                                 best_island_index = tindex;
1668                                         }
1669                                 }
1670
1671                                 if (best_island_index != -1 && isld_steps_src) {
1672                                         best_island = use_islands ? island_store.islands[best_island_index] : NULL;
1673                                         as_graph = &as_graphdata[best_island_index];
1674                                         poly_island_index_map = (int *)as_graph->custom_data;
1675                                         BLI_astar_solution_init(as_graph, &as_solution, NULL);
1676                                 }
1677
1678                                 for (plidx_dst = 0; plidx_dst < mp_dst->totloop; plidx_dst++) {
1679                                         IslandResult *isld_res;
1680                                         lidx_dst = plidx_dst + mp_dst->loopstart;
1681
1682                                         if (best_island_index == -1) {
1683                                                 /* No source for any loops of our dest poly in any source islands. */
1684                                                 BKE_mesh_remap_item_define_invalid(r_map, lidx_dst);
1685                                                 continue;
1686                                         }
1687
1688                                         as_solution.custom_data = SET_INT_IN_POINTER(false);
1689
1690                                         isld_res = &islands_res[best_island_index][plidx_dst];
1691                                         if (use_from_vert) {
1692                                                 /* Indices stored in islands_res are those of loops, one per dest loop. */
1693                                                 lidx_src = isld_res->index_src;
1694                                                 if (lidx_src >= 0) {
1695                                                         pidx_src = loop_to_poly_map_src[lidx_src];
1696                                                         /* If prev and curr poly are the same, no need to do anything more!!! */
1697                                                         if (!ELEM(pidx_src_prev, -1, pidx_src) && isld_steps_src) {
1698                                                                 int pidx_isld_src, pidx_isld_src_prev;
1699                                                                 if (poly_island_index_map) {
1700                                                                         pidx_isld_src = poly_island_index_map[pidx_src];
1701                                                                         pidx_isld_src_prev = poly_island_index_map[pidx_src_prev];
1702                                                                 }
1703                                                                 else {
1704                                                                         pidx_isld_src = pidx_src;
1705                                                                         pidx_isld_src_prev = pidx_src_prev;
1706                                                                 }
1707
1708                                                                 BLI_astar_graph_solve(
1709                                                                         as_graph, pidx_isld_src_prev, pidx_isld_src,
1710                                                                         mesh_remap_calc_loops_astar_f_cost, &as_solution, isld_steps_src);
1711                                                                 if (GET_INT_FROM_POINTER(as_solution.custom_data) && (as_solution.steps > 0)) {
1712                                                                         /* Find first 'cutting edge' on path, and bring back lidx_src on poly just
1713                                                                          * before that edge.
1714                                                                          * Note we could try to be much smarter (like e.g. storing a whole poly's indices,
1715                                                                          * and making decision (on which side of cutting edge(s!) to be) on the end,
1716                                                                          * but this is one more level of complexity, better to first see if
1717                                                                          * simple solution works!
1718                                                                          */
1719                                                                         int last_valid_pidx_isld_src = -1;
1720                                                                         /* Note we go backward here, from dest to src poly. */
1721                                                                         for (i = as_solution.steps - 1; i--;) {
1722                                                                                 BLI_AStarGNLink *as_link = as_solution.prev_links[pidx_isld_src];
1723                                                                                 const int eidx = GET_INT_FROM_POINTER(as_link->custom_data);
1724                                                                                 pidx_isld_src = as_solution.prev_nodes[pidx_isld_src];
1725                                                                                 BLI_assert(pidx_isld_src != -1);
1726                                                                                 if (eidx != -1) {
1727                                                                                         /* we are 'crossing' a cutting edge. */
1728                                                                                         last_valid_pidx_isld_src = pidx_isld_src;
1729                                                                                 }
1730                                                                         }
1731                                                                         if (last_valid_pidx_isld_src != -1) {
1732                                                                                 /* Find a new valid loop in that new poly (nearest one for now).
1733                                                                                  * Note we could be much more subtle here, again that's for later... */
1734                                                                                 int j;
1735                                                                                 float best_dist_sq = FLT_MAX;
1736
1737                                                                                 ml_dst = &loops_dst[lidx_dst];
1738                                                                                 copy_v3_v3(tmp_co, verts_dst[ml_dst->v].co);
1739
1740                                                                                 /* We do our transform here, since we may do several raycast/nearest queries. */
1741                                                                                 if (space_transform) {
1742                                                                                         BLI_space_transform_apply(space_transform, tmp_co);
1743                                                                                 }
1744
1745                                                                                 pidx_src = use_islands ? best_island->indices[last_valid_pidx_isld_src] :
1746                                                                                                          last_valid_pidx_isld_src;
1747                                                                                 mp_src = &polys_src[pidx_src];
1748                                                                                 ml_src = &loops_src[mp_src->loopstart];
1749                                                                                 for (j = 0; j < mp_src->totloop; j++, ml_src++) {
1750                                                                                         const float dist_sq = len_squared_v3v3(verts_src[ml_src->v].co, tmp_co);
1751                                                                                         if (dist_sq < best_dist_sq) {
1752                                                                                                 best_dist_sq = dist_sq;
1753                                                                                                 lidx_src = mp_src->loopstart + j;
1754                                                                                         }
1755                                                                                 }
1756                                                                         }
1757                                                                 }
1758                                                         }
1759                                                         mesh_remap_item_define(
1760                                                                 r_map, lidx_dst, isld_res->hit_dist,
1761                                                                 best_island_index, 1, &lidx_src, &full_weight);
1762                                                         pidx_src_prev = pidx_src;
1763                                                 }
1764                                                 else {
1765                                                         /* No source for this loop in this island. */
1766                                                         /* TODO: would probably be better to get a source at all cost in best island anyway? */
1767                                                         mesh_remap_item_define(
1768                                                                 r_map, lidx_dst, FLT_MAX,
1769                                                                 best_island_index, 0, NULL, NULL);
1770                                                 }
1771                                         }
1772                                         else {
1773                                                 /* Else, we use source poly, indices stored in islands_res are those of polygons. */
1774                                                 pidx_src = isld_res->index_src;
1775                                                 if (pidx_src >= 0) {
1776                                                         float *hit_co = isld_res->hit_point;
1777                                                         int best_loop_index_src;
1778
1779                                                         mp_src = &polys_src[pidx_src];
1780                                                         /* If prev and curr poly are the same, no need to do anything more!!! */
1781                                                         if (!ELEM(pidx_src_prev, -1, pidx_src) && isld_steps_src) {
1782                                                                 int pidx_isld_src, pidx_isld_src_prev;
1783                                                                 if (poly_island_index_map) {
1784                                                                         pidx_isld_src = poly_island_index_map[pidx_src];
1785                                                                         pidx_isld_src_prev = poly_island_index_map[pidx_src_prev];
1786                                                                 }
1787                                                                 else {
1788                                                                         pidx_isld_src = pidx_src;
1789                                                                         pidx_isld_src_prev = pidx_src_prev;
1790                                                                 }
1791
1792                                                                 BLI_astar_graph_solve(
1793                                                                         as_graph, pidx_isld_src_prev, pidx_isld_src,
1794                                                                         mesh_remap_calc_loops_astar_f_cost, &as_solution, isld_steps_src);
1795                                                                 if (GET_INT_FROM_POINTER(as_solution.custom_data) && (as_solution.steps > 0)) {
1796                                                                         /* Find first 'cutting edge' on path, and bring back lidx_src on poly just
1797                                                                          * before that edge.
1798                                                                          * Note we could try to be much smarter (like e.g. storing a whole poly's indices,
1799                                                                          * and making decision (one which side of cutting edge(s!) to be on the end,
1800                                                                          * but this is one more level of complexity, better to first see if
1801                                                                          * simple solution works!
1802                                                                          */
1803                                                                         int last_valid_pidx_isld_src = -1;
1804                                                                         /* Note we go backward here, from dest to src poly. */
1805                                                                         for (i = as_solution.steps - 1; i--;) {
1806                                                                                 BLI_AStarGNLink *as_link = as_solution.prev_links[pidx_isld_src];
1807                                                                                 int eidx = GET_INT_FROM_POINTER(as_link->custom_data);
1808
1809                                                                                 pidx_isld_src = as_solution.prev_nodes[pidx_isld_src];
1810                                                                                 BLI_assert(pidx_isld_src != -1);
1811                                                                                 if (eidx != -1) {
1812                                                                                         /* we are 'crossing' a cutting edge. */
1813                                                                                         last_valid_pidx_isld_src = pidx_isld_src;
1814                                                                                 }
1815                                                                         }
1816                                                                         if (last_valid_pidx_isld_src != -1) {
1817                                                                                 /* Find a new valid loop in that new poly (nearest point on poly for now).
1818                                                                                  * Note we could be much more subtle here, again that's for later... */
1819                                                                                 float best_dist_sq = FLT_MAX;
1820                                                                                 int j;
1821
1822                                                                                 ml_dst = &loops_dst[lidx_dst];
1823                                                                                 copy_v3_v3(tmp_co, verts_dst[ml_dst->v].co);
1824
1825                                                                                 /* We do our transform here, since we may do several raycast/nearest queries. */
1826                                                                                 if (space_transform) {
1827                                                                                         BLI_space_transform_apply(space_transform, tmp_co);
1828                                                                                 }
1829
1830                                                                                 pidx_src = use_islands ? best_island->indices[last_valid_pidx_isld_src] :
1831                                                                                                          last_valid_pidx_isld_src;
1832                                                                                 mp_src = &polys_src[pidx_src];
1833
1834                                                                                 /* Create that one on demand. */
1835                                                                                 if (poly_to_looptri_map_src == NULL) {
1836                                                                                         BKE_mesh_origindex_map_create_looptri(
1837                                                                                                 &poly_to_looptri_map_src, &poly_to_looptri_map_src_buff,
1838                                                                                                 polys_src, num_polys_src,
1839                                                                                                 looptri_src, num_looptri_src);
1840                                                                                 }
1841
1842                                                                                 for (j = poly_to_looptri_map_src[pidx_src].count; j--;) {
1843                                                                                         float h[3];
1844                                                                                         const MLoopTri *lt = &looptri_src[poly_to_looptri_map_src[pidx_src].indices[j]];
1845                                                                                         float dist_sq;
1846
1847                                                                                         closest_on_tri_to_point_v3(
1848                                                                                                 h, tmp_co,
1849                                                                                                 vcos_src[loops_src[lt->tri[0]].v],
1850                                                                                                 vcos_src[loops_src[lt->tri[1]].v],
1851                                                                                                 vcos_src[loops_src[lt->tri[2]].v]);
1852                                                                                         dist_sq = len_squared_v3v3(tmp_co, h);
1853                                                                                         if (dist_sq < best_dist_sq) {
1854                                                                                                 copy_v3_v3(hit_co, h);
1855                                                                                                 best_dist_sq = dist_sq;
1856                                                                                         }
1857                                                                                 }
1858                                                                         }
1859                                                                 }
1860                                                         }
1861
1862                                                         if (mode == MREMAP_MODE_LOOP_POLY_NEAREST) {
1863                                                                 mesh_remap_interp_poly_data_get(
1864                                                                         mp_src, loops_src, (const float (*)[3])vcos_src, hit_co,
1865                                                                         &buff_size_interp, &vcos_interp, true, &indices_interp,
1866                                                                         &weights_interp, false, &best_loop_index_src);
1867
1868                                                                 mesh_remap_item_define(
1869                                                                         r_map, lidx_dst, isld_res->hit_dist,
1870                                                                         best_island_index, 1, &best_loop_index_src, &full_weight);
1871                                                         }
1872                                                         else {
1873                                                                 const int sources_num = mesh_remap_interp_poly_data_get(
1874                                                                         mp_src, loops_src, (const float (*)[3])vcos_src, hit_co,
1875                                                                         &buff_size_interp, &vcos_interp, true, &indices_interp,
1876                                                                         &weights_interp, true, NULL);
1877
1878                                                                 mesh_remap_item_define(
1879                                                                         r_map, lidx_dst,
1880                                                                         isld_res->hit_dist, best_island_index,
1881                                                                         sources_num, indices_interp, weights_interp);
1882                                                         }
1883
1884                                                         pidx_src_prev = pidx_src;
1885                                                 }
1886                                                 else {
1887                                                         /* No source for this loop in this island. */
1888                                                         /* TODO: would probably be better to get a source at all cost in best island anyway? */
1889                                                         mesh_remap_item_define(r_map, lidx_dst, FLT_MAX, best_island_index, 0, NULL, NULL);
1890                                                 }
1891                                         }
1892                                 }
1893
1894                                 BLI_astar_solution_clear(&as_solution);
1895                         }
1896                 }
1897
1898                 for (tindex = 0; tindex < num_trees; tindex++) {
1899                         MEM_freeN(islands_res[tindex]);
1900                         free_bvhtree_from_mesh(&treedata[tindex]);
1901                         if (isld_steps_src) {
1902                                 BLI_astar_graph_free(&as_graphdata[tindex]);
1903                         }
1904                 }
1905                 MEM_freeN(islands_res);
1906                 BKE_mesh_loop_islands_free(&island_store);
1907                 MEM_freeN(treedata);
1908                 if (isld_steps_src) {
1909                         MEM_freeN(as_graphdata);
1910                         BLI_astar_solution_free(&as_solution);
1911                 }
1912
1913                 if (verts_allocated_src) {
1914                         MEM_freeN(verts_src);
1915                 }
1916                 if (vcos_src) {
1917                         MEM_freeN(vcos_src);
1918                 }
1919                 if (edges_allocated_src) {
1920                         MEM_freeN(edges_src);
1921                 }
1922                 if (loops_allocated_src) {
1923                         MEM_freeN(loops_src);
1924                 }
1925                 if (polys_allocated_src) {
1926                         MEM_freeN(polys_src);
1927                 }
1928                 if (vert_to_loop_map_src) {
1929                         MEM_freeN(vert_to_loop_map_src);
1930                 }
1931                 if (vert_to_loop_map_src_buff) {
1932                         MEM_freeN(vert_to_loop_map_src_buff);
1933                 }
1934                 if (vert_to_poly_map_src) {
1935                         MEM_freeN(vert_to_poly_map_src);
1936                 }
1937                 if (vert_to_poly_map_src_buff) {
1938                         MEM_freeN(vert_to_poly_map_src_buff);
1939                 }
1940                 if (edge_to_poly_map_src) {
1941                         MEM_freeN(edge_to_poly_map_src);
1942                 }
1943                 if (edge_to_poly_map_src_buff) {
1944                         MEM_freeN(edge_to_poly_map_src_buff);
1945                 }
1946                 if (poly_to_looptri_map_src) {
1947                         MEM_freeN(poly_to_looptri_map_src);
1948                 }
1949                 if (poly_to_looptri_map_src_buff) {
1950                         MEM_freeN(poly_to_looptri_map_src_buff);
1951                 }
1952                 if (loop_to_poly_map_src) {
1953                         MEM_freeN(loop_to_poly_map_src);
1954                 }
1955                 if (poly_cents_src) {
1956                         MEM_freeN(poly_cents_src);
1957                 }
1958                 if (vcos_interp) {
1959                         MEM_freeN(vcos_interp);
1960                 }
1961                 if (indices_interp) {
1962                         MEM_freeN(indices_interp);
1963                 }
1964                 if (weights_interp) {
1965                         MEM_freeN(weights_interp);
1966                 }
1967         }
1968 }
1969
1970 void BKE_mesh_remap_calc_polys_from_dm(
1971         const int mode, const SpaceTransform *space_transform, const float max_dist, const float ray_radius,
1972         MVert *verts_dst, const int numverts_dst, MLoop *loops_dst, const int numloops_dst,
1973         MPoly *polys_dst, const int numpolys_dst, CustomData *pdata_dst, const bool dirty_nors_dst,
1974         DerivedMesh *dm_src, MeshPairRemap *r_map)
1975 {
1976         const float full_weight = 1.0f;
1977         const float max_dist_sq = max_dist * max_dist;
1978         float (*poly_nors_dst)[3] = NULL;
1979         float tmp_co[3], tmp_no[3];
1980         int i;
1981
1982         BLI_assert(mode & MREMAP_MODE_POLY);
1983
1984         if (mode & (MREMAP_USE_NORMAL | MREMAP_USE_NORPROJ)) {
1985                 /* Cache poly nors into a temp CDLayer. */
1986                 poly_nors_dst = CustomData_get_layer(pdata_dst, CD_NORMAL);
1987                 if (!poly_nors_dst) {
1988                         poly_nors_dst = CustomData_add_layer(pdata_dst, CD_NORMAL, CD_CALLOC, NULL, numpolys_dst);
1989                         CustomData_set_layer_flag(pdata_dst, CD_NORMAL, CD_FLAG_TEMPORARY);
1990                 }
1991                 if (dirty_nors_dst) {
1992                         BKE_mesh_calc_normals_poly(verts_dst, NULL, numverts_dst, loops_dst, polys_dst, numloops_dst, numpolys_dst,
1993                                                    poly_nors_dst, true);
1994                 }
1995         }
1996
1997         BKE_mesh_remap_init(r_map, numpolys_dst);
1998
1999         if (mode == MREMAP_MODE_TOPOLOGY) {
2000                 BLI_assert(numpolys_dst == dm_src->getNumPolys(dm_src));
2001                 for (i = 0; i < numpolys_dst; i++) {
2002                         mesh_remap_item_define(r_map, i, FLT_MAX, 0, 1, &i, &full_weight);
2003                 }
2004         }
2005         else {
2006                 BVHTreeFromMesh treedata = {NULL};
2007                 BVHTreeNearest nearest = {0};
2008                 BVHTreeRayHit rayhit = {0};
2009                 float hit_dist;
2010
2011                 bvhtree_from_mesh_get(&treedata, dm_src, BVHTREE_FROM_LOOPTRI, 2);
2012
2013                 if (mode == MREMAP_MODE_POLY_NEAREST) {
2014                         nearest.index = -1;
2015
2016                         for (i = 0; i < numpolys_dst; i++) {
2017                                 MPoly *mp = &polys_dst[i];
2018
2019                                 BKE_mesh_calc_poly_center(mp, &loops_dst[mp->loopstart], verts_dst, tmp_co);
2020
2021                                 /* Convert the vertex to tree coordinates, if needed. */
2022                                 if (space_transform) {
2023                                         BLI_space_transform_apply(space_transform, tmp_co);
2024                                 }
2025
2026                                 if (mesh_remap_bvhtree_query_nearest(&treedata, &nearest, tmp_co, max_dist_sq, &hit_dist)) {
2027                                         const MLoopTri *lt = &treedata.looptri[nearest.index];
2028                                         const int poly_index = (int)lt->poly;
2029                                         mesh_remap_item_define(
2030                                                 r_map, i, hit_dist, 0,
2031                                                 1, &poly_index, &full_weight);
2032                                 }
2033                                 else {
2034                                         /* No source for this dest poly! */
2035                                         BKE_mesh_remap_item_define_invalid(r_map, i);
2036                                 }
2037                         }
2038                 }
2039                 else if (mode == MREMAP_MODE_POLY_NOR) {
2040                         BLI_assert(poly_nors_dst);
2041
2042                         for (i = 0; i < numpolys_dst; i++) {
2043                                 MPoly *mp = &polys_dst[i];
2044
2045                                 BKE_mesh_calc_poly_center(mp, &loops_dst[mp->loopstart], verts_dst, tmp_co);
2046                                 copy_v3_v3(tmp_no, poly_nors_dst[i]);
2047
2048                                 /* Convert the vertex to tree coordinates, if needed. */
2049                                 if (space_transform) {
2050                                         BLI_space_transform_apply(space_transform, tmp_co);
2051                                         BLI_space_transform_apply_normal(space_transform, tmp_no);
2052                                 }
2053
2054                                 treedata.sphere_radius = ray_radius;
2055                                 if (mesh_remap_bvhtree_query_raycast(
2056                                         &treedata, &rayhit, tmp_co, tmp_no, ray_radius, max_dist, &hit_dist))
2057                                 {
2058                                         const MLoopTri *lt = &treedata.looptri[rayhit.index];
2059                                         const int poly_index = (int)lt->poly;
2060
2061                                         mesh_remap_item_define(
2062                                                 r_map, i, hit_dist, 0,
2063                                                 1, &poly_index, &full_weight);
2064                                 }
2065                                 else {
2066                                         /* No source for this dest poly! */
2067                                         BKE_mesh_remap_item_define_invalid(r_map, i);
2068                                 }
2069                         }
2070                 }
2071                 else if (mode == MREMAP_MODE_POLY_POLYINTERP_PNORPROJ) {
2072                         /* We cast our rays randomly, with a pseudo-even distribution (since we spread across tessellated tris,
2073                          * with additional weighting based on each tri's relative area).
2074                          */
2075                         RNG *rng = BLI_rng_new(0);
2076
2077                         const size_t numpolys_src = (size_t)dm_src->getNumPolys(dm_src);
2078
2079                         /* Here it's simpler to just allocate for all polys :/ */
2080                         int *indices = MEM_mallocN(sizeof(*indices) * numpolys_src, __func__);
2081                         float *weights = MEM_mallocN(sizeof(*weights) * numpolys_src, __func__);
2082
2083                         size_t tmp_poly_size = MREMAP_DEFAULT_BUFSIZE;
2084                         float (*poly_vcos_2d)[2] = MEM_mallocN(sizeof(*poly_vcos_2d) * tmp_poly_size, __func__);
2085                         /* Tessellated 2D poly, always (num_loops - 2) triangles. */
2086                         int (*tri_vidx_2d)[3] = MEM_mallocN(sizeof(*tri_vidx_2d) * (tmp_poly_size - 2), __func__);
2087
2088                         for (i = 0; i < numpolys_dst; i++) {
2089                                 /* For each dst poly, we sample some rays from it (2D grid in pnor space)
2090                                  * and use their hits to interpolate from source polys. */
2091                                 /* Note: dst poly is early-converted into src space! */
2092                                 MPoly *mp = &polys_dst[i];
2093
2094                                 int tot_rays, done_rays = 0;
2095                                 float poly_area_2d_inv, done_area = 0.0f;
2096
2097                                 float pcent_dst[3];
2098                                 float to_pnor_2d_mat[3][3], from_pnor_2d_mat[3][3];
2099                                 float poly_dst_2d_min[2], poly_dst_2d_max[2], poly_dst_2d_z;
2100                                 float poly_dst_2d_size[2];
2101
2102                                 float totweights = 0.0f;
2103                                 float hit_dist_accum = 0.0f;
2104                                 int sources_num = 0;
2105                                 const int tris_num = mp->totloop - 2;
2106                                 int j;
2107
2108                                 BKE_mesh_calc_poly_center(mp, &loops_dst[mp->loopstart], verts_dst, pcent_dst);
2109                                 copy_v3_v3(tmp_no, poly_nors_dst[i]);
2110
2111                                 /* We do our transform here, else it'd be redone by raycast helper for each ray, ugh! */
2112                                 if (space_transform) {
2113                                         BLI_space_transform_apply(space_transform, pcent_dst);
2114                                         BLI_space_transform_apply_normal(space_transform, tmp_no);
2115                                 }
2116
2117                                 copy_vn_fl(weights, (int)numpolys_src, 0.0f);
2118
2119                                 if (UNLIKELY((size_t)mp->totloop > tmp_poly_size)) {
2120                                         tmp_poly_size = (size_t)mp->totloop;
2121                                         poly_vcos_2d = MEM_reallocN(poly_vcos_2d, sizeof(*poly_vcos_2d) * tmp_poly_size);
2122                                         tri_vidx_2d = MEM_reallocN(tri_vidx_2d, sizeof(*tri_vidx_2d) * (tmp_poly_size - 2));
2123                                 }
2124
2125                                 axis_dominant_v3_to_m3(to_pnor_2d_mat, tmp_no);
2126                                 invert_m3_m3(from_pnor_2d_mat, to_pnor_2d_mat);
2127
2128                                 mul_m3_v3(to_pnor_2d_mat, pcent_dst);
2129                                 poly_dst_2d_z = pcent_dst[2];
2130
2131                                 /* Get (2D) bounding square of our poly. */
2132                                 INIT_MINMAX2(poly_dst_2d_min, poly_dst_2d_max);
2133
2134                                 for (j = 0; j < mp->totloop; j++) {
2135                                         MLoop *ml = &loops_dst[j + mp->loopstart];
2136                                         copy_v3_v3(tmp_co, verts_dst[ml->v].co);
2137                                         if (space_transform) {
2138                                                 BLI_space_transform_apply(space_transform, tmp_co);
2139                                         }
2140                                         mul_v2_m3v3(poly_vcos_2d[j], to_pnor_2d_mat, tmp_co);
2141                                         minmax_v2v2_v2(poly_dst_2d_min, poly_dst_2d_max, poly_vcos_2d[j]);
2142                                 }
2143
2144                                 /* We adjust our ray-casting grid to ray_radius (the smaller, the more rays are cast),
2145                                  * with lower/upper bounds. */
2146                                 sub_v2_v2v2(poly_dst_2d_size, poly_dst_2d_max, poly_dst_2d_min);
2147
2148                                 if (ray_radius) {
2149                                         tot_rays = (int)((max_ff(poly_dst_2d_size[0], poly_dst_2d_size[1]) / ray_radius) + 0.5f);
2150                                         CLAMP(tot_rays, MREMAP_RAYCAST_TRI_SAMPLES_MIN, MREMAP_RAYCAST_TRI_SAMPLES_MAX);
2151                                 }
2152                                 else {
2153                                         /* If no radius (pure rays), give max number of rays! */
2154                                         tot_rays = MREMAP_RAYCAST_TRI_SAMPLES_MIN;
2155                                 }
2156                                 tot_rays *= tot_rays;
2157
2158                                 poly_area_2d_inv = area_poly_v2((const float(*)[2])poly_vcos_2d, (unsigned int)mp->totloop);
2159                                 /* In case we have a null-area degenerated poly... */
2160                                 poly_area_2d_inv = 1.0f / max_ff(poly_area_2d_inv, 1e-9f);
2161
2162                                 /* Tessellate our poly. */
2163                                 if (mp->totloop == 3) {
2164                                         tri_vidx_2d[0][0] = 0;
2165                                         tri_vidx_2d[0][1] = 1;
2166                                         tri_vidx_2d[0][2] = 2;
2167                                 }
2168                                 if (mp->totloop == 4) {
2169                                         tri_vidx_2d[0][0] = 0;
2170                                         tri_vidx_2d[0][1] = 1;
2171                                         tri_vidx_2d[0][2] = 2;
2172                                         tri_vidx_2d[1][0] = 0;
2173                                         tri_vidx_2d[1][1] = 2;
2174                                         tri_vidx_2d[1][2] = 3;
2175                                 }
2176                                 else {
2177                                         BLI_polyfill_calc(poly_vcos_2d, (unsigned int)mp->totloop, -1,
2178                                                           (unsigned int (*)[3])tri_vidx_2d);
2179                                 }
2180
2181                                 for (j = 0; j < tris_num; j++) {
2182                                         float *v1 = poly_vcos_2d[tri_vidx_2d[j][0]];
2183                                         float *v2 = poly_vcos_2d[tri_vidx_2d[j][1]];
2184                                         float *v3 = poly_vcos_2d[tri_vidx_2d[j][2]];
2185                                         int rays_num;
2186
2187                                         /* All this allows us to get 'absolute' number of rays for each tri, avoiding accumulating
2188                                          * errors over iterations, and helping better even distribution. */
2189                                         done_area += area_tri_v2(v1, v2, v3);
2190                                         rays_num = max_ii((int)((float)tot_rays * done_area * poly_area_2d_inv + 0.5f) - done_rays, 0);
2191                                         done_rays += rays_num;
2192
2193                                         while (rays_num--) {
2194                                                 int n = (ray_radius > 0.0f) ? MREMAP_RAYCAST_APPROXIMATE_NR : 1;
2195                                                 float w = 1.0f;
2196
2197                                                 BLI_rng_get_tri_sample_float_v2(rng, v1, v2, v3, tmp_co);
2198
2199                                                 tmp_co[2] = poly_dst_2d_z;
2200                                                 mul_m3_v3(from_pnor_2d_mat, tmp_co);
2201
2202                                                 /* At this point, tmp_co is a point on our poly surface, in mesh_src space! */
2203                                                 while (n--) {
2204                                                         treedata.sphere_radius = ray_radius / w;
2205                                                         if (mesh_remap_bvhtree_query_raycast(
2206                                                                 &treedata, &rayhit, tmp_co, tmp_no, treedata.sphere_radius, max_dist, &hit_dist))
2207                                                         {
2208                                                                 const MLoopTri *lt = &treedata.looptri[rayhit.index];
2209
2210                                                                 weights[lt->poly] += w;
2211                                                                 totweights += w;
2212                                                                 hit_dist_accum += hit_dist;
2213                                                                 break;
2214                                                         }
2215                                                         /* Next iteration will get bigger radius but smaller weight! */
2216                                                         w /= MREMAP_RAYCAST_APPROXIMATE_FAC;
2217                                                 }
2218                                         }
2219                                 }
2220
2221                                 if (totweights > 0.0f) {
2222                                         for (j = 0; j < (int)numpolys_src; j++) {
2223                                                 if (!weights[j]) {
2224                                                         continue;
2225                                                 }
2226                                                 /* Note: sources_num is always <= j! */
2227                                                 weights[sources_num] = weights[j] / totweights;
2228                                                 indices[sources_num] = j;
2229                                                 sources_num++;
2230                                         }
2231                                         mesh_remap_item_define(r_map, i, hit_dist_accum / totweights, 0, sources_num, indices, weights);
2232                                 }
2233                                 else {
2234                                         /* No source for this dest poly! */
2235                                         BKE_mesh_remap_item_define_invalid(r_map, i);
2236                                 }
2237                         }
2238
2239                         MEM_freeN(tri_vidx_2d);
2240                         MEM_freeN(poly_vcos_2d);
2241                         MEM_freeN(indices);
2242                         MEM_freeN(weights);
2243                         BLI_rng_free(rng);
2244                 }
2245                 else {
2246                         printf("WARNING! Unsupported mesh-to-mesh poly mapping mode (%d)!\n", mode);
2247                         memset(r_map->items, 0, sizeof(*r_map->items) * (size_t)numpolys_dst);
2248                 }
2249
2250                 free_bvhtree_from_mesh(&treedata);
2251         }
2252 }
2253
2254 #undef MREMAP_RAYCAST_APPROXIMATE_NR
2255 #undef MREMAP_RAYCAST_APPROXIMATE_FAC
2256 #undef MREMAP_RAYCAST_TRI_SAMPLES_MIN
2257 #undef MREMAP_RAYCAST_TRI_SAMPLES_MAX
2258 #undef MREMAP_DEFAULT_BUFSIZE
2259
2260 /** \} */