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