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