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