2744212848361e4f66480f2b0ea23cf4bb876290
[blender.git] / source / blender / blenkernel / intern / shrinkwrap.c
1 /**
2  * shrinkwrap.c
3  *
4  * ***** BEGIN GPL LICENSE BLOCK *****
5  *
6  * This program is free software; you can redistribute it and/or
7  * modify it under the terms of the GNU General Public License
8  * as published by the Free Software Foundation; either version 2
9  * of the License, or (at your option) any later version.
10  *
11  * This program is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with this program; if not, write to the Free Software Foundation,
18  * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
19  *
20  * The Original Code is Copyright (C) Blender Foundation.
21  * All rights reserved.
22  *
23  * The Original Code is: all of this file.
24  *
25  * Contributor(s): AndrĂ© Pinto
26  *
27  * ***** END GPL LICENSE BLOCK *****
28  */
29 #include <string.h>
30 #include <float.h>
31 #include <math.h>
32 #include <memory.h>
33 #include <stdio.h>
34 #include <time.h>
35
36 #include "DNA_object_types.h"
37 #include "DNA_modifier_types.h"
38 #include "DNA_meshdata_types.h"
39
40 #include "BKE_shrinkwrap.h"
41 #include "BKE_DerivedMesh.h"
42 #include "BKE_utildefines.h"
43 #include "BKE_deform.h"
44 #include "BKE_cdderivedmesh.h"
45 #include "BKE_displist.h"
46 #include "BKE_global.h"
47
48 #include "BLI_arithb.h"
49 #include "BLI_kdtree.h"
50 #include "BLI_kdopbvh.h"
51
52 #include "RE_raytrace.h"
53 #include "MEM_guardedalloc.h"
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 1
64
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 #define BENCH_VAR(name)         clock_t JOIN(_bench_step,name) = 0, JOIN(_bench_total,name) = 0
82 #define BENCH_BEGIN(name)       JOIN(_bench_step, name) = clock()
83 #define BENCH_END(name)         JOIN(_bench_total,name) += clock() - JOIN(_bench_step,name)
84 #define BENCH_RESET(name)       JOIN(_bench_total, name) = 0
85 #define BENCH_REPORT(name)      printf("%s: %fms (cpu) \n", TO_STR(name), JOIN(_bench_total,name)*1000.0f/CLOCKS_PER_SEC)
86
87 #else
88
89 #define BENCH(a)        (a)
90 #define BENCH_VAR(name)
91 #define BENCH_BEGIN(name)
92 #define BENCH_END(name)
93 #define BENCH_RESET(name)
94 #define BENCH_REPORT(name)
95
96 #endif
97
98 typedef void ( *Shrinkwrap_ForeachVertexCallback) (DerivedMesh *target, float *co, float *normal);
99
100 /* get derived mesh */
101 //TODO is anyfunction that does this? returning the derivedFinal witouth we caring if its in edit mode or not?
102 DerivedMesh *object_get_derived_final(Object *ob, CustomDataMask dataMask)
103 {
104         if (ob==G.obedit)
105         {
106                 DerivedMesh *final = NULL;
107                 editmesh_get_derived_cage_and_final(&final, dataMask);
108                 return final;
109         }
110         else
111                 return mesh_get_derived_final(ob, dataMask);
112 }
113
114 /* Space transform */
115 void space_transform_from_matrixs(SpaceTransform *data, float local[4][4], float target[4][4])
116 {
117         float itarget[4][4];
118         Mat4Invert(itarget, target); //Invserse might be outdated
119         Mat4MulSerie(data->local2target, itarget, local, 0, 0, 0, 0, 0, 0);
120         Mat4Invert(data->target2local, data->local2target);
121 }
122
123 void space_transform_apply(const SpaceTransform *data, float *co)
124 {
125         VecMat4MulVecfl(co, data->local2target, co);
126 }
127
128 void space_transform_invert(const SpaceTransform *data, float *co)
129 {
130         VecMat4MulVecfl(co, data->target2local, co);
131 }
132
133 void space_transform_apply_normal(const SpaceTransform *data, float *no)
134 {
135         Mat4Mul3Vecfl(data->local2target, no);
136         Normalize(no); // TODO: could we just determine de scale value from the matrix?
137 }
138
139 void space_transform_invert_normal(const SpaceTransform *data, float *no)
140 {
141         Mat4Mul3Vecfl(data->target2local, no);
142         Normalize(no); // TODO: could we just determine de scale value from the matrix?
143 }
144
145 /*
146  * Returns the squared distance between two given points
147  */
148 static float squared_dist(const float *a, const float *b)
149 {
150         float tmp[3];
151         VECSUB(tmp, a, b);
152         return INPR(tmp, tmp);
153 }
154
155 /*
156  *
157  */
158 static void derivedmesh_mergeNearestPoints(DerivedMesh *dm, float mdist, BitSet skipVert)
159 {
160         if(mdist > 0.0f)
161         {
162                 int i, j, merged;
163                 int     numVerts = dm->getNumVerts(dm);
164                 int *translate_vert = MEM_mallocN( sizeof(int)*numVerts, "merge points array");
165
166                 MVert *vert = dm->getVertDataArray(dm, CD_MVERT);
167
168                 if(!translate_vert) return;
169
170                 merged = 0;
171                 for(i=0; i<numVerts; i++)
172                 {
173                         translate_vert[i] = i;
174
175                         if(skipVert && bitset_get(skipVert, i)) continue;
176
177                         for(j = 0; j<i; j++)
178                         {
179                                 if(skipVert && bitset_get(skipVert, j)) continue;
180                                 if(squared_dist(vert[i].co, vert[j].co) < mdist)
181                                 {
182                                         translate_vert[i] = j;
183                                         merged++;
184                                         break;
185                                 }
186                         }
187                 }
188
189                 //some vertexs were merged.. recalculate structure (edges and faces)
190                 if(merged > 0)
191                 {
192                         int     numFaces = dm->getNumFaces(dm);
193                         int freeVert;
194                         MFace *face = dm->getFaceDataArray(dm, CD_MFACE);
195
196
197                         //Adjust vertexs using the translation_table.. only translations to back indexs are allowed
198                         //which means t[i] <= i must always verify
199                         for(i=0, freeVert = 0; i<numVerts; i++)
200                         {
201                                 if(translate_vert[i] == i)
202                                 {
203                                         memcpy(&vert[freeVert], &vert[i], sizeof(*vert));
204                                         translate_vert[i] = freeVert++;
205                                 }
206                                 else translate_vert[i] = translate_vert[ translate_vert[i] ];
207                         }
208
209                         CDDM_lower_num_verts(dm, numVerts - merged);
210
211                         for(i=0; i<numFaces; i++)
212                         {
213                                 MFace *f = face+i;
214                                 f->v1 = translate_vert[f->v1];
215                                 f->v2 = translate_vert[f->v2];
216                                 f->v3 = translate_vert[f->v3];
217                                 //TODO be carefull with vertexs v4 being translated to 0
218                                 f->v4 = translate_vert[f->v4];
219                         }
220
221                         //TODO: maybe update edges could be done outside this function
222                         CDDM_calc_edges(dm);
223                         //CDDM_calc_normals(dm);
224                 }
225
226                 if(translate_vert) MEM_freeN( translate_vert );
227         }
228 }
229
230
231 /*
232  * This function removes Unused faces, vertexs and edges from calc->target
233  *
234  * This function may modify calc->final. As so no data retrieved from
235  * it before the call to this function  can be considered valid
236  * In case it creates a new DerivedMesh, the old calc->final is freed
237  */
238 //TODO memory checks on allocs
239 /*
240 static void shrinkwrap_removeUnused(ShrinkwrapCalcData *calc)
241 {
242         int i, t;
243
244         DerivedMesh *old = calc->final, *new = NULL;
245         MFace *new_face = NULL;
246         MVert *new_vert  = NULL;
247
248         int numVerts= old->getNumVerts(old);
249         MVert *vert = old->getVertDataArray(old, CD_MVERT);
250
251         int     numFaces= old->getNumFaces(old);
252         MFace *face = old->getFaceDataArray(old, CD_MFACE);
253
254         BitSet moved_verts = calc->moved;
255
256         //Arrays to translate to new vertexs indexs
257         int *vert_index = (int*)MEM_callocN(sizeof(int)*(numVerts), "shrinkwrap used verts");
258         BitSet used_faces = bitset_new(numFaces, "shrinkwrap used faces");
259         int numUsedFaces = 0;
260
261
262         //calculate which vertexs need to be used
263         //even unmoved vertices might need to be used if theres a face that needs it
264         //calc real number of faces, and vertices
265         //Count used faces
266         for(i=0; i<numFaces; i++)
267         {
268                 char res = 0;
269                 if(bitset_get(moved_verts, face[i].v1)) res++;
270                 if(bitset_get(moved_verts, face[i].v2)) res++;
271                 if(bitset_get(moved_verts, face[i].v3)) res++;
272                 if(face[i].v4 && bitset_get(moved_verts, face[i].v4)) res++;
273
274                 //Ignore a face were not a single vertice moved
275                 if(res == 0) continue;
276
277
278                 //Only 1 vertice moved.. (if its a quad.. remove the vertice oposite to it)
279                 if(res == 1 && face[i].v4)
280                 {
281                         if(bitset_get(moved_verts, face[i].v1))
282                         {
283                                 //remove vertex 3
284                                 face[i].v3 = face[i].v4;
285                         }
286                         else if(bitset_get(moved_verts, face[i].v2))
287                         {
288                                 //remove vertex 4
289                         }
290                         else if(bitset_get(moved_verts, face[i].v3))
291                         {
292                                 //remove vertex 1
293                                 face[i].v1 = face[i].v4;
294                         }
295                         else if(bitset_get(moved_verts, face[i].v4))
296                         {
297                                 //remove vertex 2
298                                 face[i].v2 = face[i].v3;
299                                 face[i].v3 = face[i].v4;
300                         }
301
302                         face[i].v4 = 0; //this quad turned on a tri
303                 }
304                 
305 #if 0
306                 if(face[i].v4 && res == 3)
307                 {
308                         if(!bitset_get(moved_verts, face[i].v1))
309                         {
310                                 face[i].v1 = face[i].v4;
311                         }
312                         else if(!bitset_get(moved_verts, face[i].v2))
313                         {
314                                 face[i].v2 = face[i].v3;
315                                 face[i].v3 = face[i].v4;
316                         }
317                         else if(!bitset_get(moved_verts, face[i].v3))
318                         {
319                                 face[i].v3 = face[i].v4;
320                         }
321
322                         face[i].v4 = 0; //this quad turned on a tri
323                 }
324 #endif
325
326                 bitset_set(used_faces, i);      //Mark face to maintain
327                 numUsedFaces++;
328
329                 //Mark vertices are needed
330                 vert_index[face[i].v1] = 1;
331                 vert_index[face[i].v2] = 1;
332                 vert_index[face[i].v3] = 1;
333                 if(face[i].v4) vert_index[face[i].v4] = 1;
334         }
335
336         //DP: Accumulate vertexs indexs.. (will calculate the new vertex index with a 1 offset)
337         for(i=1; i<numVerts; i++)
338                 vert_index[i] += vert_index[i-1];
339                 
340         
341         //Start creating the clean mesh
342         new = CDDM_new(vert_index[numVerts-1], 0, numUsedFaces);
343
344         //Copy vertexs (unused are are removed)
345         new_vert  = new->getVertDataArray(new, CD_MVERT);
346         for(i=0, t=0; i<numVerts; i++)
347         {
348                 
349                 if(vert_index[i] != t)
350                 {
351                         t = vert_index[i];
352                         memcpy(new_vert++, vert+i, sizeof(MVert));
353
354                         if(bitset_get(moved_verts, i))
355                                 bitset_set(moved_verts, t-1);
356                         else
357                                 bitset_unset(moved_verts, t-1);
358                 }
359         }
360
361         //Copy faces
362         new_face = new->getFaceDataArray(new, CD_MFACE);
363         for(i=0, t=0; i<numFaces; i++)
364         {
365                 if(bitset_get(used_faces, i))
366                 {
367                         memcpy(new_face, face+i, sizeof(MFace));
368                         //update vertices indexs
369                         new_face->v1 = vert_index[new_face->v1]-1;
370                         new_face->v2 = vert_index[new_face->v2]-1;
371                         new_face->v3 = vert_index[new_face->v3]-1;
372                         if(new_face->v4)
373                         {
374                                 new_face->v4 = vert_index[new_face->v4]-1;
375
376                                 //Ups translated vertex ended on 0 .. TODO fix this
377                                 if(new_face->v4 == 0)
378                                 {
379                                 }
380                         }                       
381                         new_face++;
382                 }
383         }
384
385         //Free memory
386         bitset_free(used_faces);
387         MEM_freeN(vert_index);
388         old->release(old);
389
390         //Update edges
391         CDDM_calc_edges(new);
392         CDDM_calc_normals(new);
393
394         calc->final = new;
395 }
396
397
398 void shrinkwrap_projectToCutPlane(ShrinkwrapCalcData *calc_data)
399 {
400         if(calc_data->smd->cutPlane && calc_data->moved)
401         {
402                 int i;
403                 int unmoved = 0;
404                 int numVerts= 0;
405                 MVert *vert = NULL;
406                 MVert *vert_unmoved = NULL;
407
408                 ShrinkwrapCalcData calc;
409                 memcpy(&calc, calc_data, sizeof(calc));
410
411                 calc.moved = 0;
412
413                 if(calc.smd->cutPlane)
414                 {
415                         //TODO currently we need a copy in case object_get_derived_final returns an emDM that does not defines getVertArray or getFace array
416                         calc.target = CDDM_copy( object_get_derived_final(calc.smd->cutPlane, CD_MASK_BAREMESH) );
417
418                         if(!calc.target)
419                         {
420                                 return;
421                         }
422
423                         space_transform_setup(&calc.local2target, calc.ob, calc.smd->cutPlane);
424                         calc.keptDist = 0;
425                 }
426
427
428                 //Make a mesh with the points we want to project
429                 numVerts = calc_data->final->getNumVerts(calc_data->final);
430
431                 unmoved = 0;
432                 for(i=0; i<numVerts; i++)
433                         if(!bitset_get(calc_data->moved, i))
434                                 unmoved++;
435
436                 calc.final = CDDM_new(unmoved, 0, 0);
437                 if(!calc.final) return;
438
439
440                 vert = calc_data->final->getVertDataArray(calc_data->final, CD_MVERT);
441                 vert_unmoved = calc.final->getVertDataArray(calc.final, CD_MVERT);
442
443                 for(i=0; i<numVerts; i++)
444                         if(!bitset_get(calc_data->moved, i))
445                                 memcpy(vert_unmoved++, vert+i, sizeof(*vert_unmoved));
446
447                 //use shrinkwrap projection
448                 shrinkwrap_calc_normal_projection(&calc);
449
450                 //Copy the points back to the mesh
451                 vert = calc_data->final->getVertDataArray(calc_data->final, CD_MVERT);
452                 vert_unmoved = calc.final->getVertDataArray(calc.final, CD_MVERT);
453                 for(i=0; i<numVerts; i++)
454                         if(!bitset_get(calc_data->moved, i))
455                                 memcpy(vert+i, vert_unmoved++, sizeof(*vert_unmoved) );
456
457                 //free memory
458                 calc.final->release(calc.final);
459                 calc.target->release(calc.target);
460         }       
461
462 }
463 */
464
465 /* Main shrinkwrap function */
466 /*
467 DerivedMesh *shrinkwrapModifier_do(ShrinkwrapModifierData *smd, Object *ob, DerivedMesh *dm, int useRenderParams, int isFinalCalc)
468 {
469
470         ShrinkwrapCalcData calc;
471         memset(&calc, 0, sizeof(calc));
472
473         //Init Shrinkwrap calc data
474         calc.smd = smd;
475
476         calc.ob = ob;
477         calc.original = dm;
478         calc.final = CDDM_copy(calc.original);
479
480         if(!calc.final)
481         {
482                 OUT_OF_MEMORY();
483                 return dm;
484         }
485
486         CDDM_calc_normals(calc.final);  //Normals maybe not be calculated yet
487
488         //remove loop dependencies on derived meshs (TODO should this be done elsewhere?)
489         if(smd->target == ob) smd->target = NULL;
490         if(smd->cutPlane == ob) smd->cutPlane = NULL;
491
492
493         if(smd->target)
494         {
495                 //TODO currently we need a copy in case object_get_derived_final returns an emDM that does not defines getVertArray or getFace array
496                 calc.target = CDDM_copy( object_get_derived_final(smd->target, CD_MASK_BAREMESH) );
497
498                 if(!calc.target)
499                 {
500                         printf("Target derived mesh is null! :S\n");
501                 }
502
503                 //TODO there might be several "bugs" on non-uniform scales matrixs.. because it will no longer be nearest surface, not sphere projection
504                 //because space has been deformed
505                 space_transform_setup(&calc.local2target, ob, smd->target);
506
507                 calc.keptDist = smd->keptDist;  //TODO: smd->keptDist is in global units.. must change to local
508         }
509
510         //Projecting target defined - lets work!
511         if(calc.target)
512         {
513
514                 printf("Shrinkwrap (%s)%d over (%s)%d\n",
515                         calc.ob->id.name,                       calc.final->getNumVerts(calc.final),
516                         calc.smd->target->id.name,      calc.target->getNumVerts(calc.target)
517                 );
518
519
520                 switch(smd->shrinkType)
521                 {
522                         case MOD_SHRINKWRAP_NEAREST_SURFACE:
523                                 BENCH(shrinkwrap_calc_nearest_surface_point(&calc));
524 //                              BENCH(shrinkwrap_calc_foreach_vertex(&calc, bruteforce_shrinkwrap_calc_nearest_surface_point));
525                         break;
526
527                         case MOD_SHRINKWRAP_NORMAL:
528
529                                 if(calc.smd->shrinkOpts & MOD_SHRINKWRAP_REMOVE_UNPROJECTED_FACES)
530                                 calc.moved = bitset_new( calc.final->getNumVerts(calc.final), "shrinkwrap bitset data");
531
532
533 //                              BENCH(shrinkwrap_calc_normal_projection_raytree(&calc));
534 //                              calc.final->release( calc.final );
535 //                              calc.final = CDDM_copy(calc.original);
536
537                                 BENCH(shrinkwrap_calc_normal_projection(&calc));
538 //                              BENCH(shrinkwrap_calc_foreach_vertex(&calc, bruteforce_shrinkwrap_calc_normal_projection));
539
540                                 if(calc.moved)
541                                 {
542                                         //Adjust vertxs that didn't moved (project to cut plane)
543                                         shrinkwrap_projectToCutPlane(&calc);
544
545                                         //Destroy faces, edges and stuff
546                                         shrinkwrap_removeUnused(&calc);
547
548                                         //Merge points that didn't moved
549                                         derivedmesh_mergeNearestPoints(calc.final, calc.smd->mergeDist, calc.moved);
550                                         bitset_free(calc.moved);
551                                 }
552                         break;
553
554                         case MOD_SHRINKWRAP_NEAREST_VERTEX:
555
556                                 BENCH(shrinkwrap_calc_nearest_vertex(&calc));
557 //                              BENCH(shrinkwrap_calc_foreach_vertex(&calc, bruteforce_shrinkwrap_calc_nearest_vertex));
558                         break;
559                 }
560
561                 //free derived mesh
562                 calc.target->release( calc.target );
563                 calc.target = NULL;
564         }
565
566         CDDM_calc_normals(calc.final);
567
568         return calc.final;
569 }
570 */
571
572 /* Main shrinkwrap function */
573 void shrinkwrapModifier_deform(ShrinkwrapModifierData *smd, Object *ob, DerivedMesh *dm, float (*vertexCos)[3], int numVerts)
574 {
575
576         ShrinkwrapCalcData calc = NULL_ShrinkwrapCalcData;
577
578         //Init Shrinkwrap calc data
579         calc.smd = smd;
580
581         calc.ob = ob;
582         calc.original = dm;
583
584         calc.numVerts = numVerts;
585         calc.vertexCos = vertexCos;
586
587         //remove loop dependencies on derived meshs (TODO should this be done elsewhere?)
588         if(smd->target == ob) smd->target = NULL;
589         if(smd->cutPlane == ob) smd->cutPlane = NULL;
590
591
592         if(smd->target)
593         {
594                 //TODO currently we need a copy in case object_get_derived_final returns an emDM that does not defines getVertArray or getFace array
595                 calc.target = CDDM_copy( object_get_derived_final(smd->target, CD_MASK_BAREMESH) );
596
597                 if(!calc.target)
598                 {
599                         printf("Target derived mesh is null! :S\n");
600                 }
601
602                 //TODO there might be several "bugs" on non-uniform scales matrixs.. because it will no longer be nearest surface, not sphere projection
603                 //because space has been deformed
604                 space_transform_setup(&calc.local2target, ob, smd->target);
605
606                 calc.keptDist = smd->keptDist;  //TODO: smd->keptDist is in global units.. must change to local
607         }
608
609         //Projecting target defined - lets work!
610         if(calc.target)
611         {
612
613                 printf("Shrinkwrap (%s)%d over (%s)%d\n",
614                         calc.ob->id.name,                       calc.numVerts,
615                         calc.smd->target->id.name,      calc.target->getNumVerts(calc.target)
616                 );
617
618
619                 switch(smd->shrinkType)
620                 {
621                         case MOD_SHRINKWRAP_NEAREST_SURFACE:
622                                 BENCH(shrinkwrap_calc_nearest_surface_point(&calc));
623                         break;
624
625                         case MOD_SHRINKWRAP_NORMAL:
626                                 BENCH(shrinkwrap_calc_normal_projection(&calc));
627                         break;
628
629                         case MOD_SHRINKWRAP_NEAREST_VERTEX:
630                                 BENCH(shrinkwrap_calc_nearest_vertex(&calc));
631                         break;
632                 }
633
634                 //free derived mesh
635                 calc.target->release( calc.target );
636                 calc.target = NULL;
637         }
638 }
639
640 /*
641  * Shrinkwrap to the nearest vertex
642  *
643  * it builds a kdtree of vertexs we can attach to and then
644  * for each vertex on performs a nearest vertex search on the tree
645  */
646 void shrinkwrap_calc_nearest_vertex(ShrinkwrapCalcData *calc)
647 {
648         int i;
649         int vgroup              = get_named_vertexgroup_num(calc->ob, calc->smd->vgroup_name);
650
651         BVHTreeFromMesh treeData = NULL_BVHTreeFromMesh;
652         BVHTreeNearest  nearest  = NULL_BVHTreeNearest;
653
654         MDeformVert *dvert = calc->original ? calc->original->getVertDataArray(calc->original, CD_MDEFORMVERT) : NULL;
655
656
657         BENCH(bvhtree_from_mesh_verts(&treeData, calc->target, 0.0, 2, 6));
658         if(treeData.tree == NULL) return OUT_OF_MEMORY();
659
660         //Setup nearest
661         nearest.index = -1;
662         nearest.dist = FLT_MAX;
663
664 //#pragma omp parallel for private(i) private(nearest) schedule(static)
665         for(i = 0; i<calc->numVerts; ++i)
666         {
667                 float *co = calc->vertexCos[i];
668                 int index;
669                 float tmp_co[3];
670                 float weight = vertexgroup_get_vertex_weight(dvert, i, vgroup);
671                 if(weight == 0.0f) continue;
672
673                 VECCOPY(tmp_co, co);
674                 space_transform_apply(&calc->local2target, tmp_co);
675
676                 //Use local proximity heuristics (to reduce the nearest search)
677                 if(nearest.index != -1)
678                         nearest.dist = squared_dist(tmp_co, nearest.co);
679                 else
680                         nearest.dist = FLT_MAX;
681
682                 index = BLI_bvhtree_find_nearest(treeData.tree, tmp_co, &nearest, treeData.nearest_callback, &treeData);
683
684                 if(index != -1)
685                 {
686                         float dist = nearest.dist;
687                         if(dist > 1e-5) weight *= (dist - calc->keptDist)/dist;
688
689                         VECCOPY(tmp_co, nearest.co);
690                         space_transform_invert(&calc->local2target, tmp_co);
691                         VecLerpf(co, co, tmp_co, weight);       //linear interpolation
692                 }
693         }
694
695         free_bvhtree_from_mesh(&treeData);
696 }
697
698 /*
699  * This function raycast a single vertex and updates the hit if the "hit" is considered valid.
700  * Returns TRUE if "hit" was updated.
701  * Opts control whether an hit is valid or not
702  * Supported options are:
703  *      MOD_SHRINKWRAP_CULL_TARGET_FRONTFACE (front faces hits are ignored)
704  *      MOD_SHRINKWRAP_CULL_TARGET_BACKFACE (back faces hits are ignored)
705  */
706 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)
707 {
708         float tmp_co[3], tmp_no[3];
709         const float *co, *no;
710         BVHTreeRayHit hit_tmp;
711
712         //Copy from hit (we need to convert hit rays from one space coordinates to the other
713         memcpy( &hit_tmp, hit, sizeof(hit_tmp) );
714
715         //Apply space transform (TODO readjust dist)
716         if(transf)
717         {
718                 VECCOPY( tmp_co, vert );
719                 space_transform_apply( transf, tmp_co );
720                 co = tmp_co;
721
722                 VECCOPY( tmp_no, dir );
723                 space_transform_apply_normal( transf, tmp_no );
724                 no = tmp_no;
725
726                 hit_tmp.dist *= Mat4ToScalef( transf->local2target );
727         }
728         else
729         {
730                 co = vert;
731                 no = dir;
732         }
733
734         hit_tmp.index = -1;
735
736         BLI_bvhtree_ray_cast(tree, co, no, &hit_tmp, callback, userdata);
737
738         if(hit_tmp.index != -1)
739         {
740                 float dot = INPR( dir, hit_tmp.no);
741
742                 if(((options & MOD_SHRINKWRAP_CULL_TARGET_FRONTFACE) && dot < 0)
743                 || ((options & MOD_SHRINKWRAP_CULL_TARGET_BACKFACE) && dot > 0))
744                         return FALSE; //Ignore hit
745
746
747                 //Inverting space transform (TODO make coeherent with the initial dist readjust)
748                 if(transf)
749                 {
750                         space_transform_invert( transf, hit_tmp.co );
751                         space_transform_invert_normal( transf, hit_tmp.no );
752
753                         hit_tmp.dist = VecLenf( vert, hit_tmp.co );
754                 }
755
756                 memcpy(hit, &hit_tmp, sizeof(hit_tmp) );
757                 return TRUE;
758         }
759         return FALSE;
760 }
761
762
763 void shrinkwrap_calc_normal_projection(ShrinkwrapCalcData *calc)
764 {
765         int i;
766         int vgroup              = get_named_vertexgroup_num(calc->ob, calc->smd->vgroup_name);
767         char use_normal = calc->smd->shrinkOpts;
768
769         //setup raytracing
770         BVHTreeFromMesh treeData = NULL_BVHTreeFromMesh;
771         BVHTreeRayHit   hit      = NULL_BVHTreeRayHit;
772
773         //cutTree
774         DerivedMesh * limit_mesh = NULL;
775         BVHTreeFromMesh limitData= NULL_BVHTreeFromMesh;
776         SpaceTransform local2cut;
777
778         MVert        *vert = calc->original ? calc->original->getVertDataArray(calc->original, CD_MVERT) : NULL;                //Needed because of vertex normal
779         MDeformVert *dvert = calc->original ? calc->original->getVertDataArray(calc->original, CD_MDEFORMVERT) : NULL;
780
781         if(vert == NULL)
782         {
783                 printf("Shrinkwrap cant normal project witouth normal information");
784                 return;
785         }
786         if((use_normal & (MOD_SHRINKWRAP_ALLOW_INVERTED_NORMAL | MOD_SHRINKWRAP_ALLOW_DEFAULT_NORMAL)) == 0)
787                 return; //Nothing todo
788
789         CDDM_calc_normals(calc->original);      //Normals maybe arent yet calculated
790
791         BENCH(bvhtree_from_mesh_faces(&treeData, calc->target, calc->keptDist, 4, 6));
792         if(treeData.tree == NULL) return OUT_OF_MEMORY();
793
794
795         if(calc->smd->cutPlane)
796         {
797                 space_transform_setup( &local2cut, calc->ob, calc->smd->cutPlane);
798
799                 //TODO currently we need a copy in case object_get_derived_final returns an emDM that does not defines getVertArray or getFace array
800                 limit_mesh = CDDM_copy( object_get_derived_final(calc->smd->cutPlane, CD_MASK_BAREMESH) );
801                 if(limit_mesh)
802                         BENCH(bvhtree_from_mesh_faces(&limitData, limit_mesh, 0.0, 4, 6));
803                 else
804                         printf("CutPlane finalDerived mesh is null\n");
805         }
806
807 //#pragma omp parallel for private(i) private(hit) schedule(static)
808         for(i = 0; i<calc->numVerts; ++i)
809         {
810                 float *co = calc->vertexCos[i];
811                 float tmp_co[3], tmp_no[3];
812                 float lim = 1000;               //TODO: we should use FLT_MAX here, but sweepsphere code isnt prepared for that
813                 float weight = vertexgroup_get_vertex_weight(dvert, i, vgroup);
814                 char moved = FALSE;
815
816                 if(weight == 0.0f) continue;
817
818                 VECCOPY(tmp_co, co);
819                 NormalShortToFloat(tmp_no, vert[i].no);
820
821
822                 hit.index = -1;
823                 hit.dist = lim;
824
825
826                 if(use_normal & MOD_SHRINKWRAP_ALLOW_DEFAULT_NORMAL)
827                 {
828
829                         if(limitData.tree)
830                                 normal_projection_project_vertex(0, tmp_co, tmp_no, &local2cut, limitData.tree, &hit, limitData.raycast_callback, &limitData);
831
832                         if(normal_projection_project_vertex(calc->smd->shrinkOpts, tmp_co, tmp_no, &calc->local2target, treeData.tree, &hit, treeData.raycast_callback, &treeData))
833                                 moved = TRUE;
834                 }
835
836
837                 if(use_normal & MOD_SHRINKWRAP_ALLOW_INVERTED_NORMAL)
838                 {
839                         float inv_no[3] = { -tmp_no[0], -tmp_no[1], -tmp_no[2] };
840
841
842                         if(limitData.tree)
843                                 normal_projection_project_vertex(0, tmp_co, inv_no, &local2cut, limitData.tree, &hit, limitData.raycast_callback, &limitData);
844
845                         if(normal_projection_project_vertex(calc->smd->shrinkOpts, tmp_co, inv_no, &calc->local2target, treeData.tree, &hit, treeData.raycast_callback, &treeData))
846                                 moved = TRUE;
847                 }
848
849
850                 if(hit.index != -1)
851                 {
852                         VecLerpf(co, co, hit.co, weight);
853
854                         if(moved && calc->moved)
855                                 bitset_set(calc->moved, i);
856                 }
857         }
858
859         free_bvhtree_from_mesh(&treeData);
860         free_bvhtree_from_mesh(&limitData);
861
862         if(limit_mesh)
863                 limit_mesh->release(limit_mesh);
864 }
865
866 /*
867  * Shrinkwrap moving vertexs to the nearest surface point on the target
868  *
869  * it builds a BVHTree from the target mesh and then performs a
870  * NN matchs for each vertex
871  */
872 void shrinkwrap_calc_nearest_surface_point(ShrinkwrapCalcData *calc)
873 {
874         int i;
875
876         const int vgroup                 = get_named_vertexgroup_num(calc->ob, calc->smd->vgroup_name);
877         const MDeformVert *dvert = calc->original ? calc->original->getVertDataArray(calc->original, CD_MDEFORMVERT) : NULL;
878
879         BVHTreeFromMesh treeData = NULL_BVHTreeFromMesh;
880         BVHTreeNearest  nearest  = NULL_BVHTreeNearest;
881
882
883
884         //Create a bvh-tree of the given target
885         BENCH(bvhtree_from_mesh_faces( &treeData, calc->target, 0.0, 2, 6));
886         if(treeData.tree == NULL) return OUT_OF_MEMORY();
887
888         //Setup nearest
889         nearest.index = -1;
890         nearest.dist = FLT_MAX;
891
892
893         //Find the nearest vertex 
894 //#pragma omp parallel for private(i) private(nearest) schedule(static)
895         for(i = 0; i<calc->numVerts; ++i)
896         {
897                 float *co = calc->vertexCos[i];
898                 int index;
899                 float tmp_co[3];
900                 float weight = vertexgroup_get_vertex_weight(dvert, i, vgroup);
901                 if(weight == 0.0f) continue;
902
903                 VECCOPY(tmp_co, co);
904                 space_transform_apply(&calc->local2target, tmp_co);
905
906                 if(nearest.index != -1)
907                 {
908                         nearest.dist = squared_dist(tmp_co, nearest.co);
909                 }
910                 else nearest.dist = FLT_MAX;
911
912                 index = BLI_bvhtree_find_nearest(treeData.tree, tmp_co, &nearest, treeData.nearest_callback, &treeData);
913
914                 if(index != -1)
915                 {
916                         if(calc->smd->shrinkOpts & MOD_SHRINKWRAP_KEPT_ABOVE_SURFACE)
917                         {
918                                 VECADDFAC(tmp_co, nearest.co, nearest.no, calc->keptDist);
919                         }
920                         else
921                         {
922                                 float dist = sasqrt( nearest.dist );
923                                 VecLerpf(tmp_co, tmp_co, nearest.co, (dist - calc->keptDist)/dist);     //linear interpolation
924                         }
925                         space_transform_invert(&calc->local2target, tmp_co);
926                         VecLerpf(co, co, tmp_co, weight);       //linear interpolation
927                 }
928         }
929
930         free_bvhtree_from_mesh(&treeData);
931 }
932