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