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