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