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