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