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