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