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