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