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