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