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