Threaded object update and EvaluationContext
[blender.git] / source / blender / blenkernel / intern / shrinkwrap.c
1 /*
2  * ***** BEGIN GPL LICENSE BLOCK *****
3  *
4  * This program is free software; you can redistribute it and/or
5  * modify it under the terms of the GNU General Public License
6  * as published by the Free Software Foundation; either version 2
7  * of the License, or (at your option) any later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software Foundation,
16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  *
18  * The Original Code is Copyright (C) Blender Foundation.
19  * All rights reserved.
20  *
21  * The Original Code is: all of this file.
22  *
23  * Contributor(s): Andr Pinto
24  *
25  * ***** END GPL LICENSE BLOCK *****
26  */
27
28 /** \file blender/blenkernel/intern/shrinkwrap.c
29  *  \ingroup bke
30  */
31
32 #include <string.h>
33 #include <float.h>
34 #include <math.h>
35 #include <memory.h>
36 #include <stdio.h>
37 #include <time.h>
38 #include <assert.h>
39
40 #include "DNA_object_types.h"
41 #include "DNA_modifier_types.h"
42 #include "DNA_meshdata_types.h"
43 #include "DNA_mesh_types.h"
44 #include "DNA_scene_types.h"
45
46 #include "BLI_math.h"
47 #include "BLI_utildefines.h"
48
49 #include "BKE_shrinkwrap.h"
50 #include "BKE_DerivedMesh.h"
51 #include "BKE_lattice.h"
52
53 #include "BKE_deform.h"
54 #include "BKE_mesh.h"
55 #include "BKE_subsurf.h"
56 #include "BKE_mesh.h"
57 #include "BKE_editmesh.h"
58
59 /* for timing... */
60 #if 0
61 #  include "PIL_time.h"
62 #else
63 #  define TIMEIT_BENCH(expr, id) (expr)
64 #endif
65
66 /* Util macros */
67 #define OUT_OF_MEMORY() ((void)printf("Shrinkwrap: Out of memory\n"))
68
69 /* get derived mesh */
70 /* TODO is anyfunction that does this? returning the derivedFinal without we caring if its in edit mode or not? */
71 DerivedMesh *object_get_derived_final(Object *ob, bool for_render)
72 {
73         Mesh *me = ob->data;
74         BMEditMesh *em = me->edit_btmesh;
75
76         if (for_render) {
77                 /* TODO(sergey): use proper derived render here in the future. */
78                 return ob->derivedFinal;
79         }
80
81         if (em) {
82                 DerivedMesh *dm = em->derivedFinal;
83                 return dm;
84         }
85
86         return ob->derivedFinal;
87 }
88
89 /* Space transform */
90 void space_transform_from_matrixs(SpaceTransform *data, float local[4][4], float target[4][4])
91 {
92         float itarget[4][4];
93         invert_m4_m4(itarget, target);
94         mul_m4_m4m4(data->local2target, itarget, local);
95         invert_m4_m4(data->target2local, data->local2target);
96 }
97
98 void space_transform_apply(const SpaceTransform *data, float co[3])
99 {
100         mul_v3_m4v3(co, ((SpaceTransform *)data)->local2target, co);
101 }
102
103 void space_transform_invert(const SpaceTransform *data, float co[3])
104 {
105         mul_v3_m4v3(co, ((SpaceTransform *)data)->target2local, co);
106 }
107
108 static void space_transform_apply_normal(const SpaceTransform *data, float no[3])
109 {
110         mul_mat3_m4_v3(((SpaceTransform *)data)->local2target, no);
111         normalize_v3(no); /* TODO: could we just determine de scale value from the matrix? */
112 }
113
114 static void space_transform_invert_normal(const SpaceTransform *data, float no[3])
115 {
116         mul_mat3_m4_v3(((SpaceTransform *)data)->target2local, no);
117         normalize_v3(no); /* TODO: could we just determine de scale value from the matrix? */
118 }
119
120 /*
121  * Shrinkwrap to the nearest vertex
122  *
123  * it builds a kdtree of vertexs we can attach to and then
124  * for each vertex performs a nearest vertex search on the tree
125  */
126 static void shrinkwrap_calc_nearest_vertex(ShrinkwrapCalcData *calc)
127 {
128         int i;
129
130         BVHTreeFromMesh treeData = NULL_BVHTreeFromMesh;
131         BVHTreeNearest nearest  = NULL_BVHTreeNearest;
132
133
134         TIMEIT_BENCH(bvhtree_from_mesh_verts(&treeData, calc->target, 0.0, 2, 6), bvhtree_verts);
135         if (treeData.tree == NULL) {
136                 OUT_OF_MEMORY();
137                 return;
138         }
139
140         /* Setup nearest */
141         nearest.index = -1;
142         nearest.dist = FLT_MAX;
143 #ifndef __APPLE__
144 #pragma omp parallel for default(none) private(i) firstprivate(nearest) shared(treeData, calc) schedule(static)
145 #endif
146         for (i = 0; i < calc->numVerts; ++i) {
147                 float *co = calc->vertexCos[i];
148                 float tmp_co[3];
149                 float weight = defvert_array_find_weight_safe(calc->dvert, i, calc->vgroup);
150                 if (weight == 0.0f) {
151                         continue;
152                 }
153
154
155                 /* Convert the vertex to tree coordinates */
156                 if (calc->vert) {
157                         copy_v3_v3(tmp_co, calc->vert[i].co);
158                 }
159                 else {
160                         copy_v3_v3(tmp_co, co);
161                 }
162                 space_transform_apply(&calc->local2target, tmp_co);
163
164                 /* Use local proximity heuristics (to reduce the nearest search)
165                  *
166                  * If we already had an hit before.. we assume this vertex is going to have a close hit to that other vertex
167                  * so we can initiate the "nearest.dist" with the expected value to that last hit.
168                  * This will lead in prunning of the search tree. */
169                 if (nearest.index != -1)
170                         nearest.dist = len_squared_v3v3(tmp_co, nearest.co);
171                 else
172                         nearest.dist = FLT_MAX;
173
174                 BLI_bvhtree_find_nearest(treeData.tree, tmp_co, &nearest, treeData.nearest_callback, &treeData);
175
176
177                 /* Found the nearest vertex */
178                 if (nearest.index != -1) {
179                         /* Adjusting the vertex weight,
180                          * so that after interpolating it keeps a certain distance from the nearest position */
181                         if (nearest.dist > FLT_EPSILON) {
182                                 const float dist = sqrtf(nearest.dist);
183                                 weight *= (dist - calc->keepDist) / dist;
184                         }
185
186                         /* Convert the coordinates back to mesh coordinates */
187                         copy_v3_v3(tmp_co, nearest.co);
188                         space_transform_invert(&calc->local2target, tmp_co);
189
190                         interp_v3_v3v3(co, co, tmp_co, weight);  /* linear interpolation */
191                 }
192         }
193
194         free_bvhtree_from_mesh(&treeData);
195 }
196
197
198 /*
199  * This function raycast a single vertex and updates the hit if the "hit" is considered valid.
200  * Returns TRUE if "hit" was updated.
201  * Opts control whether an hit is valid or not
202  * Supported options are:
203  *      MOD_SHRINKWRAP_CULL_TARGET_FRONTFACE (front faces hits are ignored)
204  *      MOD_SHRINKWRAP_CULL_TARGET_BACKFACE (back faces hits are ignored)
205  */
206 int BKE_shrinkwrap_project_normal(char options, const float vert[3],
207                                   const float dir[3], const SpaceTransform *transf,
208                                   BVHTree *tree, BVHTreeRayHit *hit,
209                                   BVHTree_RayCastCallback callback, void *userdata)
210 {
211         /* don't use this because this dist value could be incompatible
212          * this value used by the callback for comparing prev/new dist values.
213          * also, at the moment there is no need to have a corrected 'dist' value */
214 // #define USE_DIST_CORRECT
215
216         float tmp_co[3], tmp_no[3];
217         const float *co, *no;
218         BVHTreeRayHit hit_tmp;
219
220         /* Copy from hit (we need to convert hit rays from one space coordinates to the other */
221         memcpy(&hit_tmp, hit, sizeof(hit_tmp));
222
223         /* Apply space transform (TODO readjust dist) */
224         if (transf) {
225                 copy_v3_v3(tmp_co, vert);
226                 space_transform_apply(transf, tmp_co);
227                 co = tmp_co;
228
229                 copy_v3_v3(tmp_no, dir);
230                 space_transform_apply_normal(transf, tmp_no);
231                 no = tmp_no;
232
233 #ifdef USE_DIST_CORRECT
234                 hit_tmp.dist *= mat4_to_scale(((SpaceTransform *)transf)->local2target);
235 #endif
236         }
237         else {
238                 co = vert;
239                 no = dir;
240         }
241
242         hit_tmp.index = -1;
243
244         BLI_bvhtree_ray_cast(tree, co, no, 0.0f, &hit_tmp, callback, userdata);
245
246         if (hit_tmp.index != -1) {
247                 /* invert the normal first so face culling works on rotated objects */
248                 if (transf) {
249                         space_transform_invert_normal(transf, hit_tmp.no);
250                 }
251
252                 if (options & (MOD_SHRINKWRAP_CULL_TARGET_FRONTFACE | MOD_SHRINKWRAP_CULL_TARGET_BACKFACE)) {
253                         /* apply backface */
254                         const float dot = dot_v3v3(dir, hit_tmp.no);
255                         if (((options & MOD_SHRINKWRAP_CULL_TARGET_FRONTFACE) && dot <= 0.0f) ||
256                             ((options & MOD_SHRINKWRAP_CULL_TARGET_BACKFACE)  && dot >= 0.0f))
257                         {
258                                 return FALSE; /* Ignore hit */
259                         }
260                 }
261
262                 if (transf) {
263                         /* Inverting space transform (TODO make coeherent with the initial dist readjust) */
264                         space_transform_invert(transf, hit_tmp.co);
265 #ifdef USE_DIST_CORRECT
266                         hit_tmp.dist = len_v3v3(vert, hit_tmp.co);
267 #endif
268                 }
269
270                 BLI_assert(hit_tmp.dist <= hit->dist);
271
272                 memcpy(hit, &hit_tmp, sizeof(hit_tmp));
273                 return TRUE;
274         }
275         return FALSE;
276 }
277
278
279 static void shrinkwrap_calc_normal_projection(ShrinkwrapCalcData *calc, bool forRender)
280 {
281         int i;
282
283         /* Options about projection direction */
284         const char use_normal   = calc->smd->shrinkOpts;
285         const float proj_limit_squared = calc->smd->projLimit * calc->smd->projLimit;
286         float proj_axis[3]      = {0.0f, 0.0f, 0.0f};
287
288         /* Raycast and tree stuff */
289
290         /** \note 'hit.dist' is kept in the targets space, this is only used
291          * for finding the best hit, to get the real dist,
292          * measure the len_v3v3() from the input coord to hit.co */
293         BVHTreeRayHit hit;
294         BVHTreeFromMesh treeData = NULL_BVHTreeFromMesh;
295
296         /* auxiliary target */
297         DerivedMesh *auxMesh    = NULL;
298         BVHTreeFromMesh auxData = NULL_BVHTreeFromMesh;
299         SpaceTransform local2aux;
300
301         /* If the user doesn't allows to project in any direction of projection axis
302          * then theres nothing todo. */
303         if ((use_normal & (MOD_SHRINKWRAP_PROJECT_ALLOW_POS_DIR | MOD_SHRINKWRAP_PROJECT_ALLOW_NEG_DIR)) == 0)
304                 return;
305
306
307         /* Prepare data to retrieve the direction in which we should project each vertex */
308         if (calc->smd->projAxis == MOD_SHRINKWRAP_PROJECT_OVER_NORMAL) {
309                 if (calc->vert == NULL) return;
310         }
311         else {
312                 /* The code supports any axis that is a combination of X,Y,Z
313                  * although currently UI only allows to set the 3 different axis */
314                 if (calc->smd->projAxis & MOD_SHRINKWRAP_PROJECT_OVER_X_AXIS) proj_axis[0] = 1.0f;
315                 if (calc->smd->projAxis & MOD_SHRINKWRAP_PROJECT_OVER_Y_AXIS) proj_axis[1] = 1.0f;
316                 if (calc->smd->projAxis & MOD_SHRINKWRAP_PROJECT_OVER_Z_AXIS) proj_axis[2] = 1.0f;
317
318                 normalize_v3(proj_axis);
319
320                 /* Invalid projection direction */
321                 if (len_squared_v3(proj_axis) < FLT_EPSILON) {
322                         return;
323                 }
324         }
325
326         if (calc->smd->auxTarget) {
327                 auxMesh = object_get_derived_final(calc->smd->auxTarget, forRender);
328                 if (!auxMesh)
329                         return;
330                 SPACE_TRANSFORM_SETUP(&local2aux, calc->ob, calc->smd->auxTarget);
331         }
332
333         /* After sucessufuly build the trees, start projection vertexs */
334         if (bvhtree_from_mesh_faces(&treeData, calc->target, 0.0, 4, 6) &&
335             (auxMesh == NULL || bvhtree_from_mesh_faces(&auxData, auxMesh, 0.0, 4, 6)))
336         {
337
338 #ifndef __APPLE__
339 #pragma omp parallel for private(i, hit) schedule(static)
340 #endif
341                 for (i = 0; i < calc->numVerts; ++i) {
342                         float *co = calc->vertexCos[i];
343                         float tmp_co[3], tmp_no[3];
344                         const float weight = defvert_array_find_weight_safe(calc->dvert, i, calc->vgroup);
345
346                         if (weight == 0.0f) {
347                                 continue;
348                         }
349
350                         if (calc->vert) {
351                                 /* calc->vert contains verts from derivedMesh  */
352                                 /* this coordinated are deformed by vertexCos only for normal projection (to get correct normals) */
353                                 /* for other cases calc->varts contains undeformed coordinates and vertexCos should be used */
354                                 if (calc->smd->projAxis == MOD_SHRINKWRAP_PROJECT_OVER_NORMAL) {
355                                         copy_v3_v3(tmp_co, calc->vert[i].co);
356                                         normal_short_to_float_v3(tmp_no, calc->vert[i].no);
357                                 }
358                                 else {
359                                         copy_v3_v3(tmp_co, co);
360                                         copy_v3_v3(tmp_no, proj_axis);
361                                 }
362                         }
363                         else {
364                                 copy_v3_v3(tmp_co, co);
365                                 copy_v3_v3(tmp_no, proj_axis);
366                         }
367
368
369                         hit.index = -1;
370                         hit.dist = 10000.0f; /* TODO: we should use FLT_MAX here, but sweepsphere code isn't prepared for that */
371
372                         /* Project over positive direction of axis */
373                         if (use_normal & MOD_SHRINKWRAP_PROJECT_ALLOW_POS_DIR) {
374
375                                 if (auxData.tree) {
376                                         BKE_shrinkwrap_project_normal(0, tmp_co, tmp_no,
377                                                                       &local2aux, auxData.tree, &hit,
378                                                                       auxData.raycast_callback, &auxData);
379                                 }
380
381                                 BKE_shrinkwrap_project_normal(calc->smd->shrinkOpts, tmp_co, tmp_no,
382                                                               &calc->local2target, treeData.tree, &hit,
383                                                               treeData.raycast_callback, &treeData);
384                         }
385
386                         /* Project over negative direction of axis */
387                         if (use_normal & MOD_SHRINKWRAP_PROJECT_ALLOW_NEG_DIR) {
388                                 float inv_no[3];
389                                 negate_v3_v3(inv_no, tmp_no);
390
391                                 if (auxData.tree) {
392                                         BKE_shrinkwrap_project_normal(0, tmp_co, inv_no,
393                                                                       &local2aux, auxData.tree, &hit,
394                                                                       auxData.raycast_callback, &auxData);
395                                 }
396
397                                 BKE_shrinkwrap_project_normal(calc->smd->shrinkOpts, tmp_co, inv_no,
398                                                               &calc->local2target, treeData.tree, &hit,
399                                                               treeData.raycast_callback, &treeData);
400                         }
401
402                         /* don't set the initial dist (which is more efficient),
403                          * because its calculated in the targets space, we want the dist in our own space */
404                         if (proj_limit_squared != 0.0f) {
405                                 if (len_squared_v3v3(hit.co, co) > proj_limit_squared) {
406                                         hit.index = -1;
407                                 }
408                         }
409
410                         if (hit.index != -1) {
411                                 madd_v3_v3v3fl(hit.co, hit.co, tmp_no, calc->keepDist);
412                                 interp_v3_v3v3(co, co, hit.co, weight);
413                         }
414                 }
415         }
416
417         /* free data structures */
418         free_bvhtree_from_mesh(&treeData);
419         free_bvhtree_from_mesh(&auxData);
420 }
421
422 /*
423  * Shrinkwrap moving vertexs to the nearest surface point on the target
424  *
425  * it builds a BVHTree from the target mesh and then performs a
426  * NN matches for each vertex
427  */
428 static void shrinkwrap_calc_nearest_surface_point(ShrinkwrapCalcData *calc)
429 {
430         int i;
431
432         BVHTreeFromMesh treeData = NULL_BVHTreeFromMesh;
433         BVHTreeNearest nearest  = NULL_BVHTreeNearest;
434
435         /* Create a bvh-tree of the given target */
436         bvhtree_from_mesh_faces(&treeData, calc->target, 0.0, 2, 6);
437         if (treeData.tree == NULL) {
438                 OUT_OF_MEMORY();
439                 return;
440         }
441
442         /* Setup nearest */
443         nearest.index = -1;
444         nearest.dist = FLT_MAX;
445
446
447         /* Find the nearest vertex */
448 #ifndef __APPLE__
449 #pragma omp parallel for default(none) private(i) firstprivate(nearest) shared(calc, treeData) schedule(static)
450 #endif
451         for (i = 0; i < calc->numVerts; ++i) {
452                 float *co = calc->vertexCos[i];
453                 float tmp_co[3];
454                 float weight = defvert_array_find_weight_safe(calc->dvert, i, calc->vgroup);
455                 if (weight == 0.0f) continue;
456
457                 /* Convert the vertex to tree coordinates */
458                 if (calc->vert) {
459                         copy_v3_v3(tmp_co, calc->vert[i].co);
460                 }
461                 else {
462                         copy_v3_v3(tmp_co, co);
463                 }
464                 space_transform_apply(&calc->local2target, tmp_co);
465
466                 /* Use local proximity heuristics (to reduce the nearest search)
467                  *
468                  * If we already had an hit before.. we assume this vertex is going to have a close hit to that other vertex
469                  * so we can initiate the "nearest.dist" with the expected value to that last hit.
470                  * This will lead in prunning of the search tree. */
471                 if (nearest.index != -1)
472                         nearest.dist = len_squared_v3v3(tmp_co, nearest.co);
473                 else
474                         nearest.dist = FLT_MAX;
475
476                 BLI_bvhtree_find_nearest(treeData.tree, tmp_co, &nearest, treeData.nearest_callback, &treeData);
477
478                 /* Found the nearest vertex */
479                 if (nearest.index != -1) {
480                         if (calc->smd->shrinkOpts & MOD_SHRINKWRAP_KEEP_ABOVE_SURFACE) {
481                                 /* Make the vertex stay on the front side of the face */
482                                 madd_v3_v3v3fl(tmp_co, nearest.co, nearest.no, calc->keepDist);
483                         }
484                         else {
485                                 /* Adjusting the vertex weight,
486                                  * so that after interpolating it keeps a certain distance from the nearest position */
487                                 float dist = sasqrt(nearest.dist);
488                                 if (dist > FLT_EPSILON) {
489                                         /* linear interpolation */
490                                         interp_v3_v3v3(tmp_co, tmp_co, nearest.co, (dist - calc->keepDist) / dist);
491                                 }
492                                 else {
493                                         copy_v3_v3(tmp_co, nearest.co);
494                                 }
495                         }
496
497                         /* Convert the coordinates back to mesh coordinates */
498                         space_transform_invert(&calc->local2target, tmp_co);
499                         interp_v3_v3v3(co, co, tmp_co, weight);  /* linear interpolation */
500                 }
501         }
502
503         free_bvhtree_from_mesh(&treeData);
504 }
505
506 /* Main shrinkwrap function */
507 void shrinkwrapModifier_deform(ShrinkwrapModifierData *smd, Object *ob, DerivedMesh *dm,
508                                float (*vertexCos)[3], int numVerts, bool forRender)
509 {
510
511         DerivedMesh *ss_mesh    = NULL;
512         ShrinkwrapCalcData calc = NULL_ShrinkwrapCalcData;
513
514         /* remove loop dependencies on derived meshes (TODO should this be done elsewhere?) */
515         if (smd->target == ob) smd->target = NULL;
516         if (smd->auxTarget == ob) smd->auxTarget = NULL;
517
518
519         /* Configure Shrinkwrap calc data */
520         calc.smd = smd;
521         calc.ob = ob;
522         calc.numVerts = numVerts;
523         calc.vertexCos = vertexCos;
524
525         /* DeformVertex */
526         calc.vgroup = defgroup_name_index(calc.ob, calc.smd->vgroup_name);
527         if (dm) {
528                 calc.dvert = dm->getVertDataArray(dm, CD_MDEFORMVERT);
529         }
530         else if (calc.ob->type == OB_LATTICE) {
531                 calc.dvert = BKE_lattice_deform_verts_get(calc.ob);
532         }
533
534
535         if (smd->target) {
536                 calc.target = object_get_derived_final(smd->target, forRender);
537
538                 /* TODO there might be several "bugs" on non-uniform scales matrixs
539                  * because it will no longer be nearest surface, not sphere projection
540                  * because space has been deformed */
541                 SPACE_TRANSFORM_SETUP(&calc.local2target, ob, smd->target);
542
543                 /* TODO: smd->keepDist is in global units.. must change to local */
544                 calc.keepDist = smd->keepDist;
545         }
546
547
548
549         calc.vgroup = defgroup_name_index(calc.ob, smd->vgroup_name);
550
551         if (dm != NULL && smd->shrinkType == MOD_SHRINKWRAP_PROJECT) {
552                 /* Setup arrays to get vertexs positions, normals and deform weights */
553                 calc.vert   = dm->getVertDataArray(dm, CD_MVERT);
554                 calc.dvert  = dm->getVertDataArray(dm, CD_MDEFORMVERT);
555
556                 /* Using vertexs positions/normals as if a subsurface was applied */
557                 if (smd->subsurfLevels) {
558                         SubsurfModifierData ssmd = {{NULL}};
559                         ssmd.subdivType = ME_CC_SUBSURF;        /* catmull clark */
560                         ssmd.levels     = smd->subsurfLevels;   /* levels */
561
562                         ss_mesh = subsurf_make_derived_from_derived(dm, &ssmd, NULL, (ob->mode & OB_MODE_EDIT) ? SUBSURF_IN_EDIT_MODE : 0);
563
564                         if (ss_mesh) {
565                                 calc.vert = ss_mesh->getVertDataArray(ss_mesh, CD_MVERT);
566                                 if (calc.vert) {
567                                         /* TRICKY: this code assumes subsurface will have the transformed original vertices
568                                          * in their original order at the end of the vert array. */
569                                         calc.vert = calc.vert + ss_mesh->getNumVerts(ss_mesh) - dm->getNumVerts(dm);
570                                 }
571                         }
572
573                         /* Just to make sure we are not leaving any memory behind */
574                         assert(ssmd.emCache == NULL);
575                         assert(ssmd.mCache == NULL);
576                 }
577         }
578
579         /* Projecting target defined - lets work! */
580         if (calc.target) {
581                 switch (smd->shrinkType) {
582                         case MOD_SHRINKWRAP_NEAREST_SURFACE:
583                                 TIMEIT_BENCH(shrinkwrap_calc_nearest_surface_point(&calc), deform_surface);
584                                 break;
585
586                         case MOD_SHRINKWRAP_PROJECT:
587                                 TIMEIT_BENCH(shrinkwrap_calc_normal_projection(&calc, forRender), deform_project);
588                                 break;
589
590                         case MOD_SHRINKWRAP_NEAREST_VERTEX:
591                                 TIMEIT_BENCH(shrinkwrap_calc_nearest_vertex(&calc), deform_vertex);
592                                 break;
593                 }
594         }
595
596         /* free memory */
597         if (ss_mesh)
598                 ss_mesh->release(ss_mesh);
599 }