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