Fixed a typo
[blender-staging.git] / source / blender / blenkernel / intern / shrinkwrap.c
1 /**
2  * shrinkwrap.c
3  *
4  * ***** BEGIN GPL LICENSE BLOCK *****
5  *
6  * This program is free software; you can redistribute it and/or
7  * modify it under the terms of the GNU General Public License
8  * as published by the Free Software Foundation; either version 2
9  * of the License, or (at your option) any later version.
10  *
11  * This program is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with this program; if not, write to the Free Software Foundation,
18  * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
19  *
20  * The Original Code is Copyright (C) Blender Foundation.
21  * All rights reserved.
22  *
23  * The Original Code is: all of this file.
24  *
25  * Contributor(s): AndrĂ© Pinto
26  *
27  * ***** END GPL LICENSE BLOCK *****
28  */
29 #include <string.h>
30 #include <float.h>
31 #include <math.h>
32 #include <memory.h>
33 #include <stdio.h>
34 #include <time.h>
35
36 #include "DNA_object_types.h"
37 #include "DNA_modifier_types.h"
38 #include "DNA_meshdata_types.h"
39
40 #include "BKE_shrinkwrap.h"
41 #include "BKE_DerivedMesh.h"
42 #include "BKE_utildefines.h"
43 #include "BKE_deform.h"
44 #include "BKE_cdderivedmesh.h"
45 #include "BKE_displist.h"
46 #include "BKE_global.h"
47
48 #include "BLI_arithb.h"
49 #include "BLI_kdtree.h"
50 #include "BLI_kdopbvh.h"
51
52 #include "RE_raytrace.h"
53 #include "MEM_guardedalloc.h"
54
55
56 /* Util macros */
57 #define TO_STR(a)       #a
58 #define JOIN(a,b)       a##b
59
60 #define OUT_OF_MEMORY() ((void)printf("Shrinkwrap: Out of memory\n"))
61
62 /* Benchmark macros */
63 #ifndef _WIN32
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 typedef void ( *Shrinkwrap_ForeachVertexCallback) (DerivedMesh *target, float *co, float *normal);
87
88 /* get derived mesh */
89 //TODO is anyfunction that does this? returning the derivedFinal witouth we caring if its in edit mode or not?
90 DerivedMesh *object_get_derived_final(Object *ob, CustomDataMask dataMask)
91 {
92         if (ob==G.obedit)
93         {
94                 DerivedMesh *final = NULL;
95                 editmesh_get_derived_cage_and_final(&final, dataMask);
96                 return final;
97         }
98         else
99                 return mesh_get_derived_final(ob, dataMask);
100 }
101
102 /* Space transform */
103 void space_transform_from_matrixs(SpaceTransform *data, float local[4][4], float target[4][4])
104 {
105         float itarget[4][4];
106         Mat4Invert(itarget, target);
107         Mat4MulSerie(data->local2target, itarget, local, 0, 0, 0, 0, 0, 0);
108         Mat4Invert(data->target2local, data->local2target);
109 }
110
111 void space_transform_apply(const SpaceTransform *data, float *co)
112 {
113         VecMat4MulVecfl(co, data->local2target, co);
114 }
115
116 void space_transform_invert(const SpaceTransform *data, float *co)
117 {
118         VecMat4MulVecfl(co, data->target2local, co);
119 }
120
121 void space_transform_apply_normal(const SpaceTransform *data, float *no)
122 {
123         Mat4Mul3Vecfl(data->local2target, no);
124         Normalize(no); // TODO: could we just determine de scale value from the matrix?
125 }
126
127 void space_transform_invert_normal(const SpaceTransform *data, float *no)
128 {
129         Mat4Mul3Vecfl(data->target2local, no);
130         Normalize(no); // TODO: could we just determine de scale value from the matrix?
131 }
132
133 /*
134  * Returns the squared distance between two given points
135  */
136 static float squared_dist(const float *a, const float *b)
137 {
138         float tmp[3];
139         VECSUB(tmp, a, b);
140         return INPR(tmp, tmp);
141 }
142
143 /* Main shrinkwrap function */
144 void shrinkwrapModifier_deform(ShrinkwrapModifierData *smd, Object *ob, DerivedMesh *dm, float (*vertexCos)[3], int numVerts)
145 {
146
147         ShrinkwrapCalcData calc = NULL_ShrinkwrapCalcData;
148
149         //remove loop dependencies on derived meshs (TODO should this be done elsewhere?)
150         if(smd->target == ob) smd->target = NULL;
151         if(smd->auxTarget == ob) smd->auxTarget = NULL;
152
153
154         //Configure Shrinkwrap calc data
155         calc.smd = smd;
156         calc.ob = ob;
157         calc.original = dm;
158         calc.numVerts = numVerts;
159         calc.vertexCos = vertexCos;
160
161         if(smd->target)
162         {
163                 //TODO currently we need a copy in case object_get_derived_final returns an emDM that does not defines getVertArray or getFace array
164                 calc.target = CDDM_copy( object_get_derived_final(smd->target, CD_MASK_BAREMESH) );
165
166                 if(!calc.target)
167                 {
168                         printf("Target derived mesh is null! :S\n");
169                 }
170
171                 //TODO there might be several "bugs" on non-uniform scales matrixs.. because it will no longer be nearest surface, not sphere projection
172                 //because space has been deformed
173                 space_transform_setup(&calc.local2target, ob, smd->target);
174
175                 calc.keepDist = smd->keepDist;  //TODO: smd->keepDist is in global units.. must change to local
176         }
177
178
179         //Projecting target defined - lets work!
180         if(calc.target)
181         {
182
183                 printf("Shrinkwrap (%s)%d over (%s)%d\n",
184                         calc.ob->id.name,                       calc.numVerts,
185                         calc.smd->target->id.name,      calc.target->getNumVerts(calc.target)
186                 );
187
188                 switch(smd->shrinkType)
189                 {
190                         case MOD_SHRINKWRAP_NEAREST_SURFACE:
191                                 BENCH(shrinkwrap_calc_nearest_surface_point(&calc));
192                         break;
193
194                         case MOD_SHRINKWRAP_PROJECT:
195                                 BENCH(shrinkwrap_calc_normal_projection(&calc));
196                         break;
197
198                         case MOD_SHRINKWRAP_NEAREST_VERTEX:
199                                 BENCH(shrinkwrap_calc_nearest_vertex(&calc));
200                         break;
201                 }
202         }
203
204         //free memory
205         if(calc.target)
206                 calc.target->release( calc.target );
207 }
208
209 /*
210  * Shrinkwrap to the nearest vertex
211  *
212  * it builds a kdtree of vertexs we can attach to and then
213  * for each vertex performs a nearest vertex search on the tree
214  */
215 void shrinkwrap_calc_nearest_vertex(ShrinkwrapCalcData *calc)
216 {
217         int i;
218         const int vgroup                 = get_named_vertexgroup_num(calc->ob, calc->smd->vgroup_name);
219         MDeformVert *const dvert = calc->original ? calc->original->getVertDataArray(calc->original, CD_MDEFORMVERT) : NULL;
220
221         BVHTreeFromMesh treeData = NULL_BVHTreeFromMesh;
222         BVHTreeNearest  nearest  = NULL_BVHTreeNearest;
223
224
225         BENCH(bvhtree_from_mesh_verts(&treeData, calc->target, 0.0, 2, 6));
226         if(treeData.tree == NULL) return OUT_OF_MEMORY();
227
228         //Setup nearest
229         nearest.index = -1;
230         nearest.dist = FLT_MAX;
231
232 #pragma omp parallel for default(none) private(i) firstprivate(nearest) shared(treeData,calc) schedule(static)
233         for(i = 0; i<calc->numVerts; ++i)
234         {
235                 float *co = calc->vertexCos[i];
236                 float tmp_co[3];
237                 float weight = vertexgroup_get_vertex_weight(dvert, i, vgroup);
238                 if(weight == 0.0f) continue;
239
240                 VECCOPY(tmp_co, co);
241                 space_transform_apply(&calc->local2target, tmp_co); //Convert the coordinates to the tree coordinates
242
243                 //Use local proximity heuristics (to reduce the nearest search)
244                 //
245                 //If we already had an hit before.. we assume this vertex is going to have a close hit to that other vertex
246                 //so we can initiate the "nearest.dist" with the expected value to that last hit.
247                 //This will lead in prunning of the search tree.
248                 if(nearest.index != -1)
249                         nearest.dist = squared_dist(tmp_co, nearest.co);
250                 else
251                         nearest.dist = FLT_MAX;
252
253                 BLI_bvhtree_find_nearest(treeData.tree, tmp_co, &nearest, treeData.nearest_callback, &treeData);
254
255
256                 //Found the nearest vertex
257                 if(nearest.index != -1)
258                 {
259                         //Adjusting the vertex weight, so that after interpolating it keeps a certain distance from the nearest position
260                         float dist = sasqrt(nearest.dist);
261                         if(dist > FLT_EPSILON) weight *= (dist - calc->keepDist)/dist;
262
263                         //Convert the coordinates back to mesh coordinates
264                         VECCOPY(tmp_co, nearest.co);
265                         space_transform_invert(&calc->local2target, tmp_co);
266
267                         VecLerpf(co, co, tmp_co, weight);       //linear interpolation
268                 }
269         }
270
271         free_bvhtree_from_mesh(&treeData);
272 }
273
274 /*
275  * This function raycast a single vertex and updates the hit if the "hit" is considered valid.
276  * Returns TRUE if "hit" was updated.
277  * Opts control whether an hit is valid or not
278  * Supported options are:
279  *      MOD_SHRINKWRAP_CULL_TARGET_FRONTFACE (front faces hits are ignored)
280  *      MOD_SHRINKWRAP_CULL_TARGET_BACKFACE (back faces hits are ignored)
281  */
282 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)
283 {
284         float tmp_co[3], tmp_no[3];
285         const float *co, *no;
286         BVHTreeRayHit hit_tmp;
287
288         //Copy from hit (we need to convert hit rays from one space coordinates to the other
289         memcpy( &hit_tmp, hit, sizeof(hit_tmp) );
290
291         //Apply space transform (TODO readjust dist)
292         if(transf)
293         {
294                 VECCOPY( tmp_co, vert );
295                 space_transform_apply( transf, tmp_co );
296                 co = tmp_co;
297
298                 VECCOPY( tmp_no, dir );
299                 space_transform_apply_normal( transf, tmp_no );
300                 no = tmp_no;
301
302                 hit_tmp.dist *= Mat4ToScalef( transf->local2target );
303         }
304         else
305         {
306                 co = vert;
307                 no = dir;
308         }
309
310         hit_tmp.index = -1;
311
312         BLI_bvhtree_ray_cast(tree, co, no, &hit_tmp, callback, userdata);
313
314         if(hit_tmp.index != -1)
315         {
316                 float dot = INPR( dir, hit_tmp.no);
317
318                 if(((options & MOD_SHRINKWRAP_CULL_TARGET_FRONTFACE) && dot <= 0.0f)
319                 || ((options & MOD_SHRINKWRAP_CULL_TARGET_BACKFACE) && dot >= 0.0f))
320                         return FALSE; //Ignore hit
321
322
323                 //Inverting space transform (TODO make coeherent with the initial dist readjust)
324                 if(transf)
325                 {
326                         space_transform_invert( transf, hit_tmp.co );
327                         space_transform_invert_normal( transf, hit_tmp.no );
328
329                         hit_tmp.dist = VecLenf( vert, hit_tmp.co );
330                 }
331
332                 memcpy(hit, &hit_tmp, sizeof(hit_tmp) );
333                 return TRUE;
334         }
335         return FALSE;
336 }
337
338
339 void shrinkwrap_calc_normal_projection(ShrinkwrapCalcData *calc)
340 {
341         int i;
342
343         //Options about projection direction
344         const char use_normal    = calc->smd->shrinkOpts;
345         float proj_axis[3] = {0.0f, 0.0f, 0.0f};
346         MVert *vert  = NULL; //Needed in case of vertex normal
347
348         //Vertex group data
349         const int vgroup                   = get_named_vertexgroup_num(calc->ob, calc->smd->vgroup_name);
350         const MDeformVert *dvert = calc->original ? calc->original->getVertDataArray(calc->original, CD_MDEFORMVERT) : NULL;
351
352
353         //Raycast and tree stuff
354         BVHTreeRayHit hit;
355         BVHTreeFromMesh treeData = NULL_BVHTreeFromMesh;        //target
356
357         //auxiliar target
358         DerivedMesh * aux_mesh = NULL;
359         BVHTreeFromMesh auxData= NULL_BVHTreeFromMesh;
360         SpaceTransform local2aux;
361
362
363         //Prepare data to retrieve the direction in which we should project each vertex
364         if(calc->smd->projAxis == MOD_SHRINKWRAP_PROJECT_OVER_NORMAL)
365         {
366                 vert = calc->original ? calc->original->getVertDataArray(calc->original, CD_MVERT) : NULL;
367                 if(vert) CDDM_calc_normals(calc->original);     //Maybe normals aren't yet calculated
368         }
369         else
370         {
371                 //The code supports any axis that is a combination of X,Y,Z.. altought currently UI only allows to set the 3 diferent axis
372                 if(calc->smd->projAxis & MOD_SHRINKWRAP_PROJECT_OVER_X_AXIS) proj_axis[0] = 1.0f;
373                 if(calc->smd->projAxis & MOD_SHRINKWRAP_PROJECT_OVER_Y_AXIS) proj_axis[1] = 1.0f;
374                 if(calc->smd->projAxis & MOD_SHRINKWRAP_PROJECT_OVER_Z_AXIS) proj_axis[2] = 1.0f;
375
376                 Normalize(proj_axis);
377         }
378
379         if(vert == NULL && (INPR(proj_axis, proj_axis) < FLT_EPSILON))
380         {
381                 printf("Shrinkwrap can't project witouth normal information");
382                 return;
383         }
384
385         //If the user doesn't allows to project in any direction of projection axis... then theres nothing todo.
386         if((use_normal & (MOD_SHRINKWRAP_PROJECT_ALLOW_POS_DIR | MOD_SHRINKWRAP_PROJECT_ALLOW_NEG_DIR)) == 0)
387                 return;
388
389
390         //Build target tree
391         BENCH(bvhtree_from_mesh_faces(&treeData, calc->target, calc->keepDist, 4, 6));
392         if(treeData.tree == NULL) return OUT_OF_MEMORY();
393
394         //Build auxiliar target
395         if(calc->smd->auxTarget)
396         {
397                 space_transform_setup( &local2aux, calc->ob, calc->smd->auxTarget);
398
399                 aux_mesh = CDDM_copy( object_get_derived_final(calc->smd->auxTarget, CD_MASK_BAREMESH) );               //TODO currently we need a copy in case object_get_derived_final returns an emDM that does not defines getVertArray or getFace array
400                 if(aux_mesh)
401                         BENCH(bvhtree_from_mesh_faces(&auxData, aux_mesh, 0.0, 4, 6));
402                 else
403                         printf("Auxiliar target finalDerived mesh is null\n");
404         }
405
406
407         //Now, everything is ready to project the vertexs!
408 //#pragma omp parallel for private(i,hit) schedule(static)
409         for(i = 0; i<calc->numVerts; ++i)
410         {
411                 float *co = calc->vertexCos[i];
412                 float tmp_co[3], tmp_no[3];
413                 float lim = 1000; //TODO: we should use FLT_MAX here, but sweepsphere code isnt prepared for that
414                 float weight = vertexgroup_get_vertex_weight(dvert, i, vgroup);
415
416                 if(weight == 0.0f) continue;
417
418                 VECCOPY(tmp_co, co);
419
420                 if(vert)
421                         NormalShortToFloat(tmp_no, vert[i].no);
422                 else
423                         VECCOPY( tmp_no, proj_axis );
424
425
426                 hit.index = -1;
427                 hit.dist = lim;
428
429
430                 //Project over positive direction of axis
431                 if(use_normal & MOD_SHRINKWRAP_PROJECT_ALLOW_POS_DIR)
432                 {
433
434                         if(auxData.tree)
435                                 normal_projection_project_vertex(0, tmp_co, tmp_no, &local2aux, auxData.tree, &hit, auxData.raycast_callback, &auxData);
436
437                         normal_projection_project_vertex(calc->smd->shrinkOpts, tmp_co, tmp_no, &calc->local2target, treeData.tree, &hit, treeData.raycast_callback, &treeData);
438                 }
439
440                 //Project over negative direction of axis
441                 if(use_normal & MOD_SHRINKWRAP_PROJECT_ALLOW_NEG_DIR)
442                 {
443                         float inv_no[3] = { -tmp_no[0], -tmp_no[1], -tmp_no[2] };
444
445
446                         if(auxData.tree)
447                                 normal_projection_project_vertex(0, tmp_co, inv_no, &local2aux, auxData.tree, &hit, auxData.raycast_callback, &auxData);
448
449                         normal_projection_project_vertex(calc->smd->shrinkOpts, tmp_co, inv_no, &calc->local2target, treeData.tree, &hit, treeData.raycast_callback, &treeData);
450                 }
451
452
453                 if(hit.index != -1)
454                 {
455                         VecLerpf(co, co, hit.co, weight);
456                 }
457         }
458
459
460         //free data structures
461
462         free_bvhtree_from_mesh(&treeData);
463         free_bvhtree_from_mesh(&auxData);
464
465         if(aux_mesh)
466                 aux_mesh->release(aux_mesh);
467 }
468
469 /*
470  * Shrinkwrap moving vertexs to the nearest surface point on the target
471  *
472  * it builds a BVHTree from the target mesh and then performs a
473  * NN matchs for each vertex
474  */
475 void shrinkwrap_calc_nearest_surface_point(ShrinkwrapCalcData *calc)
476 {
477         int i;
478
479         const int vgroup = get_named_vertexgroup_num(calc->ob, calc->smd->vgroup_name);
480         const MDeformVert *const dvert = calc->original ? calc->original->getVertDataArray(calc->original, CD_MDEFORMVERT) : NULL;
481
482         BVHTreeFromMesh treeData = NULL_BVHTreeFromMesh;
483         BVHTreeNearest  nearest  = NULL_BVHTreeNearest;
484
485
486
487         //Create a bvh-tree of the given target
488         BENCH(bvhtree_from_mesh_faces( &treeData, calc->target, 0.0, 2, 6));
489         if(treeData.tree == NULL) return OUT_OF_MEMORY();
490
491         //Setup nearest
492         nearest.index = -1;
493         nearest.dist = FLT_MAX;
494
495
496         //Find the nearest vertex 
497 #pragma omp parallel for default(none) private(i) firstprivate(nearest) shared(calc,treeData) schedule(static)
498         for(i = 0; i<calc->numVerts; ++i)
499         {
500                 float *co = calc->vertexCos[i];
501                 float tmp_co[3];
502                 float weight = vertexgroup_get_vertex_weight(dvert, i, vgroup);
503                 if(weight == 0.0f) continue;
504
505                 //Convert the vertex to tree coordinates
506                 VECCOPY(tmp_co, co);
507                 space_transform_apply(&calc->local2target, tmp_co);
508
509                 //Use local proximity heuristics (to reduce the nearest search)
510                 //
511                 //If we already had an hit before.. we assume this vertex is going to have a close hit to that other vertex
512                 //so we can initiate the "nearest.dist" with the expected value to that last hit.
513                 //This will lead in prunning of the search tree.
514                 if(nearest.index != -1)
515                         nearest.dist = squared_dist(tmp_co, nearest.co);
516                 else
517                         nearest.dist = FLT_MAX;
518
519                 BLI_bvhtree_find_nearest(treeData.tree, tmp_co, &nearest, treeData.nearest_callback, &treeData);
520
521                 //Found the nearest vertex
522                 if(nearest.index)
523                 {
524                         if(calc->smd->shrinkOpts & MOD_SHRINKWRAP_KEEP_ABOVE_SURFACE)
525                         {
526                                 //Make the vertex stay on the front side of the face
527                                 VECADDFAC(tmp_co, nearest.co, nearest.no, calc->keepDist);
528                         }
529                         else
530                         {
531                                 //Adjusting the vertex weight, so that after interpolating it keeps a certain distance from the nearest position
532                                 float dist = sasqrt( nearest.dist );
533                                 if(dist > FLT_EPSILON)
534                                         VecLerpf(tmp_co, tmp_co, nearest.co, (dist - calc->keepDist)/dist);     //linear interpolation
535                                 else
536                                         VECCOPY( tmp_co, nearest.co );
537                         }
538
539                         //Convert the coordinates back to mesh coordinates
540                         space_transform_invert(&calc->local2target, tmp_co);
541                         VecLerpf(co, co, tmp_co, weight);       //linear interpolation
542                 }
543         }
544
545
546         free_bvhtree_from_mesh(&treeData);
547 }
548