Fix T75283: GPencil Stroke created by `Merge points` can't have fill material
[blender.git] / source / blender / blenkernel / intern / shrinkwrap.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  * The Original Code is Copyright (C) Blender Foundation.
17  * All rights reserved.
18  */
19
20 /** \file
21  * \ingroup bke
22  */
23
24 #include <assert.h>
25 #include <float.h>
26 #include <math.h>
27 #include <memory.h>
28 #include <stdio.h>
29 #include <string.h>
30 #include <time.h>
31
32 #include "DNA_mesh_types.h"
33 #include "DNA_meshdata_types.h"
34 #include "DNA_modifier_types.h"
35 #include "DNA_object_types.h"
36
37 #include "BLI_math.h"
38 #include "BLI_math_solvers.h"
39 #include "BLI_task.h"
40 #include "BLI_utildefines.h"
41
42 #include "BKE_DerivedMesh.h"
43 #include "BKE_cdderivedmesh.h"
44 #include "BKE_context.h"
45 #include "BKE_lattice.h"
46 #include "BKE_lib_id.h"
47 #include "BKE_modifier.h"
48 #include "BKE_shrinkwrap.h"
49
50 #include "BKE_deform.h"
51 #include "BKE_editmesh.h"
52 #include "BKE_mesh.h" /* for OMP limits. */
53 #include "BKE_mesh_runtime.h"
54 #include "BKE_subsurf.h"
55
56 #include "DEG_depsgraph_query.h"
57
58 #include "MEM_guardedalloc.h"
59
60 #include "BLI_strict_flags.h"
61
62 /* for timing... */
63 #if 0
64 #  include "PIL_time_utildefines.h"
65 #else
66 #  define TIMEIT_BENCH(expr, id) (expr)
67 #endif
68
69 /* Util macros */
70 #define OUT_OF_MEMORY() ((void)printf("Shrinkwrap: Out of memory\n"))
71
72 typedef struct ShrinkwrapCalcData {
73   ShrinkwrapModifierData *smd;  // shrinkwrap modifier data
74
75   struct Object *ob;  // object we are applying shrinkwrap to
76
77   struct MVert *vert;     // Array of verts being projected (to fetch normals or other data)
78   float (*vertexCos)[3];  // vertexs being shrinkwraped
79   int numVerts;
80
81   struct MDeformVert *dvert;  // Pointer to mdeform array
82   int vgroup;                 // Vertex group num
83   bool invert_vgroup;         /* invert vertex group influence */
84
85   struct Mesh *target;                 // mesh we are shrinking to
86   struct SpaceTransform local2target;  // transform to move between local and target space
87   struct ShrinkwrapTreeData *tree;     // mesh BVH tree data
88
89   struct Object *aux_target;
90
91   float keepDist;  // Distance to keep above target surface (units are in local space)
92 } ShrinkwrapCalcData;
93
94 typedef struct ShrinkwrapCalcCBData {
95   ShrinkwrapCalcData *calc;
96
97   ShrinkwrapTreeData *tree;
98   ShrinkwrapTreeData *aux_tree;
99
100   float *proj_axis;
101   SpaceTransform *local2aux;
102 } ShrinkwrapCalcCBData;
103
104 /* Checks if the modifier needs target normals with these settings. */
105 bool BKE_shrinkwrap_needs_normals(int shrinkType, int shrinkMode)
106 {
107   return (shrinkType == MOD_SHRINKWRAP_TARGET_PROJECT) ||
108          (shrinkType != MOD_SHRINKWRAP_NEAREST_VERTEX &&
109           shrinkMode == MOD_SHRINKWRAP_ABOVE_SURFACE);
110 }
111
112 /* Initializes the mesh data structure from the given mesh and settings. */
113 bool BKE_shrinkwrap_init_tree(
114     ShrinkwrapTreeData *data, Mesh *mesh, int shrinkType, int shrinkMode, bool force_normals)
115 {
116   memset(data, 0, sizeof(*data));
117
118   if (!mesh || mesh->totvert <= 0) {
119     return false;
120   }
121
122   data->mesh = mesh;
123
124   if (shrinkType == MOD_SHRINKWRAP_NEAREST_VERTEX) {
125     data->bvh = BKE_bvhtree_from_mesh_get(&data->treeData, mesh, BVHTREE_FROM_VERTS, 2);
126
127     return data->bvh != NULL;
128   }
129   else {
130     if (mesh->totpoly <= 0) {
131       return false;
132     }
133
134     data->bvh = BKE_bvhtree_from_mesh_get(&data->treeData, mesh, BVHTREE_FROM_LOOPTRI, 4);
135
136     if (data->bvh == NULL) {
137       return false;
138     }
139
140     if (force_normals || BKE_shrinkwrap_needs_normals(shrinkType, shrinkMode)) {
141       data->pnors = CustomData_get_layer(&mesh->pdata, CD_NORMAL);
142       if ((mesh->flag & ME_AUTOSMOOTH) != 0) {
143         data->clnors = CustomData_get_layer(&mesh->ldata, CD_NORMAL);
144       }
145     }
146
147     if (shrinkType == MOD_SHRINKWRAP_TARGET_PROJECT) {
148       data->boundary = mesh->runtime.shrinkwrap_data;
149     }
150
151     return true;
152   }
153 }
154
155 /* Frees the tree data if necessary. */
156 void BKE_shrinkwrap_free_tree(ShrinkwrapTreeData *data)
157 {
158   free_bvhtree_from_mesh(&data->treeData);
159 }
160
161 /* Free boundary data for target project */
162 void BKE_shrinkwrap_discard_boundary_data(struct Mesh *mesh)
163 {
164   struct ShrinkwrapBoundaryData *data = mesh->runtime.shrinkwrap_data;
165
166   if (data != NULL) {
167     MEM_freeN((void *)data->edge_is_boundary);
168     MEM_freeN((void *)data->looptri_has_boundary);
169     MEM_freeN((void *)data->vert_boundary_id);
170     MEM_freeN((void *)data->boundary_verts);
171
172     MEM_freeN(data);
173   }
174
175   mesh->runtime.shrinkwrap_data = NULL;
176 }
177
178 /* Accumulate edge for average boundary edge direction. */
179 static void merge_vert_dir(ShrinkwrapBoundaryVertData *vdata,
180                            signed char *status,
181                            int index,
182                            const float edge_dir[3],
183                            signed char side)
184 {
185   BLI_assert(index >= 0);
186   float *direction = vdata[index].direction;
187
188   /* Invert the direction vector if either:
189    * - This is the second edge and both edges have the vertex as start or end.
190    * - For third and above, if it points in the wrong direction.
191    */
192   if (status[index] >= 0 ? status[index] == side : dot_v3v3(direction, edge_dir) < 0) {
193     sub_v3_v3(direction, edge_dir);
194   }
195   else {
196     add_v3_v3(direction, edge_dir);
197   }
198
199   status[index] = (status[index] == 0) ? side : -1;
200 }
201
202 static ShrinkwrapBoundaryData *shrinkwrap_build_boundary_data(struct Mesh *mesh)
203 {
204   const MLoop *mloop = mesh->mloop;
205   const MEdge *medge = mesh->medge;
206   const MVert *mvert = mesh->mvert;
207
208   /* Count faces per edge (up to 2). */
209   char *edge_mode = MEM_calloc_arrayN((size_t)mesh->totedge, sizeof(char), __func__);
210
211   for (int i = 0; i < mesh->totloop; i++) {
212     unsigned int eidx = mloop[i].e;
213
214     if (edge_mode[eidx] < 2) {
215       edge_mode[eidx]++;
216     }
217   }
218
219   /* Build the boundary edge bitmask. */
220   BLI_bitmap *edge_is_boundary = BLI_BITMAP_NEW(mesh->totedge,
221                                                 "ShrinkwrapBoundaryData::edge_is_boundary");
222   unsigned int num_boundary_edges = 0;
223
224   for (int i = 0; i < mesh->totedge; i++) {
225     edge_mode[i] = (edge_mode[i] == 1);
226
227     if (edge_mode[i]) {
228       BLI_BITMAP_ENABLE(edge_is_boundary, i);
229       num_boundary_edges++;
230     }
231   }
232
233   /* If no boundary, return NULL. */
234   if (num_boundary_edges == 0) {
235     MEM_freeN(edge_is_boundary);
236     MEM_freeN(edge_mode);
237     return NULL;
238   }
239
240   /* Allocate the data object. */
241   ShrinkwrapBoundaryData *data = MEM_callocN(sizeof(ShrinkwrapBoundaryData),
242                                              "ShrinkwrapBoundaryData");
243
244   data->edge_is_boundary = edge_is_boundary;
245
246   /* Build the boundary looptri bitmask. */
247   const MLoopTri *mlooptri = BKE_mesh_runtime_looptri_ensure(mesh);
248   int totlooptri = BKE_mesh_runtime_looptri_len(mesh);
249
250   BLI_bitmap *looptri_has_boundary = BLI_BITMAP_NEW(totlooptri,
251                                                     "ShrinkwrapBoundaryData::looptri_is_boundary");
252
253   for (int i = 0; i < totlooptri; i++) {
254     int edges[3];
255     BKE_mesh_looptri_get_real_edges(mesh, &mlooptri[i], edges);
256
257     for (int j = 0; j < 3; j++) {
258       if (edges[j] >= 0 && edge_mode[edges[j]]) {
259         BLI_BITMAP_ENABLE(looptri_has_boundary, i);
260         break;
261       }
262     }
263   }
264
265   data->looptri_has_boundary = looptri_has_boundary;
266
267   /* Find boundary vertices and build a mapping table for compact storage of data. */
268   int *vert_boundary_id = MEM_calloc_arrayN(
269       (size_t)mesh->totvert, sizeof(int), "ShrinkwrapBoundaryData::vert_boundary_id");
270
271   for (int i = 0; i < mesh->totedge; i++) {
272     if (edge_mode[i]) {
273       const MEdge *edge = &medge[i];
274
275       vert_boundary_id[edge->v1] = 1;
276       vert_boundary_id[edge->v2] = 1;
277     }
278   }
279
280   unsigned int num_boundary_verts = 0;
281
282   for (int i = 0; i < mesh->totvert; i++) {
283     vert_boundary_id[i] = (vert_boundary_id[i] != 0) ? (int)num_boundary_verts++ : -1;
284   }
285
286   data->vert_boundary_id = vert_boundary_id;
287   data->num_boundary_verts = num_boundary_verts;
288
289   /* Compute average directions. */
290   ShrinkwrapBoundaryVertData *boundary_verts = MEM_calloc_arrayN(
291       num_boundary_verts, sizeof(*boundary_verts), "ShrinkwrapBoundaryData::boundary_verts");
292
293   signed char *vert_status = MEM_calloc_arrayN(num_boundary_verts, sizeof(char), __func__);
294
295   for (int i = 0; i < mesh->totedge; i++) {
296     if (edge_mode[i]) {
297       const MEdge *edge = &medge[i];
298
299       float dir[3];
300       sub_v3_v3v3(dir, mvert[edge->v2].co, mvert[edge->v1].co);
301       normalize_v3(dir);
302
303       merge_vert_dir(boundary_verts, vert_status, vert_boundary_id[edge->v1], dir, 1);
304       merge_vert_dir(boundary_verts, vert_status, vert_boundary_id[edge->v2], dir, 2);
305     }
306   }
307
308   MEM_freeN(vert_status);
309
310   /* Finalize average direction and compute normal. */
311   for (int i = 0; i < mesh->totvert; i++) {
312     int bidx = vert_boundary_id[i];
313
314     if (bidx >= 0) {
315       ShrinkwrapBoundaryVertData *vdata = &boundary_verts[bidx];
316       float no[3], tmp[3];
317
318       normalize_v3(vdata->direction);
319
320       normal_short_to_float_v3(no, mesh->mvert[i].no);
321       cross_v3_v3v3(tmp, no, vdata->direction);
322       cross_v3_v3v3(vdata->normal_plane, tmp, no);
323       normalize_v3(vdata->normal_plane);
324     }
325   }
326
327   data->boundary_verts = boundary_verts;
328
329   MEM_freeN(edge_mode);
330   return data;
331 }
332
333 void BKE_shrinkwrap_compute_boundary_data(struct Mesh *mesh)
334 {
335   BKE_shrinkwrap_discard_boundary_data(mesh);
336
337   mesh->runtime.shrinkwrap_data = shrinkwrap_build_boundary_data(mesh);
338 }
339
340 /*
341  * Shrinkwrap to the nearest vertex
342  *
343  * it builds a kdtree of vertexs we can attach to and then
344  * for each vertex performs a nearest vertex search on the tree
345  */
346 static void shrinkwrap_calc_nearest_vertex_cb_ex(void *__restrict userdata,
347                                                  const int i,
348                                                  const TaskParallelTLS *__restrict tls)
349 {
350   ShrinkwrapCalcCBData *data = userdata;
351
352   ShrinkwrapCalcData *calc = data->calc;
353   BVHTreeFromMesh *treeData = &data->tree->treeData;
354   BVHTreeNearest *nearest = tls->userdata_chunk;
355
356   float *co = calc->vertexCos[i];
357   float tmp_co[3];
358   float weight = BKE_defvert_array_find_weight_safe(calc->dvert, i, calc->vgroup);
359
360   if (calc->invert_vgroup) {
361     weight = 1.0f - weight;
362   }
363
364   if (weight == 0.0f) {
365     return;
366   }
367
368   /* Convert the vertex to tree coordinates */
369   if (calc->vert) {
370     copy_v3_v3(tmp_co, calc->vert[i].co);
371   }
372   else {
373     copy_v3_v3(tmp_co, co);
374   }
375   BLI_space_transform_apply(&calc->local2target, tmp_co);
376
377   /* Use local proximity heuristics (to reduce the nearest search)
378    *
379    * If we already had an hit before.. we assume this vertex is going to have a close hit to that
380    * other vertex so we can initiate the "nearest.dist" with the expected value to that last hit.
381    * This will lead in pruning of the search tree. */
382   if (nearest->index != -1) {
383     nearest->dist_sq = len_squared_v3v3(tmp_co, nearest->co);
384   }
385   else {
386     nearest->dist_sq = FLT_MAX;
387   }
388
389   BLI_bvhtree_find_nearest(treeData->tree, tmp_co, nearest, treeData->nearest_callback, treeData);
390
391   /* Found the nearest vertex */
392   if (nearest->index != -1) {
393     /* Adjusting the vertex weight,
394      * so that after interpolating it keeps a certain distance from the nearest position */
395     if (nearest->dist_sq > FLT_EPSILON) {
396       const float dist = sqrtf(nearest->dist_sq);
397       weight *= (dist - calc->keepDist) / dist;
398     }
399
400     /* Convert the coordinates back to mesh coordinates */
401     copy_v3_v3(tmp_co, nearest->co);
402     BLI_space_transform_invert(&calc->local2target, tmp_co);
403
404     interp_v3_v3v3(co, co, tmp_co, weight); /* linear interpolation */
405   }
406 }
407
408 static void shrinkwrap_calc_nearest_vertex(ShrinkwrapCalcData *calc)
409 {
410   BVHTreeNearest nearest = NULL_BVHTreeNearest;
411
412   /* Setup nearest */
413   nearest.index = -1;
414   nearest.dist_sq = FLT_MAX;
415
416   ShrinkwrapCalcCBData data = {
417       .calc = calc,
418       .tree = calc->tree,
419   };
420   TaskParallelSettings settings;
421   BLI_parallel_range_settings_defaults(&settings);
422   settings.use_threading = (calc->numVerts > BKE_MESH_OMP_LIMIT);
423   settings.userdata_chunk = &nearest;
424   settings.userdata_chunk_size = sizeof(nearest);
425   BLI_task_parallel_range(
426       0, calc->numVerts, &data, shrinkwrap_calc_nearest_vertex_cb_ex, &settings);
427 }
428
429 /*
430  * This function raycast a single vertex and updates the hit if the "hit" is considered valid.
431  * Returns true if "hit" was updated.
432  * Opts control whether an hit is valid or not
433  * Supported options are:
434  * - MOD_SHRINKWRAP_CULL_TARGET_FRONTFACE (front faces hits are ignored)
435  * - MOD_SHRINKWRAP_CULL_TARGET_BACKFACE (back faces hits are ignored)
436  */
437 bool BKE_shrinkwrap_project_normal(char options,
438                                    const float vert[3],
439                                    const float dir[3],
440                                    const float ray_radius,
441                                    const SpaceTransform *transf,
442                                    ShrinkwrapTreeData *tree,
443                                    BVHTreeRayHit *hit)
444 {
445   /* don't use this because this dist value could be incompatible
446    * this value used by the callback for comparing prev/new dist values.
447    * also, at the moment there is no need to have a corrected 'dist' value */
448   // #define USE_DIST_CORRECT
449
450   float tmp_co[3], tmp_no[3];
451   const float *co, *no;
452   BVHTreeRayHit hit_tmp;
453
454   /* Copy from hit (we need to convert hit rays from one space coordinates to the other */
455   memcpy(&hit_tmp, hit, sizeof(hit_tmp));
456
457   /* Apply space transform (TODO readjust dist) */
458   if (transf) {
459     copy_v3_v3(tmp_co, vert);
460     BLI_space_transform_apply(transf, tmp_co);
461     co = tmp_co;
462
463     copy_v3_v3(tmp_no, dir);
464     BLI_space_transform_apply_normal(transf, tmp_no);
465     no = tmp_no;
466
467 #ifdef USE_DIST_CORRECT
468     hit_tmp.dist *= mat4_to_scale(((SpaceTransform *)transf)->local2target);
469 #endif
470   }
471   else {
472     co = vert;
473     no = dir;
474   }
475
476   hit_tmp.index = -1;
477
478   BLI_bvhtree_ray_cast(
479       tree->bvh, co, no, ray_radius, &hit_tmp, tree->treeData.raycast_callback, &tree->treeData);
480
481   if (hit_tmp.index != -1) {
482     /* invert the normal first so face culling works on rotated objects */
483     if (transf) {
484       BLI_space_transform_invert_normal(transf, hit_tmp.no);
485     }
486
487     if (options & MOD_SHRINKWRAP_CULL_TARGET_MASK) {
488       /* apply backface */
489       const float dot = dot_v3v3(dir, hit_tmp.no);
490       if (((options & MOD_SHRINKWRAP_CULL_TARGET_FRONTFACE) && dot <= 0.0f) ||
491           ((options & MOD_SHRINKWRAP_CULL_TARGET_BACKFACE) && dot >= 0.0f)) {
492         return false; /* Ignore hit */
493       }
494     }
495
496     if (transf) {
497       /* Inverting space transform (TODO make coeherent with the initial dist readjust) */
498       BLI_space_transform_invert(transf, hit_tmp.co);
499 #ifdef USE_DIST_CORRECT
500       hit_tmp.dist = len_v3v3(vert, hit_tmp.co);
501 #endif
502     }
503
504     BLI_assert(hit_tmp.dist <= hit->dist);
505
506     memcpy(hit, &hit_tmp, sizeof(hit_tmp));
507     return true;
508   }
509   return false;
510 }
511
512 static void shrinkwrap_calc_normal_projection_cb_ex(void *__restrict userdata,
513                                                     const int i,
514                                                     const TaskParallelTLS *__restrict tls)
515 {
516   ShrinkwrapCalcCBData *data = userdata;
517
518   ShrinkwrapCalcData *calc = data->calc;
519   ShrinkwrapTreeData *tree = data->tree;
520   ShrinkwrapTreeData *aux_tree = data->aux_tree;
521
522   float *proj_axis = data->proj_axis;
523   SpaceTransform *local2aux = data->local2aux;
524
525   BVHTreeRayHit *hit = tls->userdata_chunk;
526
527   const float proj_limit_squared = calc->smd->projLimit * calc->smd->projLimit;
528   float *co = calc->vertexCos[i];
529   float tmp_co[3], tmp_no[3];
530   float weight = BKE_defvert_array_find_weight_safe(calc->dvert, i, calc->vgroup);
531
532   if (calc->invert_vgroup) {
533     weight = 1.0f - weight;
534   }
535
536   if (weight == 0.0f) {
537     return;
538   }
539
540   if (calc->vert != NULL && calc->smd->projAxis == MOD_SHRINKWRAP_PROJECT_OVER_NORMAL) {
541     /* calc->vert contains verts from evaluated mesh.  */
542     /* These coordinates are deformed by vertexCos only for normal projection
543      * (to get correct normals) for other cases calc->verts contains undeformed coordinates and
544      * vertexCos should be used */
545     copy_v3_v3(tmp_co, calc->vert[i].co);
546     normal_short_to_float_v3(tmp_no, calc->vert[i].no);
547   }
548   else {
549     copy_v3_v3(tmp_co, co);
550     copy_v3_v3(tmp_no, proj_axis);
551   }
552
553   hit->index = -1;
554
555   /* TODO: we should use FLT_MAX here, but sweepsphere code isn't prepared for that */
556   hit->dist = BVH_RAYCAST_DIST_MAX;
557
558   bool is_aux = false;
559
560   /* Project over positive direction of axis */
561   if (calc->smd->shrinkOpts & MOD_SHRINKWRAP_PROJECT_ALLOW_POS_DIR) {
562     if (aux_tree) {
563       if (BKE_shrinkwrap_project_normal(0, tmp_co, tmp_no, 0.0, local2aux, aux_tree, hit)) {
564         is_aux = true;
565       }
566     }
567
568     if (BKE_shrinkwrap_project_normal(
569             calc->smd->shrinkOpts, tmp_co, tmp_no, 0.0, &calc->local2target, tree, hit)) {
570       is_aux = false;
571     }
572   }
573
574   /* Project over negative direction of axis */
575   if (calc->smd->shrinkOpts & MOD_SHRINKWRAP_PROJECT_ALLOW_NEG_DIR) {
576     float inv_no[3];
577     negate_v3_v3(inv_no, tmp_no);
578
579     char options = calc->smd->shrinkOpts;
580
581     if ((options & MOD_SHRINKWRAP_INVERT_CULL_TARGET) &&
582         (options & MOD_SHRINKWRAP_CULL_TARGET_MASK)) {
583       options ^= MOD_SHRINKWRAP_CULL_TARGET_MASK;
584     }
585
586     if (aux_tree) {
587       if (BKE_shrinkwrap_project_normal(0, tmp_co, inv_no, 0.0, local2aux, aux_tree, hit)) {
588         is_aux = true;
589       }
590     }
591
592     if (BKE_shrinkwrap_project_normal(
593             options, tmp_co, inv_no, 0.0, &calc->local2target, tree, hit)) {
594       is_aux = false;
595     }
596   }
597
598   /* don't set the initial dist (which is more efficient),
599    * because its calculated in the targets space, we want the dist in our own space */
600   if (proj_limit_squared != 0.0f) {
601     if (hit->index != -1 && len_squared_v3v3(hit->co, co) > proj_limit_squared) {
602       hit->index = -1;
603     }
604   }
605
606   if (hit->index != -1) {
607     if (is_aux) {
608       BKE_shrinkwrap_snap_point_to_surface(aux_tree,
609                                            local2aux,
610                                            calc->smd->shrinkMode,
611                                            hit->index,
612                                            hit->co,
613                                            hit->no,
614                                            calc->keepDist,
615                                            tmp_co,
616                                            hit->co);
617     }
618     else {
619       BKE_shrinkwrap_snap_point_to_surface(tree,
620                                            &calc->local2target,
621                                            calc->smd->shrinkMode,
622                                            hit->index,
623                                            hit->co,
624                                            hit->no,
625                                            calc->keepDist,
626                                            tmp_co,
627                                            hit->co);
628     }
629
630     interp_v3_v3v3(co, co, hit->co, weight);
631   }
632 }
633
634 static void shrinkwrap_calc_normal_projection(ShrinkwrapCalcData *calc)
635 {
636   /* Options about projection direction */
637   float proj_axis[3] = {0.0f, 0.0f, 0.0f};
638
639   /* Raycast and tree stuff */
640
641   /** \note 'hit.dist' is kept in the targets space, this is only used
642    * for finding the best hit, to get the real dist,
643    * measure the len_v3v3() from the input coord to hit.co */
644   BVHTreeRayHit hit;
645
646   /* auxiliary target */
647   Mesh *auxMesh = NULL;
648   ShrinkwrapTreeData *aux_tree = NULL;
649   ShrinkwrapTreeData aux_tree_stack;
650   SpaceTransform local2aux;
651
652   /* If the user doesn't allows to project in any direction of projection axis
653    * then there's nothing todo. */
654   if ((calc->smd->shrinkOpts &
655        (MOD_SHRINKWRAP_PROJECT_ALLOW_POS_DIR | MOD_SHRINKWRAP_PROJECT_ALLOW_NEG_DIR)) == 0) {
656     return;
657   }
658
659   /* Prepare data to retrieve the direction in which we should project each vertex */
660   if (calc->smd->projAxis == MOD_SHRINKWRAP_PROJECT_OVER_NORMAL) {
661     if (calc->vert == NULL) {
662       return;
663     }
664   }
665   else {
666     /* The code supports any axis that is a combination of X,Y,Z
667      * although currently UI only allows to set the 3 different axis */
668     if (calc->smd->projAxis & MOD_SHRINKWRAP_PROJECT_OVER_X_AXIS) {
669       proj_axis[0] = 1.0f;
670     }
671     if (calc->smd->projAxis & MOD_SHRINKWRAP_PROJECT_OVER_Y_AXIS) {
672       proj_axis[1] = 1.0f;
673     }
674     if (calc->smd->projAxis & MOD_SHRINKWRAP_PROJECT_OVER_Z_AXIS) {
675       proj_axis[2] = 1.0f;
676     }
677
678     normalize_v3(proj_axis);
679
680     /* Invalid projection direction */
681     if (len_squared_v3(proj_axis) < FLT_EPSILON) {
682       return;
683     }
684   }
685
686   if (calc->aux_target) {
687     auxMesh = BKE_modifier_get_evaluated_mesh_from_evaluated_object(calc->aux_target, false);
688     if (!auxMesh) {
689       return;
690     }
691     BLI_SPACE_TRANSFORM_SETUP(&local2aux, calc->ob, calc->aux_target);
692   }
693
694   if (BKE_shrinkwrap_init_tree(
695           &aux_tree_stack, auxMesh, calc->smd->shrinkType, calc->smd->shrinkMode, false)) {
696     aux_tree = &aux_tree_stack;
697   }
698
699   /* After successfully build the trees, start projection vertices. */
700   ShrinkwrapCalcCBData data = {
701       .calc = calc,
702       .tree = calc->tree,
703       .aux_tree = aux_tree,
704       .proj_axis = proj_axis,
705       .local2aux = &local2aux,
706   };
707   TaskParallelSettings settings;
708   BLI_parallel_range_settings_defaults(&settings);
709   settings.use_threading = (calc->numVerts > BKE_MESH_OMP_LIMIT);
710   settings.userdata_chunk = &hit;
711   settings.userdata_chunk_size = sizeof(hit);
712   BLI_task_parallel_range(
713       0, calc->numVerts, &data, shrinkwrap_calc_normal_projection_cb_ex, &settings);
714
715   /* free data structures */
716   if (aux_tree) {
717     BKE_shrinkwrap_free_tree(aux_tree);
718   }
719 }
720
721 /*
722  * Shrinkwrap Target Surface Project mode
723  *
724  * It uses Newton's method to find a surface location with its
725  * smooth normal pointing at the original point.
726  *
727  * The equation system on barycentric weights and normal multiplier:
728  *
729  *   (w0*V0 + w1*V1 + w2*V2) + l * (w0*N0 + w1*N1 + w2*N2) - CO = 0
730  *   w0 + w1 + w2 = 1
731  *
732  * The actual solution vector is [ w0, w1, l ], with w2 eliminated.
733  */
734
735 //#define TRACE_TARGET_PROJECT
736
737 typedef struct TargetProjectTriData {
738   const float **vtri_co;
739   const float (*vtri_no)[3];
740   const float *point_co;
741
742   float n0_minus_n2[3], n1_minus_n2[3];
743   float c0_minus_c2[3], c1_minus_c2[3];
744
745   /* Current interpolated position and normal. */
746   float co_interp[3], no_interp[3];
747 } TargetProjectTriData;
748
749 /* Computes the deviation of the equation system from goal. */
750 static void target_project_tri_deviation(void *userdata, const float x[3], float r_delta[3])
751 {
752   TargetProjectTriData *data = userdata;
753
754   float w[3] = {x[0], x[1], 1.0f - x[0] - x[1]};
755   interp_v3_v3v3v3(data->co_interp, data->vtri_co[0], data->vtri_co[1], data->vtri_co[2], w);
756   interp_v3_v3v3v3(data->no_interp, data->vtri_no[0], data->vtri_no[1], data->vtri_no[2], w);
757
758   madd_v3_v3v3fl(r_delta, data->co_interp, data->no_interp, x[2]);
759   sub_v3_v3(r_delta, data->point_co);
760 }
761
762 /* Computes the Jacobian matrix of the equation system. */
763 static void target_project_tri_jacobian(void *userdata, const float x[3], float r_jacobian[3][3])
764 {
765   TargetProjectTriData *data = userdata;
766
767   madd_v3_v3v3fl(r_jacobian[0], data->c0_minus_c2, data->n0_minus_n2, x[2]);
768   madd_v3_v3v3fl(r_jacobian[1], data->c1_minus_c2, data->n1_minus_n2, x[2]);
769
770   copy_v3_v3(r_jacobian[2], data->vtri_no[2]);
771   madd_v3_v3fl(r_jacobian[2], data->n0_minus_n2, x[0]);
772   madd_v3_v3fl(r_jacobian[2], data->n1_minus_n2, x[1]);
773 }
774
775 /* Clamp barycentric weights to the triangle. */
776 static void target_project_tri_clamp(float x[3])
777 {
778   if (x[0] < 0.0f) {
779     x[0] = 0.0f;
780   }
781   if (x[1] < 0.0f) {
782     x[1] = 0.0f;
783   }
784   if (x[0] + x[1] > 1.0f) {
785     x[0] = x[0] / (x[0] + x[1]);
786     x[1] = 1.0f - x[0];
787   }
788 }
789
790 /* Correct the Newton's method step to keep the coordinates within the triangle. */
791 static bool target_project_tri_correct(void *UNUSED(userdata),
792                                        const float x[3],
793                                        float step[3],
794                                        float x_next[3])
795 {
796   /* Insignificant correction threshold */
797   const float epsilon = 1e-5f;
798   /* Dot product threshold for checking if step is 'clearly' pointing outside. */
799   const float dir_epsilon = 0.5f;
800   bool fixed = false, locked = false;
801
802   /* The barycentric coordinate domain is a triangle bounded by
803    * the X and Y axes, plus the x+y=1 diagonal. First, clamp the
804    * movement against the diagonal. Note that step is subtracted. */
805   float sum = x[0] + x[1];
806   float sstep = -(step[0] + step[1]);
807
808   if (sum + sstep > 1.0f) {
809     float ldist = 1.0f - sum;
810
811     /* If already at the boundary, slide along it. */
812     if (ldist < epsilon * (float)M_SQRT2) {
813       float step_len = len_v2(step);
814
815       /* Abort if the solution is clearly outside the domain. */
816       if (step_len > epsilon && sstep > step_len * dir_epsilon * (float)M_SQRT2) {
817         return false;
818       }
819
820       /* Project the new position onto the diagonal. */
821       add_v2_fl(step, (sum + sstep - 1.0f) * 0.5f);
822       fixed = locked = true;
823     }
824     else {
825       /* Scale a significant step down to arrive at the boundary. */
826       mul_v3_fl(step, ldist / sstep);
827       fixed = true;
828     }
829   }
830
831   /* Weight 0 and 1 boundary checks - along axis. */
832   for (int i = 0; i < 2; i++) {
833     if (step[i] > x[i]) {
834       /* If already at the boundary, slide along it. */
835       if (x[i] < epsilon) {
836         float step_len = len_v2(step);
837
838         /* Abort if the solution is clearly outside the domain. */
839         if (step_len > epsilon && (locked || step[i] > step_len * dir_epsilon)) {
840           return false;
841         }
842
843         /* Reset precision errors to stay at the boundary. */
844         step[i] = x[i];
845         fixed = true;
846       }
847       else {
848         /* Scale a significant step down to arrive at the boundary. */
849         mul_v3_fl(step, x[i] / step[i]);
850         fixed = true;
851       }
852     }
853   }
854
855   /* Recompute and clamp the new coordinates after step correction. */
856   if (fixed) {
857     sub_v3_v3v3(x_next, x, step);
858     target_project_tri_clamp(x_next);
859   }
860
861   return true;
862 }
863
864 static bool target_project_solve_point_tri(const float *vtri_co[3],
865                                            const float vtri_no[3][3],
866                                            const float point_co[3],
867                                            const float hit_co[3],
868                                            float hit_dist_sq,
869                                            float r_hit_co[3],
870                                            float r_hit_no[3])
871 {
872   float x[3], tmp[3];
873   float dist = sqrtf(hit_dist_sq);
874   float magnitude_estimate = dist + len_manhattan_v3(vtri_co[0]) + len_manhattan_v3(vtri_co[1]) +
875                              len_manhattan_v3(vtri_co[2]);
876   float epsilon = magnitude_estimate * 1.0e-6f;
877
878   /* Initial solution vector: barycentric weights plus distance along normal. */
879   interp_weights_tri_v3(x, UNPACK3(vtri_co), hit_co);
880
881   interp_v3_v3v3v3(r_hit_no, UNPACK3(vtri_no), x);
882   sub_v3_v3v3(tmp, point_co, hit_co);
883
884   x[2] = (dot_v3v3(tmp, r_hit_no) < 0) ? -dist : dist;
885
886   /* Solve the equations iteratively. */
887   TargetProjectTriData tri_data = {
888       .vtri_co = vtri_co,
889       .vtri_no = vtri_no,
890       .point_co = point_co,
891   };
892
893   sub_v3_v3v3(tri_data.n0_minus_n2, vtri_no[0], vtri_no[2]);
894   sub_v3_v3v3(tri_data.n1_minus_n2, vtri_no[1], vtri_no[2]);
895   sub_v3_v3v3(tri_data.c0_minus_c2, vtri_co[0], vtri_co[2]);
896   sub_v3_v3v3(tri_data.c1_minus_c2, vtri_co[1], vtri_co[2]);
897
898   target_project_tri_clamp(x);
899
900 #ifdef TRACE_TARGET_PROJECT
901   const bool trace = true;
902 #else
903   const bool trace = false;
904 #endif
905
906   bool ok = BLI_newton3d_solve(target_project_tri_deviation,
907                                target_project_tri_jacobian,
908                                target_project_tri_correct,
909                                &tri_data,
910                                epsilon,
911                                20,
912                                trace,
913                                x,
914                                x);
915
916   if (ok) {
917     copy_v3_v3(r_hit_co, tri_data.co_interp);
918     copy_v3_v3(r_hit_no, tri_data.no_interp);
919
920     return true;
921   }
922
923   return false;
924 }
925
926 static bool update_hit(BVHTreeNearest *nearest,
927                        int index,
928                        const float co[3],
929                        const float hit_co[3],
930                        const float hit_no[3])
931 {
932   float dist_sq = len_squared_v3v3(hit_co, co);
933
934   if (dist_sq < nearest->dist_sq) {
935 #ifdef TRACE_TARGET_PROJECT
936     printf(
937         "#=#=#> %d (%.3f,%.3f,%.3f) %g < %g\n", index, UNPACK3(hit_co), dist_sq, nearest->dist_sq);
938 #endif
939     nearest->index = index;
940     nearest->dist_sq = dist_sq;
941     copy_v3_v3(nearest->co, hit_co);
942     normalize_v3_v3(nearest->no, hit_no);
943     return true;
944   }
945
946   return false;
947 }
948
949 /* Target projection on a non-manifold boundary edge -
950  * treats it like an infinitely thin cylinder. */
951 static void target_project_edge(const ShrinkwrapTreeData *tree,
952                                 int index,
953                                 const float co[3],
954                                 BVHTreeNearest *nearest,
955                                 int eidx)
956 {
957   const BVHTreeFromMesh *data = &tree->treeData;
958   const MEdge *edge = &tree->mesh->medge[eidx];
959   const float *vedge_co[2] = {data->vert[edge->v1].co, data->vert[edge->v2].co};
960
961 #ifdef TRACE_TARGET_PROJECT
962   printf("EDGE %d (%.3f,%.3f,%.3f) (%.3f,%.3f,%.3f)\n",
963          eidx,
964          UNPACK3(vedge_co[0]),
965          UNPACK3(vedge_co[1]));
966 #endif
967
968   /* Retrieve boundary vertex IDs */
969   const int *vert_boundary_id = tree->boundary->vert_boundary_id;
970   int bid1 = vert_boundary_id[edge->v1], bid2 = vert_boundary_id[edge->v2];
971
972   if (bid1 < 0 || bid2 < 0) {
973     return;
974   }
975
976   /* Retrieve boundary vertex normals and align them to direction. */
977   const ShrinkwrapBoundaryVertData *boundary_verts = tree->boundary->boundary_verts;
978   float vedge_dir[2][3], dir[3];
979
980   copy_v3_v3(vedge_dir[0], boundary_verts[bid1].normal_plane);
981   copy_v3_v3(vedge_dir[1], boundary_verts[bid2].normal_plane);
982
983   sub_v3_v3v3(dir, vedge_co[1], vedge_co[0]);
984
985   if (dot_v3v3(boundary_verts[bid1].direction, dir) < 0) {
986     negate_v3(vedge_dir[0]);
987   }
988   if (dot_v3v3(boundary_verts[bid2].direction, dir) < 0) {
989     negate_v3(vedge_dir[1]);
990   }
991
992   /* Solve a quadratic equation: lerp(d0,d1,x) * (co - lerp(v0,v1,x)) = 0 */
993   float d0v0 = dot_v3v3(vedge_dir[0], vedge_co[0]), d0v1 = dot_v3v3(vedge_dir[0], vedge_co[1]);
994   float d1v0 = dot_v3v3(vedge_dir[1], vedge_co[0]), d1v1 = dot_v3v3(vedge_dir[1], vedge_co[1]);
995   float d0co = dot_v3v3(vedge_dir[0], co);
996
997   float a = d0v1 - d0v0 + d1v0 - d1v1;
998   float b = 2 * d0v0 - d0v1 - d0co - d1v0 + dot_v3v3(vedge_dir[1], co);
999   float c = d0co - d0v0;
1000   float det = b * b - 4 * a * c;
1001
1002   if (det >= 0) {
1003     const float epsilon = 1e-6f;
1004     float sdet = sqrtf(det);
1005     float hit_co[3], hit_no[3];
1006
1007     for (int i = (det > 0 ? 2 : 0); i >= 0; i -= 2) {
1008       float x = (-b + ((float)i - 1) * sdet) / (2 * a);
1009
1010       if (x >= -epsilon && x <= 1.0f + epsilon) {
1011         CLAMP(x, 0, 1);
1012
1013         float vedge_no[2][3];
1014         normal_short_to_float_v3(vedge_no[0], data->vert[edge->v1].no);
1015         normal_short_to_float_v3(vedge_no[1], data->vert[edge->v2].no);
1016
1017         interp_v3_v3v3(hit_co, vedge_co[0], vedge_co[1], x);
1018         interp_v3_v3v3(hit_no, vedge_no[0], vedge_no[1], x);
1019
1020         update_hit(nearest, index, co, hit_co, hit_no);
1021       }
1022     }
1023   }
1024 }
1025
1026 /* Target normal projection BVH callback - based on mesh_looptri_nearest_point. */
1027 static void mesh_looptri_target_project(void *userdata,
1028                                         int index,
1029                                         const float co[3],
1030                                         BVHTreeNearest *nearest)
1031 {
1032   const ShrinkwrapTreeData *tree = (ShrinkwrapTreeData *)userdata;
1033   const BVHTreeFromMesh *data = &tree->treeData;
1034   const MLoopTri *lt = &data->looptri[index];
1035   const MLoop *loop[3] = {
1036       &data->loop[lt->tri[0]], &data->loop[lt->tri[1]], &data->loop[lt->tri[2]]};
1037   const MVert *vtri[3] = {
1038       &data->vert[loop[0]->v], &data->vert[loop[1]->v], &data->vert[loop[2]->v]};
1039   const float *vtri_co[3] = {vtri[0]->co, vtri[1]->co, vtri[2]->co};
1040   float raw_hit_co[3], hit_co[3], hit_no[3], dist_sq, vtri_no[3][3];
1041
1042   /* First find the closest point and bail out if it's worse than the current solution. */
1043   closest_on_tri_to_point_v3(raw_hit_co, co, UNPACK3(vtri_co));
1044   dist_sq = len_squared_v3v3(co, raw_hit_co);
1045
1046 #ifdef TRACE_TARGET_PROJECT
1047   printf("TRIANGLE %d (%.3f,%.3f,%.3f) (%.3f,%.3f,%.3f) (%.3f,%.3f,%.3f) %g %g\n",
1048          index,
1049          UNPACK3(vtri_co[0]),
1050          UNPACK3(vtri_co[1]),
1051          UNPACK3(vtri_co[2]),
1052          dist_sq,
1053          nearest->dist_sq);
1054 #endif
1055
1056   if (dist_sq >= nearest->dist_sq) {
1057     return;
1058   }
1059
1060   /* Decode normals */
1061   normal_short_to_float_v3(vtri_no[0], vtri[0]->no);
1062   normal_short_to_float_v3(vtri_no[1], vtri[1]->no);
1063   normal_short_to_float_v3(vtri_no[2], vtri[2]->no);
1064
1065   /* Solve the equations for the triangle */
1066   if (target_project_solve_point_tri(vtri_co, vtri_no, co, raw_hit_co, dist_sq, hit_co, hit_no)) {
1067     update_hit(nearest, index, co, hit_co, hit_no);
1068   }
1069   /* Boundary edges */
1070   else if (tree->boundary && BLI_BITMAP_TEST(tree->boundary->looptri_has_boundary, index)) {
1071     const BLI_bitmap *is_boundary = tree->boundary->edge_is_boundary;
1072     int edges[3];
1073
1074     BKE_mesh_looptri_get_real_edges(tree->mesh, lt, edges);
1075
1076     for (int i = 0; i < 3; i++) {
1077       if (edges[i] >= 0 && BLI_BITMAP_TEST(is_boundary, edges[i])) {
1078         target_project_edge(tree, index, co, nearest, edges[i]);
1079       }
1080     }
1081   }
1082 }
1083
1084 /*
1085  * Maps the point to the nearest surface, either by simple nearest, or by target normal projection.
1086  */
1087 void BKE_shrinkwrap_find_nearest_surface(struct ShrinkwrapTreeData *tree,
1088                                          BVHTreeNearest *nearest,
1089                                          float co[3],
1090                                          int type)
1091 {
1092   BVHTreeFromMesh *treeData = &tree->treeData;
1093
1094   if (type == MOD_SHRINKWRAP_TARGET_PROJECT) {
1095 #ifdef TRACE_TARGET_PROJECT
1096     printf("\n====== TARGET PROJECT START ======\n");
1097 #endif
1098
1099     BLI_bvhtree_find_nearest_ex(
1100         tree->bvh, co, nearest, mesh_looptri_target_project, tree, BVH_NEAREST_OPTIMAL_ORDER);
1101
1102 #ifdef TRACE_TARGET_PROJECT
1103     printf("====== TARGET PROJECT END: %d %g ======\n\n", nearest->index, nearest->dist_sq);
1104 #endif
1105
1106     if (nearest->index < 0) {
1107       /* fallback to simple nearest */
1108       BLI_bvhtree_find_nearest(tree->bvh, co, nearest, treeData->nearest_callback, treeData);
1109     }
1110   }
1111   else {
1112     BLI_bvhtree_find_nearest(tree->bvh, co, nearest, treeData->nearest_callback, treeData);
1113   }
1114 }
1115
1116 /*
1117  * Shrinkwrap moving vertexs to the nearest surface point on the target
1118  *
1119  * it builds a BVHTree from the target mesh and then performs a
1120  * NN matches for each vertex
1121  */
1122 static void shrinkwrap_calc_nearest_surface_point_cb_ex(void *__restrict userdata,
1123                                                         const int i,
1124                                                         const TaskParallelTLS *__restrict tls)
1125 {
1126   ShrinkwrapCalcCBData *data = userdata;
1127
1128   ShrinkwrapCalcData *calc = data->calc;
1129   BVHTreeNearest *nearest = tls->userdata_chunk;
1130
1131   float *co = calc->vertexCos[i];
1132   float tmp_co[3];
1133   float weight = BKE_defvert_array_find_weight_safe(calc->dvert, i, calc->vgroup);
1134
1135   if (calc->invert_vgroup) {
1136     weight = 1.0f - weight;
1137   }
1138
1139   if (weight == 0.0f) {
1140     return;
1141   }
1142
1143   /* Convert the vertex to tree coordinates */
1144   if (calc->vert) {
1145     copy_v3_v3(tmp_co, calc->vert[i].co);
1146   }
1147   else {
1148     copy_v3_v3(tmp_co, co);
1149   }
1150   BLI_space_transform_apply(&calc->local2target, tmp_co);
1151
1152   /* Use local proximity heuristics (to reduce the nearest search)
1153    *
1154    * If we already had an hit before.. we assume this vertex is going to have a close hit to that
1155    * other vertex so we can initiate the "nearest.dist" with the expected value to that last hit.
1156    * This will lead in pruning of the search tree. */
1157   if (nearest->index != -1) {
1158     if (calc->smd->shrinkType == MOD_SHRINKWRAP_TARGET_PROJECT) {
1159       /* Heuristic doesn't work because of additional restrictions. */
1160       nearest->index = -1;
1161       nearest->dist_sq = FLT_MAX;
1162     }
1163     else {
1164       nearest->dist_sq = len_squared_v3v3(tmp_co, nearest->co);
1165     }
1166   }
1167   else {
1168     nearest->dist_sq = FLT_MAX;
1169   }
1170
1171   BKE_shrinkwrap_find_nearest_surface(data->tree, nearest, tmp_co, calc->smd->shrinkType);
1172
1173   /* Found the nearest vertex */
1174   if (nearest->index != -1) {
1175     BKE_shrinkwrap_snap_point_to_surface(data->tree,
1176                                          NULL,
1177                                          calc->smd->shrinkMode,
1178                                          nearest->index,
1179                                          nearest->co,
1180                                          nearest->no,
1181                                          calc->keepDist,
1182                                          tmp_co,
1183                                          tmp_co);
1184
1185     /* Convert the coordinates back to mesh coordinates */
1186     BLI_space_transform_invert(&calc->local2target, tmp_co);
1187     interp_v3_v3v3(co, co, tmp_co, weight); /* linear interpolation */
1188   }
1189 }
1190
1191 /**
1192  * Compute a smooth normal of the target (if applicable) at the hit location.
1193  *
1194  * \param tree: information about the mesh
1195  * \param transform: transform from the hit coordinate space to the object space; may be null
1196  * \param r_no: output in hit coordinate space; may be shared with inputs
1197  */
1198 void BKE_shrinkwrap_compute_smooth_normal(const struct ShrinkwrapTreeData *tree,
1199                                           const struct SpaceTransform *transform,
1200                                           int looptri_idx,
1201                                           const float hit_co[3],
1202                                           const float hit_no[3],
1203                                           float r_no[3])
1204 {
1205   const BVHTreeFromMesh *treeData = &tree->treeData;
1206   const MLoopTri *tri = &treeData->looptri[looptri_idx];
1207
1208   /* Interpolate smooth normals if enabled. */
1209   if ((tree->mesh->mpoly[tri->poly].flag & ME_SMOOTH) != 0) {
1210     const MVert *verts[] = {
1211         &treeData->vert[treeData->loop[tri->tri[0]].v],
1212         &treeData->vert[treeData->loop[tri->tri[1]].v],
1213         &treeData->vert[treeData->loop[tri->tri[2]].v],
1214     };
1215     float w[3], no[3][3], tmp_co[3];
1216
1217     /* Custom and auto smooth split normals. */
1218     if (tree->clnors) {
1219       copy_v3_v3(no[0], tree->clnors[tri->tri[0]]);
1220       copy_v3_v3(no[1], tree->clnors[tri->tri[1]]);
1221       copy_v3_v3(no[2], tree->clnors[tri->tri[2]]);
1222     }
1223     /* Ordinary vertex normals. */
1224     else {
1225       normal_short_to_float_v3(no[0], verts[0]->no);
1226       normal_short_to_float_v3(no[1], verts[1]->no);
1227       normal_short_to_float_v3(no[2], verts[2]->no);
1228     }
1229
1230     /* Barycentric weights from hit point. */
1231     copy_v3_v3(tmp_co, hit_co);
1232
1233     if (transform) {
1234       BLI_space_transform_apply(transform, tmp_co);
1235     }
1236
1237     interp_weights_tri_v3(w, verts[0]->co, verts[1]->co, verts[2]->co, tmp_co);
1238
1239     /* Interpolate using weights. */
1240     interp_v3_v3v3v3(r_no, no[0], no[1], no[2], w);
1241
1242     if (transform) {
1243       BLI_space_transform_invert_normal(transform, r_no);
1244     }
1245     else {
1246       normalize_v3(r_no);
1247     }
1248   }
1249   /* Use the polygon normal if flat. */
1250   else if (tree->pnors != NULL) {
1251     copy_v3_v3(r_no, tree->pnors[tri->poly]);
1252   }
1253   /* Finally fallback to the looptri normal. */
1254   else {
1255     copy_v3_v3(r_no, hit_no);
1256   }
1257 }
1258
1259 /* Helper for MOD_SHRINKWRAP_INSIDE, MOD_SHRINKWRAP_OUTSIDE and MOD_SHRINKWRAP_OUTSIDE_SURFACE. */
1260 static void shrinkwrap_snap_with_side(float r_point_co[3],
1261                                       const float point_co[3],
1262                                       const float hit_co[3],
1263                                       const float hit_no[3],
1264                                       float goal_dist,
1265                                       float forcesign,
1266                                       bool forcesnap)
1267 {
1268   float delta[3];
1269   sub_v3_v3v3(delta, point_co, hit_co);
1270
1271   float dist = len_v3(delta);
1272
1273   /* If exactly on the surface, push out along normal */
1274   if (dist < FLT_EPSILON) {
1275     if (forcesnap || goal_dist > 0) {
1276       madd_v3_v3v3fl(r_point_co, hit_co, hit_no, goal_dist * forcesign);
1277     }
1278     else {
1279       copy_v3_v3(r_point_co, hit_co);
1280     }
1281   }
1282   /* Move to the correct side if needed */
1283   else {
1284     float dsign = signf(dot_v3v3(delta, hit_no));
1285
1286     if (forcesign == 0.0f) {
1287       forcesign = dsign;
1288     }
1289
1290     /* If on the wrong side or too close, move to correct */
1291     if (forcesnap || dsign * dist * forcesign < goal_dist) {
1292       mul_v3_fl(delta, dsign / dist);
1293
1294       /* At very small distance, blend in the hit normal to stabilize math. */
1295       float dist_epsilon = (fabsf(goal_dist) + len_manhattan_v3(hit_co)) * 1e-4f;
1296
1297       if (dist < dist_epsilon) {
1298 #ifdef TRACE_TARGET_PROJECT
1299         printf("zero_factor %g = %g / %g\n", dist / dist_epsilon, dist, dist_epsilon);
1300 #endif
1301
1302         interp_v3_v3v3(delta, hit_no, delta, dist / dist_epsilon);
1303       }
1304
1305       madd_v3_v3v3fl(r_point_co, hit_co, delta, goal_dist * forcesign);
1306     }
1307     else {
1308       copy_v3_v3(r_point_co, point_co);
1309     }
1310   }
1311 }
1312
1313 /**
1314  * Apply the shrink to surface modes to the given original coordinates and nearest point.
1315  *
1316  * \param tree: mesh data for smooth normals
1317  * \param transform: transform from the hit coordinate space to the object space; may be null
1318  * \param r_point_co: may be the same memory location as point_co, hit_co, or hit_no.
1319  */
1320 void BKE_shrinkwrap_snap_point_to_surface(const struct ShrinkwrapTreeData *tree,
1321                                           const struct SpaceTransform *transform,
1322                                           int mode,
1323                                           int hit_idx,
1324                                           const float hit_co[3],
1325                                           const float hit_no[3],
1326                                           float goal_dist,
1327                                           const float point_co[3],
1328                                           float r_point_co[3])
1329 {
1330   float tmp[3];
1331
1332   switch (mode) {
1333     /* Offsets along the line between point_co and hit_co. */
1334     case MOD_SHRINKWRAP_ON_SURFACE:
1335       if (goal_dist != 0) {
1336         shrinkwrap_snap_with_side(r_point_co, point_co, hit_co, hit_no, goal_dist, 0, true);
1337       }
1338       else {
1339         copy_v3_v3(r_point_co, hit_co);
1340       }
1341       break;
1342
1343     case MOD_SHRINKWRAP_INSIDE:
1344       shrinkwrap_snap_with_side(r_point_co, point_co, hit_co, hit_no, goal_dist, -1, false);
1345       break;
1346
1347     case MOD_SHRINKWRAP_OUTSIDE:
1348       shrinkwrap_snap_with_side(r_point_co, point_co, hit_co, hit_no, goal_dist, +1, false);
1349       break;
1350
1351     case MOD_SHRINKWRAP_OUTSIDE_SURFACE:
1352       if (goal_dist != 0) {
1353         shrinkwrap_snap_with_side(r_point_co, point_co, hit_co, hit_no, goal_dist, +1, true);
1354       }
1355       else {
1356         copy_v3_v3(r_point_co, hit_co);
1357       }
1358       break;
1359
1360     /* Offsets along the normal */
1361     case MOD_SHRINKWRAP_ABOVE_SURFACE:
1362       if (goal_dist != 0) {
1363         BKE_shrinkwrap_compute_smooth_normal(tree, transform, hit_idx, hit_co, hit_no, tmp);
1364         madd_v3_v3v3fl(r_point_co, hit_co, tmp, goal_dist);
1365       }
1366       else {
1367         copy_v3_v3(r_point_co, hit_co);
1368       }
1369       break;
1370
1371     default:
1372       printf("Unknown Shrinkwrap surface snap mode: %d\n", mode);
1373       copy_v3_v3(r_point_co, hit_co);
1374   }
1375 }
1376
1377 static void shrinkwrap_calc_nearest_surface_point(ShrinkwrapCalcData *calc)
1378 {
1379   BVHTreeNearest nearest = NULL_BVHTreeNearest;
1380
1381   /* Setup nearest */
1382   nearest.index = -1;
1383   nearest.dist_sq = FLT_MAX;
1384
1385   /* Find the nearest vertex */
1386   ShrinkwrapCalcCBData data = {
1387       .calc = calc,
1388       .tree = calc->tree,
1389   };
1390   TaskParallelSettings settings;
1391   BLI_parallel_range_settings_defaults(&settings);
1392   settings.use_threading = (calc->numVerts > BKE_MESH_OMP_LIMIT);
1393   settings.userdata_chunk = &nearest;
1394   settings.userdata_chunk_size = sizeof(nearest);
1395   BLI_task_parallel_range(
1396       0, calc->numVerts, &data, shrinkwrap_calc_nearest_surface_point_cb_ex, &settings);
1397 }
1398
1399 /* Main shrinkwrap function */
1400 void shrinkwrapModifier_deform(ShrinkwrapModifierData *smd,
1401                                const ModifierEvalContext *ctx,
1402                                struct Scene *scene,
1403                                Object *ob,
1404                                Mesh *mesh,
1405                                MDeformVert *dvert,
1406                                const int defgrp_index,
1407                                float (*vertexCos)[3],
1408                                int numVerts)
1409 {
1410
1411   DerivedMesh *ss_mesh = NULL;
1412   ShrinkwrapCalcData calc = NULL_ShrinkwrapCalcData;
1413
1414   /* remove loop dependencies on derived meshes (TODO should this be done elsewhere?) */
1415   if (smd->target == ob) {
1416     smd->target = NULL;
1417   }
1418   if (smd->auxTarget == ob) {
1419     smd->auxTarget = NULL;
1420   }
1421
1422   /* Configure Shrinkwrap calc data */
1423   calc.smd = smd;
1424   calc.ob = ob;
1425   calc.numVerts = numVerts;
1426   calc.vertexCos = vertexCos;
1427   calc.dvert = dvert;
1428   calc.vgroup = defgrp_index;
1429   calc.invert_vgroup = (smd->shrinkOpts & MOD_SHRINKWRAP_INVERT_VGROUP) != 0;
1430
1431   if (smd->target != NULL) {
1432     Object *ob_target = DEG_get_evaluated_object(ctx->depsgraph, smd->target);
1433     calc.target = BKE_modifier_get_evaluated_mesh_from_evaluated_object(ob_target, false);
1434
1435     /* TODO there might be several "bugs" with non-uniform scales matrices
1436      * because it will no longer be nearest surface, not sphere projection
1437      * because space has been deformed */
1438     BLI_SPACE_TRANSFORM_SETUP(&calc.local2target, ob, ob_target);
1439
1440     /* TODO: smd->keepDist is in global units.. must change to local */
1441     calc.keepDist = smd->keepDist;
1442   }
1443   calc.aux_target = DEG_get_evaluated_object(ctx->depsgraph, smd->auxTarget);
1444
1445   if (mesh != NULL && smd->shrinkType == MOD_SHRINKWRAP_PROJECT) {
1446     /* Setup arrays to get vertexs positions, normals and deform weights */
1447     calc.vert = mesh->mvert;
1448
1449     /* Using vertexs positions/normals as if a subsurface was applied */
1450     if (smd->subsurfLevels) {
1451       SubsurfModifierData ssmd = {{NULL}};
1452       ssmd.subdivType = ME_CC_SUBSURF;  /* catmull clark */
1453       ssmd.levels = smd->subsurfLevels; /* levels */
1454
1455       /* TODO to be moved to Mesh once we are done with changes in subsurf code. */
1456       DerivedMesh *dm = CDDM_from_mesh(mesh);
1457
1458       ss_mesh = subsurf_make_derived_from_derived(
1459           dm, &ssmd, scene, NULL, (ob->mode & OB_MODE_EDIT) ? SUBSURF_IN_EDIT_MODE : 0);
1460
1461       if (ss_mesh) {
1462         calc.vert = ss_mesh->getVertDataArray(ss_mesh, CD_MVERT);
1463         if (calc.vert) {
1464           /* TRICKY: this code assumes subsurface will have the transformed original vertices
1465            * in their original order at the end of the vert array. */
1466           calc.vert = calc.vert + ss_mesh->getNumVerts(ss_mesh) - dm->getNumVerts(dm);
1467         }
1468       }
1469
1470       /* Just to make sure we are not leaving any memory behind */
1471       BLI_assert(ssmd.emCache == NULL);
1472       BLI_assert(ssmd.mCache == NULL);
1473
1474       dm->release(dm);
1475     }
1476   }
1477
1478   /* Projecting target defined - lets work! */
1479   ShrinkwrapTreeData tree;
1480
1481   if (BKE_shrinkwrap_init_tree(&tree, calc.target, smd->shrinkType, smd->shrinkMode, false)) {
1482     calc.tree = &tree;
1483
1484     switch (smd->shrinkType) {
1485       case MOD_SHRINKWRAP_NEAREST_SURFACE:
1486       case MOD_SHRINKWRAP_TARGET_PROJECT:
1487         TIMEIT_BENCH(shrinkwrap_calc_nearest_surface_point(&calc), deform_surface);
1488         break;
1489
1490       case MOD_SHRINKWRAP_PROJECT:
1491         TIMEIT_BENCH(shrinkwrap_calc_normal_projection(&calc), deform_project);
1492         break;
1493
1494       case MOD_SHRINKWRAP_NEAREST_VERTEX:
1495         TIMEIT_BENCH(shrinkwrap_calc_nearest_vertex(&calc), deform_vertex);
1496         break;
1497     }
1498
1499     BKE_shrinkwrap_free_tree(&tree);
1500   }
1501
1502   /* free memory */
1503   if (ss_mesh) {
1504     ss_mesh->release(ss_mesh);
1505   }
1506 }
1507
1508 void BKE_shrinkwrap_mesh_nearest_surface_deform(struct bContext *C,
1509                                                 Object *ob_source,
1510                                                 Object *ob_target)
1511 {
1512   Depsgraph *depsgraph = CTX_data_depsgraph_pointer(C);
1513   struct Scene *sce = CTX_data_scene(C);
1514   ShrinkwrapModifierData ssmd = {{0}};
1515   ModifierEvalContext ctx = {depsgraph, ob_source, 0};
1516   int totvert;
1517
1518   ssmd.target = ob_target;
1519   ssmd.shrinkType = MOD_SHRINKWRAP_NEAREST_SURFACE;
1520   ssmd.shrinkMode = MOD_SHRINKWRAP_ON_SURFACE;
1521   ssmd.keepDist = 0.0f;
1522
1523   Mesh *src_me = ob_source->data;
1524   float(*vertexCos)[3] = BKE_mesh_vert_coords_alloc(src_me, &totvert);
1525
1526   shrinkwrapModifier_deform(&ssmd, &ctx, sce, ob_source, src_me, NULL, -1, vertexCos, totvert);
1527
1528   BKE_mesh_vert_coords_apply(src_me, vertexCos);
1529
1530   MEM_freeN(vertexCos);
1531 }
1532
1533 void BKE_shrinkwrap_remesh_target_project(Mesh *src_me, Mesh *target_me, Object *ob_target)
1534 {
1535   ShrinkwrapModifierData ssmd = {{0}};
1536   int totvert;
1537
1538   ssmd.target = ob_target;
1539   ssmd.shrinkType = MOD_SHRINKWRAP_PROJECT;
1540   ssmd.shrinkMode = MOD_SHRINKWRAP_ON_SURFACE;
1541   ssmd.shrinkOpts = MOD_SHRINKWRAP_PROJECT_ALLOW_NEG_DIR | MOD_SHRINKWRAP_PROJECT_ALLOW_POS_DIR;
1542   ssmd.keepDist = 0.0f;
1543   ssmd.projLimit = target_me->remesh_voxel_size;
1544
1545   float(*vertexCos)[3] = BKE_mesh_vert_coords_alloc(src_me, &totvert);
1546
1547   ShrinkwrapCalcData calc = NULL_ShrinkwrapCalcData;
1548
1549   calc.smd = &ssmd;
1550   calc.numVerts = src_me->totvert;
1551   calc.vertexCos = vertexCos;
1552   calc.vgroup = -1;
1553   calc.target = target_me;
1554   calc.keepDist = ssmd.keepDist;
1555   calc.vert = src_me->mvert;
1556   BLI_SPACE_TRANSFORM_SETUP(&calc.local2target, ob_target, ob_target);
1557
1558   ShrinkwrapTreeData tree;
1559   if (BKE_shrinkwrap_init_tree(&tree, calc.target, ssmd.shrinkType, ssmd.shrinkMode, false)) {
1560     calc.tree = &tree;
1561     TIMEIT_BENCH(shrinkwrap_calc_normal_projection(&calc), deform_project);
1562     BKE_shrinkwrap_free_tree(&tree);
1563   }
1564
1565   BKE_mesh_vert_coords_apply(src_me, vertexCos);
1566
1567   MEM_freeN(vertexCos);
1568 }