ClangFormat: apply to source, most of intern
[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   else
384     nearest->dist_sq = FLT_MAX;
385
386   BLI_bvhtree_find_nearest(treeData->tree, tmp_co, nearest, treeData->nearest_callback, treeData);
387
388   /* Found the nearest vertex */
389   if (nearest->index != -1) {
390     /* Adjusting the vertex weight,
391      * so that after interpolating it keeps a certain distance from the nearest position */
392     if (nearest->dist_sq > FLT_EPSILON) {
393       const float dist = sqrtf(nearest->dist_sq);
394       weight *= (dist - calc->keepDist) / dist;
395     }
396
397     /* Convert the coordinates back to mesh coordinates */
398     copy_v3_v3(tmp_co, nearest->co);
399     BLI_space_transform_invert(&calc->local2target, tmp_co);
400
401     interp_v3_v3v3(co, co, tmp_co, weight); /* linear interpolation */
402   }
403 }
404
405 static void shrinkwrap_calc_nearest_vertex(ShrinkwrapCalcData *calc)
406 {
407   BVHTreeNearest nearest = NULL_BVHTreeNearest;
408
409   /* Setup nearest */
410   nearest.index = -1;
411   nearest.dist_sq = FLT_MAX;
412
413   ShrinkwrapCalcCBData data = {
414       .calc = calc,
415       .tree = calc->tree,
416   };
417   ParallelRangeSettings settings;
418   BLI_parallel_range_settings_defaults(&settings);
419   settings.use_threading = (calc->numVerts > BKE_MESH_OMP_LIMIT);
420   settings.userdata_chunk = &nearest;
421   settings.userdata_chunk_size = sizeof(nearest);
422   BLI_task_parallel_range(
423       0, calc->numVerts, &data, shrinkwrap_calc_nearest_vertex_cb_ex, &settings);
424 }
425
426 /*
427  * This function raycast a single vertex and updates the hit if the "hit" is considered valid.
428  * Returns true if "hit" was updated.
429  * Opts control whether an hit is valid or not
430  * Supported options are:
431  * - MOD_SHRINKWRAP_CULL_TARGET_FRONTFACE (front faces hits are ignored)
432  * - MOD_SHRINKWRAP_CULL_TARGET_BACKFACE (back faces hits are ignored)
433  */
434 bool BKE_shrinkwrap_project_normal(char options,
435                                    const float vert[3],
436                                    const float dir[3],
437                                    const float ray_radius,
438                                    const SpaceTransform *transf,
439                                    ShrinkwrapTreeData *tree,
440                                    BVHTreeRayHit *hit)
441 {
442   /* don't use this because this dist value could be incompatible
443    * this value used by the callback for comparing prev/new dist values.
444    * also, at the moment there is no need to have a corrected 'dist' value */
445   // #define USE_DIST_CORRECT
446
447   float tmp_co[3], tmp_no[3];
448   const float *co, *no;
449   BVHTreeRayHit hit_tmp;
450
451   /* Copy from hit (we need to convert hit rays from one space coordinates to the other */
452   memcpy(&hit_tmp, hit, sizeof(hit_tmp));
453
454   /* Apply space transform (TODO readjust dist) */
455   if (transf) {
456     copy_v3_v3(tmp_co, vert);
457     BLI_space_transform_apply(transf, tmp_co);
458     co = tmp_co;
459
460     copy_v3_v3(tmp_no, dir);
461     BLI_space_transform_apply_normal(transf, tmp_no);
462     no = tmp_no;
463
464 #ifdef USE_DIST_CORRECT
465     hit_tmp.dist *= mat4_to_scale(((SpaceTransform *)transf)->local2target);
466 #endif
467   }
468   else {
469     co = vert;
470     no = dir;
471   }
472
473   hit_tmp.index = -1;
474
475   BLI_bvhtree_ray_cast(
476       tree->bvh, co, no, ray_radius, &hit_tmp, tree->treeData.raycast_callback, &tree->treeData);
477
478   if (hit_tmp.index != -1) {
479     /* invert the normal first so face culling works on rotated objects */
480     if (transf) {
481       BLI_space_transform_invert_normal(transf, hit_tmp.no);
482     }
483
484     if (options & MOD_SHRINKWRAP_CULL_TARGET_MASK) {
485       /* apply backface */
486       const float dot = dot_v3v3(dir, hit_tmp.no);
487       if (((options & MOD_SHRINKWRAP_CULL_TARGET_FRONTFACE) && dot <= 0.0f) ||
488           ((options & MOD_SHRINKWRAP_CULL_TARGET_BACKFACE) && dot >= 0.0f)) {
489         return false; /* Ignore hit */
490       }
491     }
492
493     if (transf) {
494       /* Inverting space transform (TODO make coeherent with the initial dist readjust) */
495       BLI_space_transform_invert(transf, hit_tmp.co);
496 #ifdef USE_DIST_CORRECT
497       hit_tmp.dist = len_v3v3(vert, hit_tmp.co);
498 #endif
499     }
500
501     BLI_assert(hit_tmp.dist <= hit->dist);
502
503     memcpy(hit, &hit_tmp, sizeof(hit_tmp));
504     return true;
505   }
506   return false;
507 }
508
509 static void shrinkwrap_calc_normal_projection_cb_ex(void *__restrict userdata,
510                                                     const int i,
511                                                     const ParallelRangeTLS *__restrict tls)
512 {
513   ShrinkwrapCalcCBData *data = userdata;
514
515   ShrinkwrapCalcData *calc = data->calc;
516   ShrinkwrapTreeData *tree = data->tree;
517   ShrinkwrapTreeData *aux_tree = data->aux_tree;
518
519   float *proj_axis = data->proj_axis;
520   SpaceTransform *local2aux = data->local2aux;
521
522   BVHTreeRayHit *hit = tls->userdata_chunk;
523
524   const float proj_limit_squared = calc->smd->projLimit * calc->smd->projLimit;
525   float *co = calc->vertexCos[i];
526   float tmp_co[3], tmp_no[3];
527   float weight = defvert_array_find_weight_safe(calc->dvert, i, calc->vgroup);
528
529   if (calc->invert_vgroup) {
530     weight = 1.0f - weight;
531   }
532
533   if (weight == 0.0f) {
534     return;
535   }
536
537   if (calc->vert != NULL && calc->smd->projAxis == MOD_SHRINKWRAP_PROJECT_OVER_NORMAL) {
538     /* calc->vert contains verts from evaluated mesh.  */
539     /* These coordinates are deformed by vertexCos only for normal projection (to get correct normals) */
540     /* for other cases calc->verts contains undeformed coordinates and vertexCos should be used */
541     copy_v3_v3(tmp_co, calc->vert[i].co);
542     normal_short_to_float_v3(tmp_no, calc->vert[i].no);
543   }
544   else {
545     copy_v3_v3(tmp_co, co);
546     copy_v3_v3(tmp_no, proj_axis);
547   }
548
549   hit->index = -1;
550   hit->dist =
551       BVH_RAYCAST_DIST_MAX; /* TODO: we should use FLT_MAX here, but sweepsphere code isn't prepared for that */
552
553   bool is_aux = false;
554
555   /* Project over positive direction of axis */
556   if (calc->smd->shrinkOpts & MOD_SHRINKWRAP_PROJECT_ALLOW_POS_DIR) {
557     if (aux_tree) {
558       if (BKE_shrinkwrap_project_normal(0, tmp_co, tmp_no, 0.0, local2aux, aux_tree, hit)) {
559         is_aux = true;
560       }
561     }
562
563     if (BKE_shrinkwrap_project_normal(
564             calc->smd->shrinkOpts, tmp_co, tmp_no, 0.0, &calc->local2target, tree, hit)) {
565       is_aux = false;
566     }
567   }
568
569   /* Project over negative direction of axis */
570   if (calc->smd->shrinkOpts & MOD_SHRINKWRAP_PROJECT_ALLOW_NEG_DIR) {
571     float inv_no[3];
572     negate_v3_v3(inv_no, tmp_no);
573
574     char options = calc->smd->shrinkOpts;
575
576     if ((options & MOD_SHRINKWRAP_INVERT_CULL_TARGET) &&
577         (options & MOD_SHRINKWRAP_CULL_TARGET_MASK)) {
578       options ^= MOD_SHRINKWRAP_CULL_TARGET_MASK;
579     }
580
581     if (aux_tree) {
582       if (BKE_shrinkwrap_project_normal(0, tmp_co, inv_no, 0.0, local2aux, aux_tree, hit)) {
583         is_aux = true;
584       }
585     }
586
587     if (BKE_shrinkwrap_project_normal(
588             options, tmp_co, inv_no, 0.0, &calc->local2target, tree, hit)) {
589       is_aux = false;
590     }
591   }
592
593   /* don't set the initial dist (which is more efficient),
594    * because its calculated in the targets space, we want the dist in our own space */
595   if (proj_limit_squared != 0.0f) {
596     if (hit->index != -1 && len_squared_v3v3(hit->co, co) > proj_limit_squared) {
597       hit->index = -1;
598     }
599   }
600
601   if (hit->index != -1) {
602     if (is_aux) {
603       BKE_shrinkwrap_snap_point_to_surface(aux_tree,
604                                            local2aux,
605                                            calc->smd->shrinkMode,
606                                            hit->index,
607                                            hit->co,
608                                            hit->no,
609                                            calc->keepDist,
610                                            tmp_co,
611                                            hit->co);
612     }
613     else {
614       BKE_shrinkwrap_snap_point_to_surface(tree,
615                                            &calc->local2target,
616                                            calc->smd->shrinkMode,
617                                            hit->index,
618                                            hit->co,
619                                            hit->no,
620                                            calc->keepDist,
621                                            tmp_co,
622                                            hit->co);
623     }
624
625     interp_v3_v3v3(co, co, hit->co, weight);
626   }
627 }
628
629 static void shrinkwrap_calc_normal_projection(ShrinkwrapCalcData *calc)
630 {
631   /* Options about projection direction */
632   float proj_axis[3] = {0.0f, 0.0f, 0.0f};
633
634   /* Raycast and tree stuff */
635
636   /** \note 'hit.dist' is kept in the targets space, this is only used
637    * for finding the best hit, to get the real dist,
638    * measure the len_v3v3() from the input coord to hit.co */
639   BVHTreeRayHit hit;
640
641   /* auxiliary target */
642   Mesh *auxMesh = NULL;
643   ShrinkwrapTreeData *aux_tree = NULL;
644   ShrinkwrapTreeData aux_tree_stack;
645   SpaceTransform local2aux;
646
647   /* If the user doesn't allows to project in any direction of projection axis
648    * then there's nothing todo. */
649   if ((calc->smd->shrinkOpts &
650        (MOD_SHRINKWRAP_PROJECT_ALLOW_POS_DIR | MOD_SHRINKWRAP_PROJECT_ALLOW_NEG_DIR)) == 0)
651     return;
652
653   /* Prepare data to retrieve the direction in which we should project each vertex */
654   if (calc->smd->projAxis == MOD_SHRINKWRAP_PROJECT_OVER_NORMAL) {
655     if (calc->vert == NULL)
656       return;
657   }
658   else {
659     /* The code supports any axis that is a combination of X,Y,Z
660      * although currently UI only allows to set the 3 different axis */
661     if (calc->smd->projAxis & MOD_SHRINKWRAP_PROJECT_OVER_X_AXIS)
662       proj_axis[0] = 1.0f;
663     if (calc->smd->projAxis & MOD_SHRINKWRAP_PROJECT_OVER_Y_AXIS)
664       proj_axis[1] = 1.0f;
665     if (calc->smd->projAxis & MOD_SHRINKWRAP_PROJECT_OVER_Z_AXIS)
666       proj_axis[2] = 1.0f;
667
668     normalize_v3(proj_axis);
669
670     /* Invalid projection direction */
671     if (len_squared_v3(proj_axis) < FLT_EPSILON) {
672       return;
673     }
674   }
675
676   if (calc->aux_target) {
677     auxMesh = BKE_modifier_get_evaluated_mesh_from_evaluated_object(calc->aux_target, false);
678     if (!auxMesh)
679       return;
680     BLI_SPACE_TRANSFORM_SETUP(&local2aux, calc->ob, calc->aux_target);
681   }
682
683   if (BKE_shrinkwrap_init_tree(
684           &aux_tree_stack, auxMesh, calc->smd->shrinkType, calc->smd->shrinkMode, false)) {
685     aux_tree = &aux_tree_stack;
686   }
687
688   /* After successfully build the trees, start projection vertices. */
689   ShrinkwrapCalcCBData data = {
690       .calc = calc,
691       .tree = calc->tree,
692       .aux_tree = aux_tree,
693       .proj_axis = proj_axis,
694       .local2aux = &local2aux,
695   };
696   ParallelRangeSettings settings;
697   BLI_parallel_range_settings_defaults(&settings);
698   settings.use_threading = (calc->numVerts > BKE_MESH_OMP_LIMIT);
699   settings.userdata_chunk = &hit;
700   settings.userdata_chunk_size = sizeof(hit);
701   BLI_task_parallel_range(
702       0, calc->numVerts, &data, shrinkwrap_calc_normal_projection_cb_ex, &settings);
703
704   /* free data structures */
705   if (aux_tree) {
706     BKE_shrinkwrap_free_tree(aux_tree);
707   }
708 }
709
710 /*
711  * Shrinkwrap Target Surface Project mode
712  *
713  * It uses Newton's method to find a surface location with its
714  * smooth normal pointing at the original point.
715  *
716  * The equation system on barycentric weights and normal multiplier:
717  *
718  *   (w0*V0 + w1*V1 + w2*V2) + l * (w0*N0 + w1*N1 + w2*N2) - CO = 0
719  *   w0 + w1 + w2 = 1
720  *
721  * The actual solution vector is [ w0, w1, l ], with w2 eliminated.
722  */
723
724 //#define TRACE_TARGET_PROJECT
725
726 typedef struct TargetProjectTriData {
727   const float **vtri_co;
728   const float (*vtri_no)[3];
729   const float *point_co;
730
731   float n0_minus_n2[3], n1_minus_n2[3];
732   float c0_minus_c2[3], c1_minus_c2[3];
733
734   /* Current interpolated position and normal. */
735   float co_interp[3], no_interp[3];
736 } TargetProjectTriData;
737
738 /* Computes the deviation of the equation system from goal. */
739 static void target_project_tri_deviation(void *userdata, const float x[3], float r_delta[3])
740 {
741   TargetProjectTriData *data = userdata;
742
743   float w[3] = {x[0], x[1], 1.0f - x[0] - x[1]};
744   interp_v3_v3v3v3(data->co_interp, data->vtri_co[0], data->vtri_co[1], data->vtri_co[2], w);
745   interp_v3_v3v3v3(data->no_interp, data->vtri_no[0], data->vtri_no[1], data->vtri_no[2], w);
746
747   madd_v3_v3v3fl(r_delta, data->co_interp, data->no_interp, x[2]);
748   sub_v3_v3(r_delta, data->point_co);
749 }
750
751 /* Computes the Jacobian matrix of the equation system. */
752 static void target_project_tri_jacobian(void *userdata, const float x[3], float r_jacobian[3][3])
753 {
754   TargetProjectTriData *data = userdata;
755
756   madd_v3_v3v3fl(r_jacobian[0], data->c0_minus_c2, data->n0_minus_n2, x[2]);
757   madd_v3_v3v3fl(r_jacobian[1], data->c1_minus_c2, data->n1_minus_n2, x[2]);
758
759   copy_v3_v3(r_jacobian[2], data->vtri_no[2]);
760   madd_v3_v3fl(r_jacobian[2], data->n0_minus_n2, x[0]);
761   madd_v3_v3fl(r_jacobian[2], data->n1_minus_n2, x[1]);
762 }
763
764 /* Clamp barycentric weights to the triangle. */
765 static void target_project_tri_clamp(float x[3])
766 {
767   if (x[0] < 0.0f) {
768     x[0] = 0.0f;
769   }
770   if (x[1] < 0.0f) {
771     x[1] = 0.0f;
772   }
773   if (x[0] + x[1] > 1.0f) {
774     x[0] = x[0] / (x[0] + x[1]);
775     x[1] = 1.0f - x[0];
776   }
777 }
778
779 /* Correct the Newton's method step to keep the coordinates within the triangle. */
780 static bool target_project_tri_correct(void *UNUSED(userdata),
781                                        const float x[3],
782                                        float step[3],
783                                        float x_next[3])
784 {
785   /* Insignificant correction threshold */
786   const float epsilon = 1e-6f;
787   const float dir_epsilon = 0.05f;
788   bool fixed = false, locked = false;
789
790   /* Weight 0 and 1 boundary check. */
791   for (int i = 0; i < 2; i++) {
792     if (step[i] > x[i]) {
793       if (step[i] > dir_epsilon * fabsf(step[1 - i])) {
794         /* Abort if the solution is clearly outside the domain. */
795         if (x[i] < epsilon) {
796           return false;
797         }
798
799         /* Scale a significant step down to arrive at the boundary. */
800         mul_v3_fl(step, x[i] / step[i]);
801         fixed = true;
802       }
803       else {
804         /* Reset precision errors to stay at the boundary. */
805         step[i] = x[i];
806         fixed = locked = true;
807       }
808     }
809   }
810
811   /* Weight 2 boundary check. */
812   float sum = x[0] + x[1];
813   float sstep = step[0] + step[1];
814
815   if (sum - sstep > 1.0f) {
816     if (sstep < -dir_epsilon * (fabsf(step[0]) + fabsf(step[1]))) {
817       /* Abort if the solution is clearly outside the domain. */
818       if (sum > 1.0f - epsilon) {
819         return false;
820       }
821
822       /* Scale a significant step down to arrive at the boundary. */
823       mul_v3_fl(step, (1.0f - sum) / -sstep);
824       fixed = true;
825     }
826     else {
827       /* Reset precision errors to stay at the boundary. */
828       if (locked) {
829         step[0] = step[1] = 0.0f;
830       }
831       else {
832         step[0] -= 0.5f * sstep;
833         step[1] = -step[0];
834         fixed = true;
835       }
836     }
837   }
838
839   /* Recompute and clamp the new coordinates after step correction. */
840   if (fixed) {
841     sub_v3_v3v3(x_next, x, step);
842     target_project_tri_clamp(x_next);
843   }
844
845   return true;
846 }
847
848 static bool target_project_solve_point_tri(const float *vtri_co[3],
849                                            const float vtri_no[3][3],
850                                            const float point_co[3],
851                                            const float hit_co[3],
852                                            float hit_dist_sq,
853                                            float r_hit_co[3],
854                                            float r_hit_no[3])
855 {
856   float x[3], tmp[3];
857   float dist = sqrtf(hit_dist_sq);
858   float epsilon = dist * 1.0e-5f;
859
860   CLAMP_MIN(epsilon, 1.0e-5f);
861
862   /* Initial solution vector: barycentric weights plus distance along normal. */
863   interp_weights_tri_v3(x, UNPACK3(vtri_co), hit_co);
864
865   interp_v3_v3v3v3(r_hit_no, UNPACK3(vtri_no), x);
866   sub_v3_v3v3(tmp, point_co, hit_co);
867
868   x[2] = (dot_v3v3(tmp, r_hit_no) < 0) ? -dist : dist;
869
870   /* Solve the equations iteratively. */
871   TargetProjectTriData tri_data = {
872       .vtri_co = vtri_co,
873       .vtri_no = vtri_no,
874       .point_co = point_co,
875   };
876
877   sub_v3_v3v3(tri_data.n0_minus_n2, vtri_no[0], vtri_no[2]);
878   sub_v3_v3v3(tri_data.n1_minus_n2, vtri_no[1], vtri_no[2]);
879   sub_v3_v3v3(tri_data.c0_minus_c2, vtri_co[0], vtri_co[2]);
880   sub_v3_v3v3(tri_data.c1_minus_c2, vtri_co[1], vtri_co[2]);
881
882   target_project_tri_clamp(x);
883
884 #ifdef TRACE_TARGET_PROJECT
885   const bool trace = true;
886 #else
887   const bool trace = false;
888 #endif
889
890   bool ok = BLI_newton3d_solve(target_project_tri_deviation,
891                                target_project_tri_jacobian,
892                                target_project_tri_correct,
893                                &tri_data,
894                                epsilon,
895                                20,
896                                trace,
897                                x,
898                                x);
899
900   if (ok) {
901     copy_v3_v3(r_hit_co, tri_data.co_interp);
902     copy_v3_v3(r_hit_no, tri_data.no_interp);
903
904     return true;
905   }
906
907   return false;
908 }
909
910 static bool update_hit(BVHTreeNearest *nearest,
911                        int index,
912                        const float co[3],
913                        const float hit_co[3],
914                        const float hit_no[3])
915 {
916   float dist_sq = len_squared_v3v3(hit_co, co);
917
918   if (dist_sq < nearest->dist_sq) {
919 #ifdef TRACE_TARGET_PROJECT
920     printf(
921         "===> %d (%.3f,%.3f,%.3f) %g < %g\n", index, UNPACK3(hit_co), dist_sq, nearest->dist_sq);
922 #endif
923     nearest->index = index;
924     nearest->dist_sq = dist_sq;
925     copy_v3_v3(nearest->co, hit_co);
926     normalize_v3_v3(nearest->no, hit_no);
927     return true;
928   }
929
930   return false;
931 }
932
933 /* Target projection on a non-manifold boundary edge - treats it like an infinitely thin cylinder. */
934 static void target_project_edge(const ShrinkwrapTreeData *tree,
935                                 int index,
936                                 const float co[3],
937                                 BVHTreeNearest *nearest,
938                                 int eidx)
939 {
940   const BVHTreeFromMesh *data = &tree->treeData;
941   const MEdge *edge = &tree->mesh->medge[eidx];
942   const float *vedge_co[2] = {data->vert[edge->v1].co, data->vert[edge->v2].co};
943
944 #ifdef TRACE_TARGET_PROJECT
945   printf("EDGE %d (%.3f,%.3f,%.3f) (%.3f,%.3f,%.3f)\n",
946          eidx,
947          UNPACK3(vedge_co[0]),
948          UNPACK3(vedge_co[1]));
949 #endif
950
951   /* Retrieve boundary vertex IDs */
952   const int *vert_boundary_id = tree->boundary->vert_boundary_id;
953   int bid1 = vert_boundary_id[edge->v1], bid2 = vert_boundary_id[edge->v2];
954
955   if (bid1 < 0 || bid2 < 0) {
956     return;
957   }
958
959   /* Retrieve boundary vertex normals and align them to direction. */
960   const ShrinkwrapBoundaryVertData *boundary_verts = tree->boundary->boundary_verts;
961   float vedge_dir[2][3], dir[3];
962
963   copy_v3_v3(vedge_dir[0], boundary_verts[bid1].normal_plane);
964   copy_v3_v3(vedge_dir[1], boundary_verts[bid2].normal_plane);
965
966   sub_v3_v3v3(dir, vedge_co[1], vedge_co[0]);
967
968   if (dot_v3v3(boundary_verts[bid1].direction, dir) < 0) {
969     negate_v3(vedge_dir[0]);
970   }
971   if (dot_v3v3(boundary_verts[bid2].direction, dir) < 0) {
972     negate_v3(vedge_dir[1]);
973   }
974
975   /* Solve a quadratic equation: lerp(d0,d1,x) * (co - lerp(v0,v1,x)) = 0 */
976   float d0v0 = dot_v3v3(vedge_dir[0], vedge_co[0]), d0v1 = dot_v3v3(vedge_dir[0], vedge_co[1]);
977   float d1v0 = dot_v3v3(vedge_dir[1], vedge_co[0]), d1v1 = dot_v3v3(vedge_dir[1], vedge_co[1]);
978   float d0co = dot_v3v3(vedge_dir[0], co);
979
980   float a = d0v1 - d0v0 + d1v0 - d1v1;
981   float b = 2 * d0v0 - d0v1 - d0co - d1v0 + dot_v3v3(vedge_dir[1], co);
982   float c = d0co - d0v0;
983   float det = b * b - 4 * a * c;
984
985   if (det >= 0) {
986     const float epsilon = 1e-6f;
987     float sdet = sqrtf(det);
988     float hit_co[3], hit_no[3];
989
990     for (int i = (det > 0 ? 2 : 0); i >= 0; i -= 2) {
991       float x = (-b + ((float)i - 1) * sdet) / (2 * a);
992
993       if (x >= -epsilon && x <= 1.0f + epsilon) {
994         CLAMP(x, 0, 1);
995
996         float vedge_no[2][3];
997         normal_short_to_float_v3(vedge_no[0], data->vert[edge->v1].no);
998         normal_short_to_float_v3(vedge_no[1], data->vert[edge->v2].no);
999
1000         interp_v3_v3v3(hit_co, vedge_co[0], vedge_co[1], x);
1001         interp_v3_v3v3(hit_no, vedge_no[0], vedge_no[1], x);
1002
1003         update_hit(nearest, index, co, hit_co, hit_no);
1004       }
1005     }
1006   }
1007 }
1008
1009 /* Target normal projection BVH callback - based on mesh_looptri_nearest_point. */
1010 static void mesh_looptri_target_project(void *userdata,
1011                                         int index,
1012                                         const float co[3],
1013                                         BVHTreeNearest *nearest)
1014 {
1015   const ShrinkwrapTreeData *tree = (ShrinkwrapTreeData *)userdata;
1016   const BVHTreeFromMesh *data = &tree->treeData;
1017   const MLoopTri *lt = &data->looptri[index];
1018   const MLoop *loop[3] = {
1019       &data->loop[lt->tri[0]], &data->loop[lt->tri[1]], &data->loop[lt->tri[2]]};
1020   const MVert *vtri[3] = {
1021       &data->vert[loop[0]->v], &data->vert[loop[1]->v], &data->vert[loop[2]->v]};
1022   const float *vtri_co[3] = {vtri[0]->co, vtri[1]->co, vtri[2]->co};
1023   float raw_hit_co[3], hit_co[3], hit_no[3], dist_sq, vtri_no[3][3];
1024
1025   /* First find the closest point and bail out if it's worse than the current solution. */
1026   closest_on_tri_to_point_v3(raw_hit_co, co, UNPACK3(vtri_co));
1027   dist_sq = len_squared_v3v3(co, raw_hit_co);
1028
1029 #ifdef TRACE_TARGET_PROJECT
1030   printf("TRIANGLE %d (%.3f,%.3f,%.3f) (%.3f,%.3f,%.3f) (%.3f,%.3f,%.3f) %g %g\n",
1031          index,
1032          UNPACK3(vtri_co[0]),
1033          UNPACK3(vtri_co[1]),
1034          UNPACK3(vtri_co[2]),
1035          dist_sq,
1036          nearest->dist_sq);
1037 #endif
1038
1039   if (dist_sq >= nearest->dist_sq)
1040     return;
1041
1042   /* Decode normals */
1043   normal_short_to_float_v3(vtri_no[0], vtri[0]->no);
1044   normal_short_to_float_v3(vtri_no[1], vtri[1]->no);
1045   normal_short_to_float_v3(vtri_no[2], vtri[2]->no);
1046
1047   /* Solve the equations for the triangle */
1048   if (target_project_solve_point_tri(vtri_co, vtri_no, co, raw_hit_co, dist_sq, hit_co, hit_no)) {
1049     update_hit(nearest, index, co, hit_co, hit_no);
1050   }
1051   /* Boundary edges */
1052   else if (tree->boundary && BLI_BITMAP_TEST(tree->boundary->looptri_has_boundary, index)) {
1053     const BLI_bitmap *is_boundary = tree->boundary->edge_is_boundary;
1054     int edges[3];
1055
1056     BKE_mesh_looptri_get_real_edges(tree->mesh, lt, edges);
1057
1058     for (int i = 0; i < 3; i++) {
1059       if (edges[i] >= 0 && BLI_BITMAP_TEST(is_boundary, edges[i])) {
1060         target_project_edge(tree, index, co, nearest, edges[i]);
1061       }
1062     }
1063   }
1064 }
1065
1066 /*
1067  * Maps the point to the nearest surface, either by simple nearest, or by target normal projection.
1068  */
1069 void BKE_shrinkwrap_find_nearest_surface(struct ShrinkwrapTreeData *tree,
1070                                          BVHTreeNearest *nearest,
1071                                          float co[3],
1072                                          int type)
1073 {
1074   BVHTreeFromMesh *treeData = &tree->treeData;
1075
1076   if (type == MOD_SHRINKWRAP_TARGET_PROJECT) {
1077 #ifdef TRACE_TARGET_PROJECT
1078     printf("====== TARGET PROJECT START ======\n");
1079 #endif
1080
1081     BLI_bvhtree_find_nearest_ex(
1082         tree->bvh, co, nearest, mesh_looptri_target_project, tree, BVH_NEAREST_OPTIMAL_ORDER);
1083
1084 #ifdef TRACE_TARGET_PROJECT
1085     printf("====== TARGET PROJECT END: %d %g ======\n", nearest->index, nearest->dist_sq);
1086 #endif
1087
1088     if (nearest->index < 0) {
1089       /* fallback to simple nearest */
1090       BLI_bvhtree_find_nearest(tree->bvh, co, nearest, treeData->nearest_callback, treeData);
1091     }
1092   }
1093   else {
1094     BLI_bvhtree_find_nearest(tree->bvh, co, nearest, treeData->nearest_callback, treeData);
1095   }
1096 }
1097
1098 /*
1099  * Shrinkwrap moving vertexs to the nearest surface point on the target
1100  *
1101  * it builds a BVHTree from the target mesh and then performs a
1102  * NN matches for each vertex
1103  */
1104 static void shrinkwrap_calc_nearest_surface_point_cb_ex(void *__restrict userdata,
1105                                                         const int i,
1106                                                         const ParallelRangeTLS *__restrict tls)
1107 {
1108   ShrinkwrapCalcCBData *data = userdata;
1109
1110   ShrinkwrapCalcData *calc = data->calc;
1111   BVHTreeNearest *nearest = tls->userdata_chunk;
1112
1113   float *co = calc->vertexCos[i];
1114   float tmp_co[3];
1115   float weight = defvert_array_find_weight_safe(calc->dvert, i, calc->vgroup);
1116
1117   if (calc->invert_vgroup) {
1118     weight = 1.0f - weight;
1119   }
1120
1121   if (weight == 0.0f) {
1122     return;
1123   }
1124
1125   /* Convert the vertex to tree coordinates */
1126   if (calc->vert) {
1127     copy_v3_v3(tmp_co, calc->vert[i].co);
1128   }
1129   else {
1130     copy_v3_v3(tmp_co, co);
1131   }
1132   BLI_space_transform_apply(&calc->local2target, tmp_co);
1133
1134   /* Use local proximity heuristics (to reduce the nearest search)
1135    *
1136    * If we already had an hit before.. we assume this vertex is going to have a close hit to that other vertex
1137    * so we can initiate the "nearest.dist" with the expected value to that last hit.
1138    * This will lead in pruning of the search tree. */
1139   if (nearest->index != -1) {
1140     if (calc->smd->shrinkType == MOD_SHRINKWRAP_TARGET_PROJECT) {
1141       /* Heuristic doesn't work because of additional restrictions. */
1142       nearest->index = -1;
1143       nearest->dist_sq = FLT_MAX;
1144     }
1145     else {
1146       nearest->dist_sq = len_squared_v3v3(tmp_co, nearest->co);
1147     }
1148   }
1149   else
1150     nearest->dist_sq = FLT_MAX;
1151
1152   BKE_shrinkwrap_find_nearest_surface(data->tree, nearest, tmp_co, calc->smd->shrinkType);
1153
1154   /* Found the nearest vertex */
1155   if (nearest->index != -1) {
1156     BKE_shrinkwrap_snap_point_to_surface(data->tree,
1157                                          NULL,
1158                                          calc->smd->shrinkMode,
1159                                          nearest->index,
1160                                          nearest->co,
1161                                          nearest->no,
1162                                          calc->keepDist,
1163                                          tmp_co,
1164                                          tmp_co);
1165
1166     /* Convert the coordinates back to mesh coordinates */
1167     BLI_space_transform_invert(&calc->local2target, tmp_co);
1168     interp_v3_v3v3(co, co, tmp_co, weight); /* linear interpolation */
1169   }
1170 }
1171
1172 /**
1173  * Compute a smooth normal of the target (if applicable) at the hit location.
1174  *
1175  * \param tree: information about the mesh
1176  * \param transform: transform from the hit coordinate space to the object space; may be null
1177  * \param r_no: output in hit coordinate space; may be shared with inputs
1178  */
1179 void BKE_shrinkwrap_compute_smooth_normal(const struct ShrinkwrapTreeData *tree,
1180                                           const struct SpaceTransform *transform,
1181                                           int looptri_idx,
1182                                           const float hit_co[3],
1183                                           const float hit_no[3],
1184                                           float r_no[3])
1185 {
1186   const BVHTreeFromMesh *treeData = &tree->treeData;
1187   const MLoopTri *tri = &treeData->looptri[looptri_idx];
1188
1189   /* Interpolate smooth normals if enabled. */
1190   if ((tree->mesh->mpoly[tri->poly].flag & ME_SMOOTH) != 0) {
1191     const MVert *verts[] = {
1192         &treeData->vert[treeData->loop[tri->tri[0]].v],
1193         &treeData->vert[treeData->loop[tri->tri[1]].v],
1194         &treeData->vert[treeData->loop[tri->tri[2]].v],
1195     };
1196     float w[3], no[3][3], tmp_co[3];
1197
1198     /* Custom and auto smooth split normals. */
1199     if (tree->clnors) {
1200       copy_v3_v3(no[0], tree->clnors[tri->tri[0]]);
1201       copy_v3_v3(no[1], tree->clnors[tri->tri[1]]);
1202       copy_v3_v3(no[2], tree->clnors[tri->tri[2]]);
1203     }
1204     /* Ordinary vertex normals. */
1205     else {
1206       normal_short_to_float_v3(no[0], verts[0]->no);
1207       normal_short_to_float_v3(no[1], verts[1]->no);
1208       normal_short_to_float_v3(no[2], verts[2]->no);
1209     }
1210
1211     /* Barycentric weights from hit point. */
1212     copy_v3_v3(tmp_co, hit_co);
1213
1214     if (transform) {
1215       BLI_space_transform_apply(transform, tmp_co);
1216     }
1217
1218     interp_weights_tri_v3(w, verts[0]->co, verts[1]->co, verts[2]->co, tmp_co);
1219
1220     /* Interpolate using weights. */
1221     interp_v3_v3v3v3(r_no, no[0], no[1], no[2], w);
1222
1223     if (transform) {
1224       BLI_space_transform_invert_normal(transform, r_no);
1225     }
1226     else {
1227       normalize_v3(r_no);
1228     }
1229   }
1230   /* Use the polygon normal if flat. */
1231   else if (tree->pnors != NULL) {
1232     copy_v3_v3(r_no, tree->pnors[tri->poly]);
1233   }
1234   /* Finally fallback to the looptri normal. */
1235   else {
1236     copy_v3_v3(r_no, hit_no);
1237   }
1238 }
1239
1240 /* Helper for MOD_SHRINKWRAP_INSIDE, MOD_SHRINKWRAP_OUTSIDE and MOD_SHRINKWRAP_OUTSIDE_SURFACE. */
1241 static void shrinkwrap_snap_with_side(float r_point_co[3],
1242                                       const float point_co[3],
1243                                       const float hit_co[3],
1244                                       const float hit_no[3],
1245                                       float goal_dist,
1246                                       float forcesign,
1247                                       bool forcesnap)
1248 {
1249   float dist = len_v3v3(point_co, hit_co);
1250
1251   /* If exactly on the surface, push out along normal */
1252   if (dist < FLT_EPSILON) {
1253     if (forcesnap || goal_dist > 0) {
1254       madd_v3_v3v3fl(r_point_co, hit_co, hit_no, goal_dist * forcesign);
1255     }
1256     else {
1257       copy_v3_v3(r_point_co, hit_co);
1258     }
1259   }
1260   /* Move to the correct side if needed */
1261   else {
1262     float delta[3];
1263     sub_v3_v3v3(delta, point_co, hit_co);
1264     float dsign = signf(dot_v3v3(delta, hit_no) * forcesign);
1265
1266     /* If on the wrong side or too close, move to correct */
1267     if (forcesnap || dsign * dist < goal_dist) {
1268       interp_v3_v3v3(r_point_co, point_co, hit_co, (dist - goal_dist * dsign) / dist);
1269     }
1270     else {
1271       copy_v3_v3(r_point_co, point_co);
1272     }
1273   }
1274 }
1275
1276 /**
1277  * Apply the shrink to surface modes to the given original coordinates and nearest point.
1278  *
1279  * \param tree: mesh data for smooth normals
1280  * \param transform: transform from the hit coordinate space to the object space; may be null
1281  * \param r_point_co: may be the same memory location as point_co, hit_co, or hit_no.
1282  */
1283 void BKE_shrinkwrap_snap_point_to_surface(const struct ShrinkwrapTreeData *tree,
1284                                           const struct SpaceTransform *transform,
1285                                           int mode,
1286                                           int hit_idx,
1287                                           const float hit_co[3],
1288                                           const float hit_no[3],
1289                                           float goal_dist,
1290                                           const float point_co[3],
1291                                           float r_point_co[3])
1292 {
1293   float dist, tmp[3];
1294
1295   switch (mode) {
1296     /* Offsets along the line between point_co and hit_co. */
1297     case MOD_SHRINKWRAP_ON_SURFACE:
1298       if (goal_dist != 0 && (dist = len_v3v3(point_co, hit_co)) > FLT_EPSILON) {
1299         interp_v3_v3v3(r_point_co, point_co, hit_co, (dist - goal_dist) / dist);
1300       }
1301       else {
1302         copy_v3_v3(r_point_co, hit_co);
1303       }
1304       break;
1305
1306     case MOD_SHRINKWRAP_INSIDE:
1307       shrinkwrap_snap_with_side(r_point_co, point_co, hit_co, hit_no, goal_dist, -1, false);
1308       break;
1309
1310     case MOD_SHRINKWRAP_OUTSIDE:
1311       shrinkwrap_snap_with_side(r_point_co, point_co, hit_co, hit_no, goal_dist, +1, false);
1312       break;
1313
1314     case MOD_SHRINKWRAP_OUTSIDE_SURFACE:
1315       if (goal_dist != 0) {
1316         shrinkwrap_snap_with_side(r_point_co, point_co, hit_co, hit_no, goal_dist, +1, true);
1317       }
1318       else {
1319         copy_v3_v3(r_point_co, hit_co);
1320       }
1321       break;
1322
1323     /* Offsets along the normal */
1324     case MOD_SHRINKWRAP_ABOVE_SURFACE:
1325       if (goal_dist != 0) {
1326         BKE_shrinkwrap_compute_smooth_normal(tree, transform, hit_idx, hit_co, hit_no, tmp);
1327         madd_v3_v3v3fl(r_point_co, hit_co, tmp, goal_dist);
1328       }
1329       else {
1330         copy_v3_v3(r_point_co, hit_co);
1331       }
1332       break;
1333
1334     default:
1335       printf("Unknown Shrinkwrap surface snap mode: %d\n", mode);
1336       copy_v3_v3(r_point_co, hit_co);
1337   }
1338 }
1339
1340 static void shrinkwrap_calc_nearest_surface_point(ShrinkwrapCalcData *calc)
1341 {
1342   BVHTreeNearest nearest = NULL_BVHTreeNearest;
1343
1344   /* Setup nearest */
1345   nearest.index = -1;
1346   nearest.dist_sq = FLT_MAX;
1347
1348   /* Find the nearest vertex */
1349   ShrinkwrapCalcCBData data = {
1350       .calc = calc,
1351       .tree = calc->tree,
1352   };
1353   ParallelRangeSettings settings;
1354   BLI_parallel_range_settings_defaults(&settings);
1355   settings.use_threading = (calc->numVerts > BKE_MESH_OMP_LIMIT);
1356   settings.userdata_chunk = &nearest;
1357   settings.userdata_chunk_size = sizeof(nearest);
1358   BLI_task_parallel_range(
1359       0, calc->numVerts, &data, shrinkwrap_calc_nearest_surface_point_cb_ex, &settings);
1360 }
1361
1362 /* Main shrinkwrap function */
1363 void shrinkwrapModifier_deform(ShrinkwrapModifierData *smd,
1364                                const ModifierEvalContext *ctx,
1365                                struct Scene *scene,
1366                                Object *ob,
1367                                Mesh *mesh,
1368                                MDeformVert *dvert,
1369                                const int defgrp_index,
1370                                float (*vertexCos)[3],
1371                                int numVerts)
1372 {
1373
1374   DerivedMesh *ss_mesh = NULL;
1375   ShrinkwrapCalcData calc = NULL_ShrinkwrapCalcData;
1376
1377   /* remove loop dependencies on derived meshes (TODO should this be done elsewhere?) */
1378   if (smd->target == ob)
1379     smd->target = NULL;
1380   if (smd->auxTarget == ob)
1381     smd->auxTarget = NULL;
1382
1383   /* Configure Shrinkwrap calc data */
1384   calc.smd = smd;
1385   calc.ob = ob;
1386   calc.numVerts = numVerts;
1387   calc.vertexCos = vertexCos;
1388   calc.dvert = dvert;
1389   calc.vgroup = defgrp_index;
1390   calc.invert_vgroup = (smd->shrinkOpts & MOD_SHRINKWRAP_INVERT_VGROUP) != 0;
1391
1392   if (smd->target != NULL) {
1393     Object *ob_target = DEG_get_evaluated_object(ctx->depsgraph, smd->target);
1394     calc.target = BKE_modifier_get_evaluated_mesh_from_evaluated_object(ob_target, false);
1395
1396     /* TODO there might be several "bugs" on non-uniform scales matrixs
1397      * because it will no longer be nearest surface, not sphere projection
1398      * because space has been deformed */
1399     BLI_SPACE_TRANSFORM_SETUP(&calc.local2target, ob, ob_target);
1400
1401     /* TODO: smd->keepDist is in global units.. must change to local */
1402     calc.keepDist = smd->keepDist;
1403   }
1404   calc.aux_target = DEG_get_evaluated_object(ctx->depsgraph, smd->auxTarget);
1405
1406   if (mesh != NULL && smd->shrinkType == MOD_SHRINKWRAP_PROJECT) {
1407     /* Setup arrays to get vertexs positions, normals and deform weights */
1408     calc.vert = mesh->mvert;
1409
1410     /* Using vertexs positions/normals as if a subsurface was applied */
1411     if (smd->subsurfLevels) {
1412       SubsurfModifierData ssmd = {{NULL}};
1413       ssmd.subdivType = ME_CC_SUBSURF;  /* catmull clark */
1414       ssmd.levels = smd->subsurfLevels; /* levels */
1415
1416       /* TODO to be moved to Mesh once we are done with changes in subsurf code. */
1417       DerivedMesh *dm = CDDM_from_mesh(mesh);
1418
1419       ss_mesh = subsurf_make_derived_from_derived(
1420           dm, &ssmd, scene, NULL, (ob->mode & OB_MODE_EDIT) ? SUBSURF_IN_EDIT_MODE : 0);
1421
1422       if (ss_mesh) {
1423         calc.vert = ss_mesh->getVertDataArray(ss_mesh, CD_MVERT);
1424         if (calc.vert) {
1425           /* TRICKY: this code assumes subsurface will have the transformed original vertices
1426            * in their original order at the end of the vert array. */
1427           calc.vert = calc.vert + ss_mesh->getNumVerts(ss_mesh) - dm->getNumVerts(dm);
1428         }
1429       }
1430
1431       /* Just to make sure we are not leaving any memory behind */
1432       BLI_assert(ssmd.emCache == NULL);
1433       BLI_assert(ssmd.mCache == NULL);
1434
1435       dm->release(dm);
1436     }
1437   }
1438
1439   /* Projecting target defined - lets work! */
1440   ShrinkwrapTreeData tree;
1441
1442   if (BKE_shrinkwrap_init_tree(&tree, calc.target, smd->shrinkType, smd->shrinkMode, false)) {
1443     calc.tree = &tree;
1444
1445     switch (smd->shrinkType) {
1446       case MOD_SHRINKWRAP_NEAREST_SURFACE:
1447       case MOD_SHRINKWRAP_TARGET_PROJECT:
1448         TIMEIT_BENCH(shrinkwrap_calc_nearest_surface_point(&calc), deform_surface);
1449         break;
1450
1451       case MOD_SHRINKWRAP_PROJECT:
1452         TIMEIT_BENCH(shrinkwrap_calc_normal_projection(&calc), deform_project);
1453         break;
1454
1455       case MOD_SHRINKWRAP_NEAREST_VERTEX:
1456         TIMEIT_BENCH(shrinkwrap_calc_nearest_vertex(&calc), deform_vertex);
1457         break;
1458     }
1459
1460     BKE_shrinkwrap_free_tree(&tree);
1461   }
1462
1463   /* free memory */
1464   if (ss_mesh)
1465     ss_mesh->release(ss_mesh);
1466 }