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