doxygen: blender/blenkernel tagged.
[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
30 /** \file blender/blenkernel/intern/shrinkwrap.c
31  *  \ingroup bke
32  */
33
34 #include <string.h>
35 #include <float.h>
36 #include <math.h>
37 #include <memory.h>
38 #include <stdio.h>
39 #include <time.h>
40 #include <assert.h>
41
42 #include "DNA_object_types.h"
43 #include "DNA_modifier_types.h"
44 #include "DNA_meshdata_types.h"
45 #include "DNA_mesh_types.h"
46 #include "DNA_scene_types.h"
47
48 #include "BLI_editVert.h"
49 #include "BLI_math.h"
50 #include "BLI_utildefines.h"
51
52 #include "BKE_shrinkwrap.h"
53 #include "BKE_DerivedMesh.h"
54 #include "BKE_lattice.h"
55
56 #include "BKE_deform.h"
57 #include "BKE_mesh.h"
58 #include "BKE_subsurf.h"
59
60 /* Util macros */
61 #define OUT_OF_MEMORY() ((void)printf("Shrinkwrap: Out of memory\n"))
62
63 /* Benchmark macros */
64 #if !defined(_WIN32) && 0
65
66 #include <sys/time.h>
67
68 #define BENCH(a)        \
69         do {                    \
70                 double _t1, _t2;                                \
71                 struct timeval _tstart, _tend;  \
72                 clock_t _clock_init = clock();  \
73                 gettimeofday ( &_tstart, NULL); \
74                 (a);                                                    \
75                 gettimeofday ( &_tend, NULL);   \
76                 _t1 = ( double ) _tstart.tv_sec + ( double ) _tstart.tv_usec/ ( 1000*1000 );    \
77                 _t2 = ( double )   _tend.tv_sec + ( double )   _tend.tv_usec/ ( 1000*1000 );    \
78                 printf("%s: %fs (real) %fs (cpu)\n", #a, _t2-_t1, (float)(clock()-_clock_init)/CLOCKS_PER_SEC);\
79         } while(0)
80
81 #else
82
83 #define BENCH(a)        (a)
84
85 #endif
86
87 typedef void ( *Shrinkwrap_ForeachVertexCallback) (DerivedMesh *target, float *co, float *normal);
88
89 /* get derived mesh */
90 //TODO is anyfunction that does this? returning the derivedFinal witouth we caring if its in edit mode or not?
91 DerivedMesh *object_get_derived_final(Object *ob)
92 {
93         Mesh *me= ob->data;
94         EditMesh *em = BKE_mesh_get_editmesh(me);
95
96         if(em) {
97                 DerivedMesh *dm = em->derivedFinal;
98                 BKE_mesh_end_editmesh(me, em);
99                 return dm;
100         }
101
102         return ob->derivedFinal;
103 }
104
105 /* Space transform */
106 void space_transform_from_matrixs(SpaceTransform *data, float local[4][4], float target[4][4])
107 {
108         float itarget[4][4];
109         invert_m4_m4(itarget, target);
110         mul_serie_m4(data->local2target, itarget, local, NULL, NULL, NULL, NULL, NULL, NULL);
111         invert_m4_m4(data->target2local, data->local2target);
112 }
113
114 void space_transform_apply(const SpaceTransform *data, float *co)
115 {
116         mul_v3_m4v3(co, ((SpaceTransform*)data)->local2target, co);
117 }
118
119 void space_transform_invert(const SpaceTransform *data, float *co)
120 {
121         mul_v3_m4v3(co, ((SpaceTransform*)data)->target2local, co);
122 }
123
124 static void space_transform_apply_normal(const SpaceTransform *data, float *no)
125 {
126         mul_mat3_m4_v3( ((SpaceTransform*)data)->local2target, no);
127         normalize_v3(no); // TODO: could we just determine de scale value from the matrix?
128 }
129
130 static void space_transform_invert_normal(const SpaceTransform *data, float *no)
131 {
132         mul_mat3_m4_v3(((SpaceTransform*)data)->target2local, no);
133         normalize_v3(no); // TODO: could we just determine de scale value from the matrix?
134 }
135
136 /*
137  * Returns the squared distance between two given points
138  */
139 static float squared_dist(const float *a, const float *b)
140 {
141         float tmp[3];
142         VECSUB(tmp, a, b);
143         return INPR(tmp, tmp);
144 }
145
146 /*
147  * Shrinkwrap to the nearest vertex
148  *
149  * it builds a kdtree of vertexs we can attach to and then
150  * for each vertex performs a nearest vertex search on the tree
151  */
152 static void shrinkwrap_calc_nearest_vertex(ShrinkwrapCalcData *calc)
153 {
154         int i;
155
156         BVHTreeFromMesh treeData = NULL_BVHTreeFromMesh;
157         BVHTreeNearest  nearest  = NULL_BVHTreeNearest;
158
159
160         BENCH(bvhtree_from_mesh_verts(&treeData, calc->target, 0.0, 2, 6));
161         if(treeData.tree == NULL)
162         {
163                 OUT_OF_MEMORY();
164                 return;
165         }
166
167         //Setup nearest
168         nearest.index = -1;
169         nearest.dist = FLT_MAX;
170 #ifndef __APPLE__
171 #pragma omp parallel for default(none) private(i) firstprivate(nearest) shared(treeData,calc) schedule(static)
172 #endif
173         for(i = 0; i<calc->numVerts; ++i)
174         {
175                 float *co = calc->vertexCos[i];
176                 float tmp_co[3];
177                 float weight = defvert_array_find_weight_safe(calc->dvert, i, calc->vgroup);
178                 if(weight == 0.0f) continue;
179
180
181                 //Convert the vertex to tree coordinates
182                 if(calc->vert)
183                 {
184                         VECCOPY(tmp_co, calc->vert[i].co);
185                 }
186                 else
187                 {
188                         VECCOPY(tmp_co, co);
189                 }
190                 space_transform_apply(&calc->local2target, tmp_co);
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                         interp_v3_v3v3(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 *= mat4_to_scale( ((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                 /* invert the normal first so face culling works on rotated objects */
265                 if(transf) {
266                         space_transform_invert_normal(transf, hit_tmp.no);
267                 }
268
269                 if (options & (MOD_SHRINKWRAP_CULL_TARGET_FRONTFACE|MOD_SHRINKWRAP_CULL_TARGET_BACKFACE)) {
270                         /* apply backface */
271                         const float dot= dot_v3v3(dir, hit_tmp.no);
272                         if(     ((options & MOD_SHRINKWRAP_CULL_TARGET_FRONTFACE) && dot <= 0.0f) ||
273                                 ((options & MOD_SHRINKWRAP_CULL_TARGET_BACKFACE) && dot >= 0.0f)
274                         ) {
275                                 return FALSE; /* Ignore hit */
276                         }
277                 }
278
279                 if(transf) {
280                         /* Inverting space transform (TODO make coeherent with the initial dist readjust) */
281                         space_transform_invert(transf, hit_tmp.co);
282                         hit_tmp.dist = len_v3v3((float *)vert, hit_tmp.co);
283                 }
284
285                 memcpy(hit, &hit_tmp, sizeof(hit_tmp) );
286                 return TRUE;
287         }
288         return FALSE;
289 }
290
291
292 static void shrinkwrap_calc_normal_projection(ShrinkwrapCalcData *calc)
293 {
294         int i;
295
296         //Options about projection direction
297         const char use_normal   = calc->smd->shrinkOpts;
298         float proj_axis[3]              = {0.0f, 0.0f, 0.0f};
299
300         //Raycast and tree stuff
301         BVHTreeRayHit hit;
302         BVHTreeFromMesh treeData= NULL_BVHTreeFromMesh;
303
304         //auxiliar target
305         DerivedMesh *auxMesh    = NULL;
306         BVHTreeFromMesh auxData = NULL_BVHTreeFromMesh;
307         SpaceTransform local2aux;
308
309         //If the user doesn't allows to project in any direction of projection axis
310         //then theres nothing todo.
311         if((use_normal & (MOD_SHRINKWRAP_PROJECT_ALLOW_POS_DIR | MOD_SHRINKWRAP_PROJECT_ALLOW_NEG_DIR)) == 0)
312                 return;
313
314
315         //Prepare data to retrieve the direction in which we should project each vertex
316         if(calc->smd->projAxis == MOD_SHRINKWRAP_PROJECT_OVER_NORMAL)
317         {
318                 if(calc->vert == NULL) return;
319         }
320         else
321         {
322                 //The code supports any axis that is a combination of X,Y,Z
323                 //altought currently UI only allows to set the 3 diferent axis
324                 if(calc->smd->projAxis & MOD_SHRINKWRAP_PROJECT_OVER_X_AXIS) proj_axis[0] = 1.0f;
325                 if(calc->smd->projAxis & MOD_SHRINKWRAP_PROJECT_OVER_Y_AXIS) proj_axis[1] = 1.0f;
326                 if(calc->smd->projAxis & MOD_SHRINKWRAP_PROJECT_OVER_Z_AXIS) proj_axis[2] = 1.0f;
327
328                 normalize_v3(proj_axis);
329
330                 //Invalid projection direction
331                 if(INPR(proj_axis, proj_axis) < FLT_EPSILON)
332                         return; 
333         }
334
335         if(calc->smd->auxTarget)
336         {
337                 auxMesh = object_get_derived_final(calc->smd->auxTarget);
338                 if(!auxMesh)
339                         return;
340                 space_transform_setup( &local2aux, calc->ob, calc->smd->auxTarget);
341         }
342
343         //After sucessufuly build the trees, start projection vertexs
344         if( bvhtree_from_mesh_faces(&treeData, calc->target, 0.0, 4, 6)
345         &&  (auxMesh == NULL || bvhtree_from_mesh_faces(&auxData, auxMesh, 0.0, 4, 6)))
346         {
347
348 #ifndef __APPLE__
349 #pragma omp parallel for private(i,hit) schedule(static)
350 #endif
351                 for(i = 0; i<calc->numVerts; ++i)
352                 {
353                         float *co = calc->vertexCos[i];
354                         float tmp_co[3], tmp_no[3];
355                         float weight = defvert_array_find_weight_safe(calc->dvert, i, calc->vgroup);
356
357                         if(weight == 0.0f) continue;
358
359                         if(calc->vert)
360                         {
361                                 /* calc->vert contains verts from derivedMesh  */
362                                 /* this coordinated are deformed by vertexCos only for normal projection (to get correct normals) */
363                                 /* for other cases calc->varts contains undeformed coordinates and vertexCos should be used */
364                                 if(calc->smd->projAxis == MOD_SHRINKWRAP_PROJECT_OVER_NORMAL) {
365                                         VECCOPY(tmp_co, calc->vert[i].co);
366                                         normal_short_to_float_v3(tmp_no, calc->vert[i].no);
367                                 } else {
368                                         VECCOPY(tmp_co, co);
369                                         VECCOPY(tmp_no, proj_axis);
370                                 }
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 && hit.index == -1)
394                         {
395                                 float inv_no[3];
396                                 negate_v3_v3(inv_no, tmp_no);
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                                 madd_v3_v3v3fl(hit.co, hit.co, tmp_no, calc->keepDist);
408                                 interp_v3_v3v3(co, co, hit.co, weight);
409                         }
410                 }
411         }
412
413         //free data structures
414         free_bvhtree_from_mesh(&treeData);
415         free_bvhtree_from_mesh(&auxData);
416 }
417
418 /*
419  * Shrinkwrap moving vertexs to the nearest surface point on the target
420  *
421  * it builds a BVHTree from the target mesh and then performs a
422  * NN matchs for each vertex
423  */
424 static void shrinkwrap_calc_nearest_surface_point(ShrinkwrapCalcData *calc)
425 {
426         int i;
427
428         BVHTreeFromMesh treeData = NULL_BVHTreeFromMesh;
429         BVHTreeNearest  nearest  = NULL_BVHTreeNearest;
430
431         //Create a bvh-tree of the given target
432         BENCH(bvhtree_from_mesh_faces( &treeData, calc->target, 0.0, 2, 6));
433         if(treeData.tree == NULL)
434         {
435                 OUT_OF_MEMORY();
436                 return;
437         }
438
439         //Setup nearest
440         nearest.index = -1;
441         nearest.dist = FLT_MAX;
442
443
444         //Find the nearest vertex
445 #ifndef __APPLE__
446 #pragma omp parallel for default(none) private(i) firstprivate(nearest) shared(calc,treeData) schedule(static)
447 #endif
448         for(i = 0; i<calc->numVerts; ++i)
449         {
450                 float *co = calc->vertexCos[i];
451                 float tmp_co[3];
452                 float weight = defvert_array_find_weight_safe(calc->dvert, i, calc->vgroup);
453                 if(weight == 0.0f) continue;
454
455                 //Convert the vertex to tree coordinates
456                 if(calc->vert)
457                 {
458                         VECCOPY(tmp_co, calc->vert[i].co);
459                 }
460                 else
461                 {
462                         VECCOPY(tmp_co, co);
463                 }
464                 space_transform_apply(&calc->local2target, tmp_co);
465
466                 //Use local proximity heuristics (to reduce the nearest search)
467                 //
468                 //If we already had an hit before.. we assume this vertex is going to have a close hit to that other vertex
469                 //so we can initiate the "nearest.dist" with the expected value to that last hit.
470                 //This will lead in prunning of the search tree.
471                 if(nearest.index != -1)
472                         nearest.dist = squared_dist(tmp_co, nearest.co);
473                 else
474                         nearest.dist = FLT_MAX;
475
476                 BLI_bvhtree_find_nearest(treeData.tree, tmp_co, &nearest, treeData.nearest_callback, &treeData);
477
478                 //Found the nearest vertex
479                 if(nearest.index != -1)
480                 {
481                         if(calc->smd->shrinkOpts & MOD_SHRINKWRAP_KEEP_ABOVE_SURFACE)
482                         {
483                                 //Make the vertex stay on the front side of the face
484                                 VECADDFAC(tmp_co, nearest.co, nearest.no, calc->keepDist);
485                         }
486                         else
487                         {
488                                 //Adjusting the vertex weight, so that after interpolating it keeps a certain distance from the nearest position
489                                 float dist = sasqrt( nearest.dist );
490                                 if(dist > FLT_EPSILON)
491                                         interp_v3_v3v3(tmp_co, tmp_co, nearest.co, (dist - calc->keepDist)/dist);       //linear interpolation
492                                 else
493                                         VECCOPY( tmp_co, nearest.co );
494                         }
495
496                         //Convert the coordinates back to mesh coordinates
497                         space_transform_invert(&calc->local2target, tmp_co);
498                         interp_v3_v3v3(co, co, tmp_co, weight); //linear interpolation
499                 }
500         }
501
502         free_bvhtree_from_mesh(&treeData);
503 }
504
505 /* Main shrinkwrap function */
506 void shrinkwrapModifier_deform(ShrinkwrapModifierData *smd, Object *ob, DerivedMesh *dm, float (*vertexCos)[3], int numVerts)
507 {
508
509         DerivedMesh *ss_mesh    = NULL;
510         ShrinkwrapCalcData calc = NULL_ShrinkwrapCalcData;
511
512         //remove loop dependencies on derived meshs (TODO should this be done elsewhere?)
513         if(smd->target == ob) smd->target = NULL;
514         if(smd->auxTarget == ob) smd->auxTarget = NULL;
515
516
517         //Configure Shrinkwrap calc data
518         calc.smd = smd;
519         calc.ob = ob;
520         calc.numVerts = numVerts;
521         calc.vertexCos = vertexCos;
522
523         //DeformVertex
524         calc.vgroup = defgroup_name_index(calc.ob, calc.smd->vgroup_name);
525         if(dm)
526         {
527                 calc.dvert = dm->getVertDataArray(dm, CD_MDEFORMVERT);
528         }
529         else if(calc.ob->type == OB_LATTICE)
530         {
531                 calc.dvert = lattice_get_deform_verts(calc.ob);
532         }
533
534
535         if(smd->target)
536         {
537                 calc.target = object_get_derived_final(smd->target);
538
539                 //TODO there might be several "bugs" on non-uniform scales matrixs
540                 //because it will no longer be nearest surface, not sphere projection
541                 //because space has been deformed
542                 space_transform_setup(&calc.local2target, ob, smd->target);
543
544                 //TODO: smd->keepDist is in global units.. must change to local
545                 calc.keepDist = smd->keepDist;
546         }
547
548
549
550         calc.vgroup = defgroup_name_index(calc.ob, smd->vgroup_name);
551
552         if(dm != NULL && smd->shrinkType == MOD_SHRINKWRAP_PROJECT)
553         {
554                 //Setup arrays to get vertexs positions, normals and deform weights
555                 calc.vert   = dm->getVertDataArray(dm, CD_MVERT);
556                 calc.dvert  = dm->getVertDataArray(dm, CD_MDEFORMVERT);
557
558                 //Using vertexs positions/normals as if a subsurface was applied 
559                 if(smd->subsurfLevels)
560                 {
561                         SubsurfModifierData ssmd= {{0}};
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 leaving 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));
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