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         ShrinkwrapTreeData *aux_tree = NULL;
628         ShrinkwrapTreeData aux_tree_stack;
629         SpaceTransform local2aux;
630
631         /* If the user doesn't allows to project in any direction of projection axis
632          * then there's nothing todo. */
633         if ((calc->smd->shrinkOpts & (MOD_SHRINKWRAP_PROJECT_ALLOW_POS_DIR | MOD_SHRINKWRAP_PROJECT_ALLOW_NEG_DIR)) == 0)
634                 return;
635
636         /* Prepare data to retrieve the direction in which we should project each vertex */
637         if (calc->smd->projAxis == MOD_SHRINKWRAP_PROJECT_OVER_NORMAL) {
638                 if (calc->vert == NULL) return;
639         }
640         else {
641                 /* The code supports any axis that is a combination of X,Y,Z
642                  * although currently UI only allows to set the 3 different axis */
643                 if (calc->smd->projAxis & MOD_SHRINKWRAP_PROJECT_OVER_X_AXIS) proj_axis[0] = 1.0f;
644                 if (calc->smd->projAxis & MOD_SHRINKWRAP_PROJECT_OVER_Y_AXIS) proj_axis[1] = 1.0f;
645                 if (calc->smd->projAxis & MOD_SHRINKWRAP_PROJECT_OVER_Z_AXIS) proj_axis[2] = 1.0f;
646
647                 normalize_v3(proj_axis);
648
649                 /* Invalid projection direction */
650                 if (len_squared_v3(proj_axis) < FLT_EPSILON) {
651                         return;
652                 }
653         }
654
655         if (calc->aux_target) {
656                 auxMesh = BKE_modifier_get_evaluated_mesh_from_evaluated_object(calc->aux_target, false);
657                 if (!auxMesh)
658                         return;
659                 BLI_SPACE_TRANSFORM_SETUP(&local2aux, calc->ob, calc->aux_target);
660         }
661
662         if (BKE_shrinkwrap_init_tree(&aux_tree_stack, auxMesh, calc->smd->shrinkType, calc->smd->shrinkMode, false)) {
663                 aux_tree = &aux_tree_stack;
664         }
665
666         /* After successfully build the trees, start projection vertices. */
667         ShrinkwrapCalcCBData data = {
668                 .calc = calc, .tree = calc->tree, .aux_tree = aux_tree,
669                 .proj_axis = proj_axis, .local2aux = &local2aux,
670         };
671         ParallelRangeSettings settings;
672         BLI_parallel_range_settings_defaults(&settings);
673         settings.use_threading = (calc->numVerts > BKE_MESH_OMP_LIMIT);
674         settings.userdata_chunk = &hit;
675         settings.userdata_chunk_size = sizeof(hit);
676         BLI_task_parallel_range(0, calc->numVerts,
677                                 &data,
678                                 shrinkwrap_calc_normal_projection_cb_ex,
679                                 &settings);
680
681         /* free data structures */
682         if (aux_tree) {
683                 BKE_shrinkwrap_free_tree(aux_tree);
684         }
685 }
686
687 /*
688  * Shrinkwrap Target Surface Project mode
689  *
690  * It uses Newton's method to find a surface location with its
691  * smooth normal pointing at the original point.
692  *
693  * The equation system on barycentric weights and normal multiplier:
694  *
695  *   (w0*V0 + w1*V1 + w2*V2) + l * (w0*N0 + w1*N1 + w2*N2) - CO = 0
696  *   w0 + w1 + w2 = 1
697  *
698  * The actual solution vector is [ w0, w1, l ], with w2 eliminated.
699  */
700
701 //#define TRACE_TARGET_PROJECT
702
703 typedef struct TargetProjectTriData {
704         const float **vtri_co;
705         const float (*vtri_no)[3];
706         const float *point_co;
707
708         float n0_minus_n2[3], n1_minus_n2[3];
709         float c0_minus_c2[3], c1_minus_c2[3];
710
711         /* Current interpolated position and normal. */
712         float co_interp[3], no_interp[3];
713 } TargetProjectTriData;
714
715 /* Computes the deviation of the equation system from goal. */
716 static void target_project_tri_deviation(void *userdata, const float x[3], float r_delta[3])
717 {
718         TargetProjectTriData *data = userdata;
719
720         float w[3] = { x[0], x[1], 1.0f - x[0] - x[1] };
721         interp_v3_v3v3v3(data->co_interp, data->vtri_co[0], data->vtri_co[1], data->vtri_co[2], w);
722         interp_v3_v3v3v3(data->no_interp, data->vtri_no[0], data->vtri_no[1], data->vtri_no[2], w);
723
724         madd_v3_v3v3fl(r_delta, data->co_interp, data->no_interp, x[2]);
725         sub_v3_v3(r_delta, data->point_co);
726 }
727
728 /* Computes the Jacobian matrix of the equation system. */
729 static void target_project_tri_jacobian(void *userdata, const float x[3], float r_jacobian[3][3])
730 {
731         TargetProjectTriData *data = userdata;
732
733         madd_v3_v3v3fl(r_jacobian[0], data->c0_minus_c2, data->n0_minus_n2, x[2]);
734         madd_v3_v3v3fl(r_jacobian[1], data->c1_minus_c2, data->n1_minus_n2, x[2]);
735
736         copy_v3_v3(r_jacobian[2], data->vtri_no[2]);
737         madd_v3_v3fl(r_jacobian[2], data->n0_minus_n2, x[0]);
738         madd_v3_v3fl(r_jacobian[2], data->n1_minus_n2, x[1]);
739 }
740
741 /* Clamp barycentric weights to the triangle. */
742 static void target_project_tri_clamp(float x[3])
743 {
744         if (x[0] < 0.0f) {
745                 x[0] = 0.0f;
746         }
747         if (x[1] < 0.0f) {
748                 x[1] = 0.0f;
749         }
750         if (x[0] + x[1] > 1.0f) {
751                 x[0] = x[0] / (x[0] + x[1]);
752                 x[1] = 1.0f - x[0];
753         }
754 }
755
756 /* Correct the Newton's method step to keep the coordinates within the triangle. */
757 static bool target_project_tri_correct(void *UNUSED(userdata), const float x[3], float step[3], float x_next[3])
758 {
759         /* Insignificant correction threshold */
760         const float epsilon = 1e-6f;
761         const float dir_epsilon = 0.05f;
762         bool fixed = false, locked = false;
763
764         /* Weight 0 and 1 boundary check. */
765         for (int i = 0; i < 2; i++) {
766                 if (step[i] > x[i]) {
767                         if (step[i] > dir_epsilon * fabsf(step[1 - i])) {
768                                 /* Abort if the solution is clearly outside the domain. */
769                                 if (x[i] < epsilon) {
770                                         return false;
771                                 }
772
773                                 /* Scale a significant step down to arrive at the boundary. */
774                                 mul_v3_fl(step, x[i] / step[i]);
775                                 fixed = true;
776                         }
777                         else {
778                                 /* Reset precision errors to stay at the boundary. */
779                                 step[i] = x[i];
780                                 fixed = locked = true;
781                         }
782                 }
783         }
784
785         /* Weight 2 boundary check. */
786         float sum = x[0] + x[1];
787         float sstep = step[0] + step[1];
788
789         if (sum - sstep > 1.0f) {
790                 if (sstep < -dir_epsilon * (fabsf(step[0]) + fabsf(step[1]))) {
791                         /* Abort if the solution is clearly outside the domain. */
792                         if (sum > 1.0f - epsilon) {
793                                 return false;
794                         }
795
796                         /* Scale a significant step down to arrive at the boundary. */
797                         mul_v3_fl(step, (1.0f - sum) / -sstep);
798                         fixed = true;
799                 }
800                 else {
801                         /* Reset precision errors to stay at the boundary. */
802                         if (locked) {
803                                 step[0] = step[1] = 0.0f;
804                         }
805                         else {
806                                 step[0] -= 0.5f * sstep;
807                                 step[1] = -step[0];
808                                 fixed = true;
809                         }
810                 }
811         }
812
813         /* Recompute and clamp the new coordinates after step correction. */
814         if (fixed) {
815                 sub_v3_v3v3(x_next, x, step);
816                 target_project_tri_clamp(x_next);
817         }
818
819         return true;
820 }
821
822 static bool target_project_solve_point_tri(
823         const float *vtri_co[3], const float vtri_no[3][3], const float point_co[3],
824         const float hit_co[3], float hit_dist_sq,
825         float r_hit_co[3], float r_hit_no[3])
826 {
827         float x[3], tmp[3];
828         float dist = sqrtf(hit_dist_sq);
829         float epsilon = dist * 1.0e-5f;
830
831         CLAMP_MIN(epsilon, 1.0e-5f);
832
833         /* Initial solution vector: barycentric weights plus distance along normal. */
834         interp_weights_tri_v3(x, UNPACK3(vtri_co), hit_co);
835
836         interp_v3_v3v3v3(r_hit_no, UNPACK3(vtri_no), x);
837         sub_v3_v3v3(tmp, point_co, hit_co);
838
839         x[2] = (dot_v3v3(tmp, r_hit_no) < 0) ? -dist : dist;
840
841         /* Solve the equations iteratively. */
842         TargetProjectTriData tri_data = { .vtri_co = vtri_co, .vtri_no = vtri_no, .point_co = point_co, };
843
844         sub_v3_v3v3(tri_data.n0_minus_n2, vtri_no[0], vtri_no[2]);
845         sub_v3_v3v3(tri_data.n1_minus_n2, vtri_no[1], vtri_no[2]);
846         sub_v3_v3v3(tri_data.c0_minus_c2, vtri_co[0], vtri_co[2]);
847         sub_v3_v3v3(tri_data.c1_minus_c2, vtri_co[1], vtri_co[2]);
848
849         target_project_tri_clamp(x);
850
851 #ifdef TRACE_TARGET_PROJECT
852         const bool trace = true;
853 #else
854         const bool trace = false;
855 #endif
856
857         bool ok = BLI_newton3d_solve(target_project_tri_deviation, target_project_tri_jacobian, target_project_tri_correct, &tri_data, epsilon, 20, trace, x, x);
858
859         if (ok) {
860                 copy_v3_v3(r_hit_co, tri_data.co_interp);
861                 copy_v3_v3(r_hit_no, tri_data.no_interp);
862
863                 return true;
864         }
865
866         return false;
867 }
868
869 static bool update_hit(BVHTreeNearest *nearest, int index, const float co[3], const float hit_co[3], const float hit_no[3])
870 {
871         float dist_sq = len_squared_v3v3(hit_co, co);
872
873         if (dist_sq < nearest->dist_sq) {
874 #ifdef TRACE_TARGET_PROJECT
875                 printf("===> %d (%.3f,%.3f,%.3f) %g < %g\n", index, UNPACK3(hit_co), dist_sq, nearest->dist_sq);
876 #endif
877                 nearest->index = index;
878                 nearest->dist_sq = dist_sq;
879                 copy_v3_v3(nearest->co, hit_co);
880                 normalize_v3_v3(nearest->no, hit_no);
881                 return true;
882         }
883
884         return false;
885 }
886
887 /* Target projection on a non-manifold boundary edge - treats it like an infinitely thin cylinder. */
888 static void target_project_edge(const ShrinkwrapTreeData *tree, int index, const float co[3], BVHTreeNearest *nearest, int eidx)
889 {
890         const BVHTreeFromMesh *data = &tree->treeData;
891         const MEdge *edge = &tree->mesh->medge[eidx];
892         const float *vedge_co[2] = { data->vert[edge->v1].co, data->vert[edge->v2].co };
893
894 #ifdef TRACE_TARGET_PROJECT
895         printf("EDGE %d (%.3f,%.3f,%.3f) (%.3f,%.3f,%.3f)\n", eidx, UNPACK3(vedge_co[0]), UNPACK3(vedge_co[1]));
896 #endif
897
898         /* Retrieve boundary vertex IDs */
899         const int *vert_boundary_id = tree->boundary->vert_boundary_id;
900         int bid1 = vert_boundary_id[edge->v1], bid2 = vert_boundary_id[edge->v2];
901
902         if (bid1 < 0 || bid2 < 0) {
903                 return;
904         }
905
906         /* Retrieve boundary vertex normals and align them to direction. */
907         const ShrinkwrapBoundaryVertData *boundary_verts = tree->boundary->boundary_verts;
908         float vedge_dir[2][3], dir[3];
909
910         copy_v3_v3(vedge_dir[0], boundary_verts[bid1].normal_plane);
911         copy_v3_v3(vedge_dir[1], boundary_verts[bid2].normal_plane);
912
913         sub_v3_v3v3(dir, vedge_co[1], vedge_co[0]);
914
915         if (dot_v3v3(boundary_verts[bid1].direction, dir) < 0) {
916                 negate_v3(vedge_dir[0]);
917         }
918         if (dot_v3v3(boundary_verts[bid2].direction, dir) < 0) {
919                 negate_v3(vedge_dir[1]);
920         }
921
922         /* Solve a quadratic equation: lerp(d0,d1,x) * (co - lerp(v0,v1,x)) = 0 */
923         float d0v0 = dot_v3v3(vedge_dir[0], vedge_co[0]), d0v1 = dot_v3v3(vedge_dir[0], vedge_co[1]);
924         float d1v0 = dot_v3v3(vedge_dir[1], vedge_co[0]), d1v1 = dot_v3v3(vedge_dir[1], vedge_co[1]);
925         float d0co = dot_v3v3(vedge_dir[0], co);
926
927         float a = d0v1 - d0v0 + d1v0 - d1v1;
928         float b = 2 * d0v0 - d0v1 - d0co - d1v0 + dot_v3v3(vedge_dir[1], co);
929         float c = d0co - d0v0;
930         float det = b * b - 4 * a * c;
931
932         if (det >= 0) {
933                 const float epsilon = 1e-6f;
934                 float sdet = sqrtf(det);
935                 float hit_co[3], hit_no[3];
936
937                 for (int i = (det > 0 ? 2 : 0); i >= 0; i -= 2) {
938                         float x = (- b + ((float)i - 1) * sdet) / (2 * a);
939
940                         if (x >= -epsilon && x <= 1.0f + epsilon) {
941                                 CLAMP(x, 0, 1);
942
943                                 float vedge_no[2][3];
944                                 normal_short_to_float_v3(vedge_no[0], data->vert[edge->v1].no);
945                                 normal_short_to_float_v3(vedge_no[1], data->vert[edge->v2].no);
946
947                                 interp_v3_v3v3(hit_co, vedge_co[0], vedge_co[1], x);
948                                 interp_v3_v3v3(hit_no, vedge_no[0], vedge_no[1], x);
949
950                                 update_hit(nearest, index, co, hit_co, hit_no);
951                         }
952                 }
953         }
954 }
955
956 /* Target normal projection BVH callback - based on mesh_looptri_nearest_point. */
957 static void mesh_looptri_target_project(void *userdata, int index, const float co[3], BVHTreeNearest *nearest)
958 {
959         const ShrinkwrapTreeData *tree = (ShrinkwrapTreeData *)userdata;
960         const BVHTreeFromMesh *data = &tree->treeData;
961         const MLoopTri *lt = &data->looptri[index];
962         const MLoop *loop[3] = { &data->loop[lt->tri[0]], &data->loop[lt->tri[1]], &data->loop[lt->tri[2]] };
963         const MVert *vtri[3] = { &data->vert[loop[0]->v], &data->vert[loop[1]->v], &data->vert[loop[2]->v] };
964         const float *vtri_co[3] = { vtri[0]->co, vtri[1]->co, vtri[2]->co };
965         float raw_hit_co[3], hit_co[3], hit_no[3], dist_sq, vtri_no[3][3];
966
967         /* First find the closest point and bail out if it's worse than the current solution. */
968         closest_on_tri_to_point_v3(raw_hit_co, co, UNPACK3(vtri_co));
969         dist_sq = len_squared_v3v3(co, raw_hit_co);
970
971 #ifdef TRACE_TARGET_PROJECT
972         printf("TRIANGLE %d (%.3f,%.3f,%.3f) (%.3f,%.3f,%.3f) (%.3f,%.3f,%.3f) %g %g\n",
973                index, UNPACK3(vtri_co[0]), UNPACK3(vtri_co[1]), UNPACK3(vtri_co[2]), dist_sq, nearest->dist_sq);
974 #endif
975
976         if (dist_sq >= nearest->dist_sq)
977                 return;
978
979         /* Decode normals */
980         normal_short_to_float_v3(vtri_no[0], vtri[0]->no);
981         normal_short_to_float_v3(vtri_no[1], vtri[1]->no);
982         normal_short_to_float_v3(vtri_no[2], vtri[2]->no);
983
984         /* Solve the equations for the triangle */
985         if (target_project_solve_point_tri(vtri_co, vtri_no, co, raw_hit_co, dist_sq, hit_co, hit_no)) {
986                 update_hit(nearest, index, co, hit_co, hit_no);
987         }
988         /* Boundary edges */
989         else if (tree->boundary && BLI_BITMAP_TEST(tree->boundary->looptri_has_boundary, index)) {
990                 const BLI_bitmap *is_boundary = tree->boundary->edge_is_boundary;
991                 int edges[3];
992
993                 BKE_mesh_looptri_get_real_edges(tree->mesh, lt, edges);
994
995                 for (int i = 0; i < 3; i++) {
996                         if (edges[i] >= 0 && BLI_BITMAP_TEST(is_boundary, edges[i])) {
997                                 target_project_edge(tree, index, co, nearest, edges[i]);
998                         }
999                 }
1000         }
1001 }
1002
1003 /*
1004  * Maps the point to the nearest surface, either by simple nearest, or by target normal projection.
1005  */
1006 void BKE_shrinkwrap_find_nearest_surface(struct ShrinkwrapTreeData *tree, BVHTreeNearest *nearest, float co[3], int type)
1007 {
1008         BVHTreeFromMesh *treeData = &tree->treeData;
1009
1010         if (type == MOD_SHRINKWRAP_TARGET_PROJECT) {
1011 #ifdef TRACE_TARGET_PROJECT
1012                 printf("====== TARGET PROJECT START ======\n");
1013 #endif
1014
1015                 BLI_bvhtree_find_nearest_ex(tree->bvh, co, nearest, mesh_looptri_target_project, tree, BVH_NEAREST_OPTIMAL_ORDER);
1016
1017 #ifdef TRACE_TARGET_PROJECT
1018                 printf("====== TARGET PROJECT END: %d %g ======\n", nearest->index, nearest->dist_sq);
1019 #endif
1020
1021                 if (nearest->index < 0) {
1022                         /* fallback to simple nearest */
1023                         BLI_bvhtree_find_nearest(tree->bvh, co, nearest, treeData->nearest_callback, treeData);
1024                 }
1025         }
1026         else {
1027                 BLI_bvhtree_find_nearest(tree->bvh, co, nearest, treeData->nearest_callback, treeData);
1028         }
1029 }
1030
1031 /*
1032  * Shrinkwrap moving vertexs to the nearest surface point on the target
1033  *
1034  * it builds a BVHTree from the target mesh and then performs a
1035  * NN matches for each vertex
1036  */
1037 static void shrinkwrap_calc_nearest_surface_point_cb_ex(
1038         void *__restrict userdata,
1039         const int i,
1040         const ParallelRangeTLS *__restrict tls)
1041 {
1042         ShrinkwrapCalcCBData *data = userdata;
1043
1044         ShrinkwrapCalcData *calc = data->calc;
1045         BVHTreeNearest *nearest = tls->userdata_chunk;
1046
1047         float *co = calc->vertexCos[i];
1048         float tmp_co[3];
1049         float weight = defvert_array_find_weight_safe(calc->dvert, i, calc->vgroup);
1050
1051         if (calc->invert_vgroup) {
1052                 weight = 1.0f - weight;
1053         }
1054
1055         if (weight == 0.0f) {
1056                 return;
1057         }
1058
1059         /* Convert the vertex to tree coordinates */
1060         if (calc->vert) {
1061                 copy_v3_v3(tmp_co, calc->vert[i].co);
1062         }
1063         else {
1064                 copy_v3_v3(tmp_co, co);
1065         }
1066         BLI_space_transform_apply(&calc->local2target, tmp_co);
1067
1068         /* Use local proximity heuristics (to reduce the nearest search)
1069          *
1070          * If we already had an hit before.. we assume this vertex is going to have a close hit to that other vertex
1071          * so we can initiate the "nearest.dist" with the expected value to that last hit.
1072          * This will lead in pruning of the search tree. */
1073         if (nearest->index != -1) {
1074                 if (calc->smd->shrinkType == MOD_SHRINKWRAP_TARGET_PROJECT) {
1075                         /* Heuristic doesn't work because of additional restrictions. */
1076                         nearest->index = -1;
1077                         nearest->dist_sq = FLT_MAX;
1078                 }
1079                 else {
1080                         nearest->dist_sq = len_squared_v3v3(tmp_co, nearest->co);
1081                 }
1082         }
1083         else
1084                 nearest->dist_sq = FLT_MAX;
1085
1086         BKE_shrinkwrap_find_nearest_surface(data->tree, nearest, tmp_co, calc->smd->shrinkType);
1087
1088         /* Found the nearest vertex */
1089         if (nearest->index != -1) {
1090                 BKE_shrinkwrap_snap_point_to_surface(
1091                         data->tree, NULL, calc->smd->shrinkMode,
1092                         nearest->index, nearest->co, nearest->no, calc->keepDist, tmp_co, tmp_co);
1093
1094                 /* Convert the coordinates back to mesh coordinates */
1095                 BLI_space_transform_invert(&calc->local2target, tmp_co);
1096                 interp_v3_v3v3(co, co, tmp_co, weight);  /* linear interpolation */
1097         }
1098 }
1099
1100 /**
1101  * Compute a smooth normal of the target (if applicable) at the hit location.
1102  *
1103  * \param tree: information about the mesh
1104  * \param transform: transform from the hit coordinate space to the object space; may be null
1105  * \param r_no: output in hit coordinate space; may be shared with inputs
1106  */
1107 void BKE_shrinkwrap_compute_smooth_normal(
1108         const struct ShrinkwrapTreeData *tree, const struct SpaceTransform *transform,
1109         int looptri_idx, const float hit_co[3], const float hit_no[3], float r_no[3])
1110 {
1111         const BVHTreeFromMesh *treeData = &tree->treeData;
1112         const MLoopTri *tri = &treeData->looptri[looptri_idx];
1113
1114         /* Interpolate smooth normals if enabled. */
1115         if ((tree->mesh->mpoly[tri->poly].flag & ME_SMOOTH) != 0) {
1116                 const MVert *verts[] = {
1117                         &treeData->vert[treeData->loop[tri->tri[0]].v],
1118                         &treeData->vert[treeData->loop[tri->tri[1]].v],
1119                         &treeData->vert[treeData->loop[tri->tri[2]].v],
1120                 };
1121                 float w[3], no[3][3], tmp_co[3];
1122
1123                 /* Custom and auto smooth split normals. */
1124                 if (tree->clnors) {
1125                         copy_v3_v3(no[0], tree->clnors[tri->tri[0]]);
1126                         copy_v3_v3(no[1], tree->clnors[tri->tri[1]]);
1127                         copy_v3_v3(no[2], tree->clnors[tri->tri[2]]);
1128                 }
1129                 /* Ordinary vertex normals. */
1130                 else {
1131                         normal_short_to_float_v3(no[0], verts[0]->no);
1132                         normal_short_to_float_v3(no[1], verts[1]->no);
1133                         normal_short_to_float_v3(no[2], verts[2]->no);
1134                 }
1135
1136                 /* Barycentric weights from hit point. */
1137                 copy_v3_v3(tmp_co, hit_co);
1138
1139                 if (transform) {
1140                         BLI_space_transform_apply(transform, tmp_co);
1141                 }
1142
1143                 interp_weights_tri_v3(w, verts[0]->co, verts[1]->co, verts[2]->co, tmp_co);
1144
1145                 /* Interpolate using weights. */
1146                 interp_v3_v3v3v3(r_no, no[0], no[1], no[2], w);
1147
1148                 if (transform) {
1149                         BLI_space_transform_invert_normal(transform, r_no);
1150                 }
1151                 else {
1152                         normalize_v3(r_no);
1153                 }
1154         }
1155         /* Use the polygon normal if flat. */
1156         else if (tree->pnors != NULL) {
1157                 copy_v3_v3(r_no, tree->pnors[tri->poly]);
1158         }
1159         /* Finally fallback to the looptri normal. */
1160         else {
1161                 copy_v3_v3(r_no, hit_no);
1162         }
1163 }
1164
1165 /* Helper for MOD_SHRINKWRAP_INSIDE, MOD_SHRINKWRAP_OUTSIDE and MOD_SHRINKWRAP_OUTSIDE_SURFACE. */
1166 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)
1167 {
1168         float dist = len_v3v3(point_co, hit_co);
1169
1170         /* If exactly on the surface, push out along normal */
1171         if (dist < FLT_EPSILON) {
1172                 if (forcesnap || goal_dist > 0) {
1173                         madd_v3_v3v3fl(r_point_co, hit_co, hit_no, goal_dist * forcesign);
1174                 }
1175                 else {
1176                         copy_v3_v3(r_point_co, hit_co);
1177                 }
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) * forcesign);
1184
1185                 /* If on the wrong side or too close, move to correct */
1186                 if (forcesnap || dsign * dist < goal_dist) {
1187                         interp_v3_v3v3(r_point_co, point_co, hit_co, (dist - goal_dist * dsign) / 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(
1277         ShrinkwrapModifierData *smd, const ModifierEvalContext *ctx,
1278         struct Scene *scene, Object *ob, Mesh *mesh,
1279         MDeformVert *dvert, const int defgrp_index, float (*vertexCos)[3], int numVerts)
1280 {
1281
1282         DerivedMesh *ss_mesh    = NULL;
1283         ShrinkwrapCalcData calc = NULL_ShrinkwrapCalcData;
1284
1285         /* remove loop dependencies on derived meshes (TODO should this be done elsewhere?) */
1286         if (smd->target == ob) smd->target = NULL;
1287         if (smd->auxTarget == ob) smd->auxTarget = NULL;
1288
1289
1290         /* Configure Shrinkwrap calc data */
1291         calc.smd = smd;
1292         calc.ob = ob;
1293         calc.numVerts = numVerts;
1294         calc.vertexCos = vertexCos;
1295         calc.dvert = dvert;
1296         calc.vgroup = defgrp_index;
1297         calc.invert_vgroup = (smd->shrinkOpts & MOD_SHRINKWRAP_INVERT_VGROUP) != 0;
1298
1299         if (smd->target != NULL) {
1300                 Object *ob_target = DEG_get_evaluated_object(ctx->depsgraph, smd->target);
1301                 calc.target = BKE_modifier_get_evaluated_mesh_from_evaluated_object(ob_target, false);
1302
1303                 /* TODO there might be several "bugs" on non-uniform scales matrixs
1304                  * because it will no longer be nearest surface, not sphere projection
1305                  * because space has been deformed */
1306                 BLI_SPACE_TRANSFORM_SETUP(&calc.local2target, ob, ob_target);
1307
1308                 /* TODO: smd->keepDist is in global units.. must change to local */
1309                 calc.keepDist = smd->keepDist;
1310         }
1311         calc.aux_target = DEG_get_evaluated_object(ctx->depsgraph, smd->auxTarget);
1312
1313         if (mesh != NULL && smd->shrinkType == MOD_SHRINKWRAP_PROJECT) {
1314                 /* Setup arrays to get vertexs positions, normals and deform weights */
1315                 calc.vert = mesh->mvert;
1316
1317                 /* Using vertexs positions/normals as if a subsurface was applied */
1318                 if (smd->subsurfLevels) {
1319                         SubsurfModifierData ssmd = {{NULL}};
1320                         ssmd.subdivType = ME_CC_SUBSURF;        /* catmull clark */
1321                         ssmd.levels     = smd->subsurfLevels;   /* levels */
1322
1323                         /* TODO to be moved to Mesh once we are done with changes in subsurf code. */
1324                         DerivedMesh *dm = CDDM_from_mesh(mesh);
1325
1326                         ss_mesh = subsurf_make_derived_from_derived(dm, &ssmd, scene, NULL, (ob->mode & OB_MODE_EDIT) ? SUBSURF_IN_EDIT_MODE : 0);
1327
1328                         if (ss_mesh) {
1329                                 calc.vert = ss_mesh->getVertDataArray(ss_mesh, CD_MVERT);
1330                                 if (calc.vert) {
1331                                         /* TRICKY: this code assumes subsurface will have the transformed original vertices
1332                                          * in their original order at the end of the vert array. */
1333                                         calc.vert = calc.vert + ss_mesh->getNumVerts(ss_mesh) - dm->getNumVerts(dm);
1334                                 }
1335                         }
1336
1337                         /* Just to make sure we are not leaving any memory behind */
1338                         BLI_assert(ssmd.emCache == NULL);
1339                         BLI_assert(ssmd.mCache == NULL);
1340
1341                         dm->release(dm);
1342                 }
1343         }
1344
1345         /* Projecting target defined - lets work! */
1346         ShrinkwrapTreeData tree;
1347
1348         if (BKE_shrinkwrap_init_tree(&tree, calc.target, smd->shrinkType, smd->shrinkMode, false)) {
1349                 calc.tree = &tree;
1350
1351                 switch (smd->shrinkType) {
1352                         case MOD_SHRINKWRAP_NEAREST_SURFACE:
1353                         case MOD_SHRINKWRAP_TARGET_PROJECT:
1354                                 TIMEIT_BENCH(shrinkwrap_calc_nearest_surface_point(&calc), deform_surface);
1355                                 break;
1356
1357                         case MOD_SHRINKWRAP_PROJECT:
1358                                 TIMEIT_BENCH(shrinkwrap_calc_normal_projection(&calc), deform_project);
1359                                 break;
1360
1361                         case MOD_SHRINKWRAP_NEAREST_VERTEX:
1362                                 TIMEIT_BENCH(shrinkwrap_calc_nearest_vertex(&calc), deform_vertex);
1363                                 break;
1364                 }
1365
1366                 BKE_shrinkwrap_free_tree(&tree);
1367         }
1368
1369         /* free memory */
1370         if (ss_mesh)
1371                 ss_mesh->release(ss_mesh);
1372 }