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