soc-2008-mxcurioni: foundations for Freestyle as a render layer - new UI buttons...
[blender.git] / source / blender / src / reeb.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. The Blender
10  * Foundation also sells licenses for use in proprietary software under
11  * the Blender License.  See http://www.blender.org/BL/ for information
12  * about this.
13  *
14  * This program is distributed in the hope that it will be useful,
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
17  * GNU General Public License for more details.
18  *
19  * You should have received a copy of the GNU General Public License
20  * along with this program; if not, write to the Free Software Foundation,
21  * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
22  *
23  * Contributor(s): Martin Poirier
24  *
25  * ***** END GPL LICENSE BLOCK *****
26  */
27  
28 #include <math.h>
29 #include <string.h> // for memcpy
30 #include <stdio.h>
31 #include <stdlib.h> // for qsort
32
33 #include "DNA_listBase.h"
34 #include "DNA_scene_types.h"
35 #include "DNA_space_types.h"
36 #include "DNA_meshdata_types.h"
37
38 #include "MEM_guardedalloc.h"
39
40 #include "BLI_blenlib.h"
41 #include "BLI_arithb.h"
42 #include "BLI_editVert.h"
43 #include "BLI_edgehash.h"
44
45 #include "BDR_editobject.h"
46
47 #include "BIF_editmesh.h"
48 #include "BIF_editarmature.h"
49 #include "BIF_interface.h"
50 #include "BIF_toolbox.h"
51 #include "BIF_graphics.h"
52
53 #include "BKE_global.h"
54 #include "BKE_utildefines.h"
55 #include "BKE_customdata.h"
56
57 #include "blendef.h"
58
59 #include "ONL_opennl.h"
60
61 #include "reeb.h"
62
63 /*
64  * Skeleton generation algorithm based on: 
65  * "Harmonic Skeleton for Realistic Character Animation"
66  * Gregoire Aujay, Franck Hetroy, Francis Lazarus and Christine Depraz
67  * SIGGRAPH 2007
68  * 
69  * Reeb graph generation algorithm based on: 
70  * "Robust On-line Computation of Reeb Graphs: Simplicity and Speed"
71  * Valerio Pascucci, Giorgio Scorzelli, Peer-Timo Bremer and Ajith Mascarenhas
72  * SIGGRAPH 2007
73  * 
74  * */
75
76 int mergeArcs(ReebGraph *rg, ReebArc *a0, ReebArc *a1);
77 int mergeConnectedArcs(ReebGraph *rg, ReebArc *a0, ReebArc *a1);
78 EditEdge * NextEdgeForVert(EditMesh *em, EditVert *v);
79
80 /***************************************** BUCKET UTILS **********************************************/
81
82 void addVertToBucket(EmbedBucket *b, float co[3])
83 {
84         b->nv++;
85         VecLerpf(b->p, b->p, co, 1.0f / b->nv);
86 }
87
88 void removeVertFromBucket(EmbedBucket *b, float co[3])
89 {
90         VecMulf(b->p, (float)b->nv);
91         VecSubf(b->p, b->p, co);
92         b->nv--;
93         VecMulf(b->p, 1.0f / (float)b->nv);
94 }
95
96 void mergeBuckets(EmbedBucket *bDst, EmbedBucket *bSrc)
97 {
98         if (bDst->nv > 0 && bSrc->nv > 0)
99         {
100                 bDst->nv += bSrc->nv;
101                 VecLerpf(bDst->p, bDst->p, bSrc->p, (float)bSrc->nv / (float)(bDst->nv));
102         }
103         else if (bSrc->nv > 0)
104         {
105                 bDst->nv = bSrc->nv;
106                 VECCOPY(bDst->p, bSrc->p);
107         }
108 }
109
110 void mergeArcBuckets(ReebArc *aDst, ReebArc *aSrc, float start, float end)
111 {
112         if (aDst->bcount > 0 && aSrc->bcount > 0)
113         {
114                 int indexDst = 0, indexSrc = 0;
115                 
116                 start = MAX3(start, aDst->buckets[0].val, aSrc->buckets[0].val);
117                 
118                 while(indexDst < aDst->bcount && aDst->buckets[indexDst].val < start)
119                 {
120                         indexDst++;
121                 }
122
123                 while(indexSrc < aSrc->bcount && aSrc->buckets[indexSrc].val < start)
124                 {
125                         indexSrc++;
126                 }
127                 
128                 for( ;  indexDst < aDst->bcount &&
129                                 indexSrc < aSrc->bcount &&
130                                 aDst->buckets[indexDst].val <= end &&
131                                 aSrc->buckets[indexSrc].val <= end
132                                 
133                          ;      indexDst++, indexSrc++)
134                 {
135                         mergeBuckets(aDst->buckets + indexDst, aSrc->buckets + indexSrc);
136                 }
137         }
138 }
139
140 void allocArcBuckets(ReebArc *arc)
141 {
142         int i;
143         float start = ceil(arc->v1->weight);
144         arc->bcount = (int)(floor(arc->v2->weight) - start) + 1;
145         
146         if (arc->bcount > 0)
147         {
148                 arc->buckets = MEM_callocN(sizeof(EmbedBucket) * arc->bcount, "embed bucket");
149                 
150                 for(i = 0; i < arc->bcount; i++)
151                 {
152                         arc->buckets[i].val = start + i;
153                 }
154         }
155         else
156         {
157                 arc->buckets = NULL;
158         }
159         
160 }
161
162 void resizeArcBuckets(ReebArc *arc)
163 {
164         EmbedBucket *oldBuckets = arc->buckets;
165         int oldBCount = arc->bcount;
166         
167         allocArcBuckets(arc);
168         
169         if (oldBCount != 0 && arc->bcount != 0)
170         {
171                 int oldStart = (int)oldBuckets[0].val;
172                 int oldEnd = (int)oldBuckets[oldBCount - 1].val;
173                 int newStart = (int)arc->buckets[0].val;
174                 int newEnd = (int)arc->buckets[arc->bcount - 1].val;
175                 int oldOffset = 0;
176                 int newOffset = 0;
177                 int len;
178                 
179                 if (oldStart < newStart)
180                 {
181                         oldOffset = newStart - oldStart;
182                 }
183                 else
184                 {
185                         newOffset = oldStart - newStart;
186                 }
187                 
188                 len = MIN2(oldEnd - (oldStart + oldOffset) + 1, newEnd - (newStart - newOffset) + 1);
189                 
190                 memcpy(arc->buckets + newOffset, oldBuckets + oldOffset, len * sizeof(EmbedBucket)); 
191         }
192
193         if (oldBuckets != NULL)
194         {
195                 MEM_freeN(oldBuckets);
196         }
197 }
198 /***************************************** UTILS **********************************************/
199
200 ReebEdge * copyEdge(ReebEdge *edge)
201 {
202         ReebEdge *newEdge = NULL;
203         
204         newEdge = MEM_callocN(sizeof(ReebEdge), "reeb edge");
205         memcpy(newEdge, edge, sizeof(ReebEdge));
206         
207         newEdge->next = NULL;
208         newEdge->prev = NULL;
209         
210         return newEdge;
211 }
212
213 void printArc(ReebArc *arc)
214 {
215         ReebEdge *edge;
216         printf("arc: (%i)%f -> (%i)%f\n", arc->v1->index, arc->v1->weight, arc->v2->index, arc->v2->weight);
217         
218         for(edge = arc->edges.first; edge ; edge = edge->next)
219         {
220                 printf("\tedge (%i, %i)\n", edge->v1->index, edge->v2->index);
221         }
222 }
223
224 void freeArc(ReebArc *arc)
225 {
226         BLI_freelistN(&arc->edges);
227         
228         if (arc->buckets)
229                 MEM_freeN(arc->buckets);
230         
231         MEM_freeN(arc);
232 }
233
234 void freeGraph(ReebGraph *rg)
235 {
236         ReebArc *arc;
237         ReebNode *node;
238         
239         // free nodes
240         for( node = rg->nodes.first; node; node = node->next )
241         {
242                 // Free adjacency lists
243                 if (node->arcs != NULL)
244                 {
245                         MEM_freeN(node->arcs);
246                 }
247         }
248         BLI_freelistN(&rg->nodes);
249         
250         // free arcs
251         arc = rg->arcs.first;
252         while( arc )
253         {
254                 ReebArc *next = arc->next;
255                 freeArc(arc);
256                 arc = next;
257         }
258         
259         // free edge map
260         BLI_edgehash_free(rg->emap, NULL);
261         
262         MEM_freeN(rg);
263 }
264
265 void repositionNodes(ReebGraph *rg)
266 {
267         ReebArc *arc = NULL;
268         ReebNode *node = NULL;
269         
270         // Reset node positions
271         for(node = rg->nodes.first; node; node = node->next)
272         {
273                 node->p[0] = node->p[1] = node->p[2] = 0;
274         }
275         
276         for(arc = rg->arcs.first; arc; arc = arc->next)
277         {
278                 if (arc->bcount > 0)
279                 {
280                         float p[3];
281                         
282                         VECCOPY(p, arc->buckets[0].p);
283                         VecMulf(p, 1.0f / arc->v1->degree);
284                         VecAddf(arc->v1->p, arc->v1->p, p);
285                         
286                         VECCOPY(p, arc->buckets[arc->bcount - 1].p);
287                         VecMulf(p, 1.0f / arc->v2->degree);
288                         VecAddf(arc->v2->p, arc->v2->p, p);
289                 }
290         }
291 }
292
293 void verifyNodeDegree(ReebGraph *rg)
294 {
295         ReebNode *node = NULL;
296         ReebArc *arc = NULL;
297
298         for(node = rg->nodes.first; node; node = node->next)
299         {
300                 int count = 0;
301                 for(arc = rg->arcs.first; arc; arc = arc->next)
302                 {
303                         if (arc->v1 == node || arc->v2 == node)
304                         {
305                                 count++;
306                         }
307                 }
308                 if (count != node->degree)
309                 {
310                         printf("degree error in node %i: expected %i got %i\n", node->index, count, node->degree);
311                 }
312         }
313 }
314
315 void verifyBuckets(ReebGraph *rg)
316 {
317 #ifdef DEBUG_REEB
318         ReebArc *arc = NULL;
319         for(arc = rg->arcs.first; arc; arc = arc->next)
320         {
321                 if (arc->bcount > 0)
322                 {
323                         int i;
324                         for(i = 0; i < arc->bcount; i++)
325                         {
326                                 if (arc->buckets[i].nv == 0)
327                                 {
328                                         printArc(arc);
329                                         printf("count error in bucket %i/%i\n", i+1, arc->bcount);
330                                 }
331                         }
332                         
333                         if (ceil(arc->v1->weight) < arc->buckets[0].val)
334                         {
335                                 printArc(arc);
336                                 printf("alloc error in first bucket: %f should be %f \n", arc->buckets[0].val, ceil(arc->v1->weight));
337                         }
338                         if (floor(arc->v2->weight) < arc->buckets[arc->bcount - 1].val)
339                         {
340                                 printArc(arc);
341                                 printf("alloc error in last bucket: %f should be %f \n", arc->buckets[arc->bcount - 1].val, floor(arc->v2->weight));
342                         }
343                 }
344         }
345 #endif
346 }
347
348 /************************************** ADJACENCY LIST *************************************************/
349
350 void addArcToNodeAdjacencyList(ReebNode *node, ReebArc *arc)
351 {
352         ReebArc **arclist;
353
354         for(arclist = node->arcs; *arclist; arclist++)
355         {       }
356         
357         *arclist = arc;
358 }
359
360 void buildAdjacencyList(ReebGraph *rg)
361 {
362         ReebNode *node = NULL;
363         ReebArc *arc = NULL;
364
365         for(node = rg->nodes.first; node; node = node->next)
366         {
367                 if (node->arcs != NULL)
368                 {
369                         MEM_freeN(node->arcs);
370                 }
371                 
372                 node->arcs = MEM_callocN((node->degree + 1) * sizeof(ReebArc*), "adjacency list");
373         }
374
375         for(arc = rg->arcs.first; arc; arc= arc->next)
376         {
377                 addArcToNodeAdjacencyList(arc->v1, arc);
378                 addArcToNodeAdjacencyList(arc->v2, arc);
379         }
380 }
381
382 int hasAdjacencyList(ReebGraph *rg)
383 {
384         ReebNode *node;
385         
386         for(node = rg->nodes.first; node; node = node->next)
387         {
388                 if (node->arcs == NULL)
389                 {
390                         return 0;
391                 }
392         }
393         
394         return 1;
395 }
396
397 int countConnectedArcs(ReebGraph *rg, ReebNode *node)
398 {
399         int count = 0;
400         
401         /* use adjacency list if present */
402         if (node->arcs)
403         {
404                 ReebArc **arcs;
405         
406                 for(arcs = node->arcs; *arcs; arcs++)
407                 {
408                         count++;
409                 }
410         }
411         else
412         {
413                 ReebArc *arc;
414                 for(arc = rg->arcs.first; arc; arc = arc->next)
415                 {
416                         if (arc->v1 == node || arc->v2 == node)
417                         {
418                                 count++;
419                         }
420                 }
421         }
422         
423         return count;
424 }
425
426 /****************************************** SMOOTHING **************************************************/
427
428 void postprocessGraph(ReebGraph *rg, char mode)
429 {
430         ReebArc *arc;
431         float fac1 = 0, fac2 = 1, fac3 = 0;
432
433         switch(mode)
434         {
435         case SKGEN_AVERAGE:
436                 fac1 = fac2 = fac3 = 1.0f / 3.0f;
437                 break;
438         case SKGEN_SMOOTH:
439                 fac1 = fac3 = 0.25f;
440                 fac2 = 0.5f;
441                 break;
442         case SKGEN_SHARPEN:
443                 fac1 = fac2 = -0.25f;
444                 fac2 = 1.5f;
445                 break;
446         default:
447                 error("Unknown post processing mode");
448                 return;
449         }
450         
451         for(arc = rg->arcs.first; arc; arc = arc->next)
452         {
453                 EmbedBucket *buckets = arc->buckets;
454                 int bcount = arc->bcount;
455                 int index;
456
457                 for(index = 1; index < bcount - 1; index++)
458                 {
459                         VecLerpf(buckets[index].p, buckets[index].p, buckets[index - 1].p, fac1 / (fac1 + fac2));
460                         VecLerpf(buckets[index].p, buckets[index].p, buckets[index + 1].p, fac3 / (fac1 + fac2 + fac3));
461                 }
462         }
463 }
464
465 /********************************************SORTING****************************************************/
466
467 int compareNodesWeight(void *vnode1, void *vnode2)
468 {
469         ReebNode *node1 = (ReebNode*)vnode1;
470         ReebNode *node2 = (ReebNode*)vnode2;
471         
472         if (node1->weight < node2->weight)
473         {
474                 return -1;
475         }
476         if (node1->weight > node2->weight)
477         {
478                 return 1;
479         }
480         else
481         {
482                 return 0;
483         }
484 }
485
486 void sortNodes(ReebGraph *rg)
487 {
488         BLI_sortlist(&rg->nodes, compareNodesWeight);
489 }
490
491 int compareArcsWeight(void *varc1, void *varc2)
492 {
493         ReebArc *arc1 = (ReebArc*)varc1;
494         ReebArc *arc2 = (ReebArc*)varc2;
495         
496         if (arc1->v1->weight < arc2->v1->weight)
497         {
498                 return -1;
499         }
500         if (arc1->v1->weight > arc2->v1->weight)
501         {
502                 return 1;
503         }
504         else
505         {
506                 return 0;
507         }
508 }
509
510 void sortArcs(ReebGraph *rg)
511 {
512         BLI_sortlist(&rg->arcs, compareArcsWeight);
513 }
514
515 /****************************************** FILTERING **************************************************/
516
517 int compareArcs(void *varc1, void *varc2)
518 {
519         ReebArc *arc1 = (ReebArc*)varc1;
520         ReebArc *arc2 = (ReebArc*)varc2;
521         float len1 = arc1->v2->weight - arc1->v1->weight;
522         float len2 = arc2->v2->weight - arc2->v1->weight;
523         
524         if (len1 < len2)
525         {
526                 return -1;
527         }
528         if (len1 > len2)
529         {
530                 return 1;
531         }
532         else
533         {
534                 return 0;
535         }
536 }
537
538 void filterArc(ReebGraph *rg, ReebNode *newNode, ReebNode *removedNode, ReebArc * srcArc, int merging)
539 {
540         ReebArc *arc = NULL, *nextArc = NULL;
541
542         /* first pass, merge buckets for arcs that spawned the two nodes into the source arc*/
543         for(arc = rg->arcs.first; arc; arc = arc->next)
544         {
545                 if (arc->v1 == srcArc->v1 && arc->v2 == srcArc->v2 && arc != srcArc)
546                 {
547                         mergeArcBuckets(srcArc, arc, srcArc->v1->weight, srcArc->v2->weight);
548                 }
549         }
550
551         /* second pass, replace removedNode by newNode, remove arcs that are collapsed in a loop */
552         arc = rg->arcs.first;
553         while(arc)
554         {
555                 nextArc = arc->next;
556                 
557                 if (arc->v1 == removedNode || arc->v2 == removedNode)
558                 {
559                         if (arc->v1 == removedNode)
560                         {
561                                 arc->v1 = newNode;
562                         }
563                         else
564                         {
565                                 arc->v2 = newNode;
566                         }
567
568                         // Remove looped arcs                   
569                         if (arc->v1 == arc->v2)
570                         {
571                                 // v1 or v2 was already newNode, since we're removing an arc, decrement degree
572                                 newNode->degree--;
573                                 
574                                 // If it's safeArc, it'll be removed later, so keep it for now
575                                 if (arc != srcArc)
576                                 {
577                                         BLI_remlink(&rg->arcs, arc);
578                                         freeArc(arc);
579                                 }
580                         }
581                         // Remove flipped arcs
582                         else if (arc->v1->weight > arc->v2->weight)
583                         {
584                                 // Decrement degree from the other node
585                                 OTHER_NODE(arc, newNode)->degree--;
586                                 
587                                 BLI_remlink(&rg->arcs, arc);
588                                 freeArc(arc);
589                         }
590                         else
591                         {
592                                 newNode->degree++; // incrementing degree since we're adding an arc
593
594                                 if (merging)
595                                 {
596                                         // resize bucket list
597                                         resizeArcBuckets(arc);
598                                         mergeArcBuckets(arc, srcArc, arc->v1->weight, arc->v2->weight);
599                                 }
600                         }
601                 }
602                 
603                 arc = nextArc;
604         }
605 }
606
607 void filterNullReebGraph(ReebGraph *rg)
608 {
609         ReebArc *arc = NULL, *nextArc = NULL;
610         
611         arc = rg->arcs.first;
612         while(arc)
613         {
614                 nextArc = arc->next;
615                 // Only collapse arcs too short to have any embed bucket
616                 if (arc->bcount == 0)
617                 {
618                         ReebNode *newNode = arc->v1;
619                         ReebNode *removedNode = arc->v2;
620                         float blend;
621                         
622                         blend = (float)newNode->degree / (float)(newNode->degree + removedNode->degree); // blending factors
623                         
624                         //newNode->weight = FloatLerpf(newNode->weight, removedNode->weight, blend);
625                         VecLerpf(newNode->p, newNode->p, removedNode->p, blend);
626                         
627                         filterArc(rg, newNode, removedNode, arc, 0);
628
629                         // Reset nextArc, it might have changed
630                         nextArc = arc->next;
631                         
632                         BLI_remlink(&rg->arcs, arc);
633                         freeArc(arc);
634                         
635                         BLI_freelinkN(&rg->nodes, removedNode);
636                 }
637                 
638                 arc = nextArc;
639         }
640 }
641
642 int filterInternalReebGraph(ReebGraph *rg, float threshold)
643 {
644         ReebArc *arc = NULL, *nextArc = NULL;
645         int value = 0;
646         
647         BLI_sortlist(&rg->arcs, compareArcs);
648
649         arc = rg->arcs.first;
650         while(arc)
651         {
652                 nextArc = arc->next;
653
654                 // Only collapse non-terminal arcs that are shorter than threshold
655                 if ((arc->v1->degree > 1 && arc->v2->degree > 1 && arc->v2->weight - arc->v1->weight < threshold))
656                 {
657                         ReebNode *newNode = NULL;
658                         ReebNode *removedNode = NULL;
659                         
660                         /* Keep the node with the highestn number of connected arcs */
661                         if (arc->v1->degree >= arc->v2->degree)
662                         {
663                                 newNode = arc->v1;
664                                 removedNode = arc->v2;
665                         }
666                         else
667                         {
668                                 newNode = arc->v2;
669                                 removedNode = arc->v1;
670                         }
671                         
672                         filterArc(rg, newNode, removedNode, arc, 1);
673
674                         // Reset nextArc, it might have changed
675                         nextArc = arc->next;
676                         
677                         BLI_remlink(&rg->arcs, arc);
678                         freeArc(arc);
679                         
680                         BLI_freelinkN(&rg->nodes, removedNode);
681                         value = 1;
682                 }
683                 
684                 arc = nextArc;
685         }
686         
687         return value;
688 }
689
690 int filterExternalReebGraph(ReebGraph *rg, float threshold)
691 {
692         ReebArc *arc = NULL, *nextArc = NULL;
693         int value = 0;
694         
695         BLI_sortlist(&rg->arcs, compareArcs);
696
697         arc = rg->arcs.first;
698         while(arc)
699         {
700                 nextArc = arc->next;
701
702                 // Only collapse terminal arcs that are shorter than threshold
703                 if ((arc->v1->degree == 1 || arc->v2->degree == 1) && arc->v2->weight - arc->v1->weight < threshold)
704                 {
705                         ReebNode *terminalNode = NULL;
706                         ReebNode *middleNode = NULL;
707                         ReebNode *newNode = NULL;
708                         ReebNode *removedNode = NULL;
709                         int merging = 0;
710                         
711                         // Assign terminal and middle nodes
712                         if (arc->v1->degree == 1)
713                         {
714                                 terminalNode = arc->v1;
715                                 middleNode = arc->v2;
716                         }
717                         else
718                         {
719                                 terminalNode = arc->v2;
720                                 middleNode = arc->v1;
721                         }
722                         
723                         // If middle node is a normal node, merge to terminal node
724                         if (middleNode->degree == 2)
725                         {
726                                 merging = 1;
727                                 newNode = terminalNode;
728                                 removedNode = middleNode;
729                         }
730                         // Otherwise, just plain remove of the arc
731                         else
732                         {
733                                 merging = 0;
734                                 newNode = middleNode;
735                                 removedNode = terminalNode;
736                         }
737                         
738                         // Merging arc
739                         if (merging)
740                         {
741                                 filterArc(rg, newNode, removedNode, arc, 1);
742                         }
743                         else
744                         {
745                                 // removing arc, so we need to decrease the degree of the remaining node
746                                 newNode->degree--;
747                         }
748
749                         // Reset nextArc, it might have changed
750                         nextArc = arc->next;
751                         
752                         BLI_remlink(&rg->arcs, arc);
753                         freeArc(arc);
754                         
755                         BLI_freelinkN(&rg->nodes, removedNode);
756                         value = 1;
757                 }
758                 
759                 arc = nextArc;
760         }
761         
762         return value;
763 }
764
765 /************************************** WEIGHT SPREADING ***********************************************/
766
767 int compareVerts( const void* a, const void* b )
768 {
769         EditVert *va = *(EditVert**)a;
770         EditVert *vb = *(EditVert**)b;
771         int value = 0;
772         
773         if (va->tmp.fp < vb->tmp.fp)
774         {
775                 value = -1;
776         }
777         else if (va->tmp.fp > vb->tmp.fp)
778         {
779                 value = 1;
780         }
781
782         return value;           
783 }
784
785 void spreadWeight(EditMesh *em)
786 {
787         EditVert **verts, *eve;
788         float lastWeight = 0.0f;
789         int totvert = BLI_countlist(&em->verts);
790         int i;
791         int work_needed = 1;
792         
793         verts = MEM_callocN(sizeof(EditVert*) * totvert, "verts array");
794         
795         for(eve = em->verts.first, i = 0; eve; eve = eve->next, i++)
796         {
797                 verts[i] = eve;
798         }
799         
800         while(work_needed == 1)
801         {
802                 work_needed = 0;
803                 qsort(verts, totvert, sizeof(EditVert*), compareVerts);
804                 
805                 for(i = 0; i < totvert; i++)
806                 {
807                         eve = verts[i];
808                         
809                         if (i == 0 || (eve->tmp.fp - lastWeight) > FLT_EPSILON)
810                         {
811                                 lastWeight = eve->tmp.fp;
812                         }
813                         else
814                         {
815                                 work_needed = 1;
816                                 eve->tmp.fp = lastWeight + FLT_EPSILON * 2;
817                                 lastWeight = eve->tmp.fp;
818                         }
819                 }
820         }
821         
822         MEM_freeN(verts);
823 }
824 /*********************************** GRAPH AS TREE FUNCTIONS *******************************************/
825
826 int subtreeDepth(ReebNode *node, ReebArc *rootArc)
827 {
828         int depth = 0;
829         
830         /* Base case, no arcs leading away */
831         if (node->arcs == NULL || *(node->arcs) == NULL)
832         {
833                 return 0;
834         }
835         else
836         {
837                 ReebArc ** pArc;
838
839                 for(pArc = node->arcs; *pArc; pArc++)
840                 {
841                         ReebArc *arc = *pArc;
842                         
843                         /* only arcs that go down the tree */
844                         if (arc != rootArc)
845                         {
846                                 ReebNode *newNode = OTHER_NODE(arc, node);
847                                 depth = MAX2(depth, subtreeDepth(newNode, arc));
848                         }
849                 }
850         }
851         
852         return depth + 1;
853 }
854
855 /*************************************** CYCLE DETECTION ***********************************************/
856
857 int detectCycle(ReebNode *node, ReebArc *srcArc)
858 {
859         int value = 0;
860         
861         if (node->flags == 0)
862         {
863                 ReebArc ** pArc;
864
865                 /* mark node as visited */
866                 node->flags = 1;
867
868                 for(pArc = node->arcs; *pArc && value == 0; pArc++)
869                 {
870                         ReebArc *arc = *pArc;
871                         
872                         /* don't go back on the source arc */
873                         if (arc != srcArc)
874                         {
875                                 value = detectCycle(OTHER_NODE(arc, node), arc);
876                         }
877                 }
878         }
879         else
880         {
881                 value = 1;
882         }
883         
884         return value;
885 }
886
887 int     isGraphCyclic(ReebGraph *rg)
888 {
889         ReebNode *node;
890         int value = 0;
891         
892         /* NEED TO CHECK IF ADJACENCY LIST EXIST */
893         
894         /* Mark all nodes as not visited */
895         for(node = rg->nodes.first; node; node = node->next)
896         {
897                 node->flags = 0;
898         }
899
900         /* detectCycles in subgraphs */ 
901         for(node = rg->nodes.first; node && value == 0; node = node->next)
902         {
903                 /* only for nodes in subgraphs that haven't been visited yet */
904                 if (node->flags == 0)
905                 {
906                         value = value || detectCycle(node, NULL);
907                 }               
908         }
909         
910         return value;
911 }
912
913 /******************************************** EXPORT ***************************************************/
914
915 void exportNode(FILE *f, char *text, ReebNode *node)
916 {
917         fprintf(f, "%s i:%i w:%f d:%i %f %f %f\n", text, node->index, node->weight, node->degree, node->p[0], node->p[1], node->p[2]);
918 }
919
920 void exportGraph(ReebGraph *rg, int count)
921 {
922 #ifdef DEBUG_REEB
923         ReebArc *arc;
924         char filename[128];
925         FILE *f;
926         
927         if (count == -1)
928         {
929                 sprintf(filename, "test.txt");
930         }
931         else
932         {
933                 sprintf(filename, "test%05i.txt", count);
934         }
935         f = fopen(filename, "w");
936
937         for(arc = rg->arcs.first; arc; arc = arc->next)
938         {
939                 int i;
940                 
941                 exportNode(f, "v1", arc->v1);
942                 
943                 for(i = 0; i < arc->bcount; i++)
944                 {
945                         fprintf(f, "b nv:%i %f %f %f\n", arc->buckets[i].nv, arc->buckets[i].p[0], arc->buckets[i].p[1], arc->buckets[i].p[2]);
946                 }
947                 
948                 exportNode(f, "v2", arc->v2);
949         }       
950         
951         fclose(f);
952 #endif
953 }
954
955 /***************************************** MAIN ALGORITHM **********************************************/
956
957 ReebArc * findConnectedArc(ReebGraph *rg, ReebArc *arc, ReebNode *v)
958 {
959         ReebArc *nextArc = arc->next;
960         
961         for(nextArc = rg->arcs.first; nextArc; nextArc = nextArc->next)
962         {
963                 if (arc != nextArc && (nextArc->v1 == v || nextArc->v2 == v))
964                 {
965                         break;
966                 }
967         }
968         
969         return nextArc;
970 }
971
972 void removeNormalNodes(ReebGraph *rg)
973 {
974         ReebArc *arc;
975         
976         // Merge degree 2 nodes
977         for(arc = rg->arcs.first; arc; arc = arc->next)
978         {
979                 while (arc->v1->degree == 2 || arc->v2->degree == 2)
980                 {
981                         // merge at v1
982                         if (arc->v1->degree == 2)
983                         {
984                                 ReebArc *nextArc = findConnectedArc(rg, arc, arc->v1);
985
986                                 // Merge arc only if needed
987                                 if (arc->v1 == nextArc->v2)
988                                 {                               
989                                         mergeConnectedArcs(rg, arc, nextArc);
990                                 }
991                                 // Otherwise, mark down vert
992                                 else
993                                 {
994                                         arc->v1->degree = 3;
995                                 }
996                         }
997                         
998                         // merge at v2
999                         if (arc->v2->degree == 2)
1000                         {
1001                                 ReebArc *nextArc = findConnectedArc(rg, arc, arc->v2);
1002                                 
1003                                 // Merge arc only if needed
1004                                 if (arc->v2 == nextArc->v1)
1005                                 {                               
1006                                         mergeConnectedArcs(rg, arc, nextArc);
1007                                 }
1008                                 // Otherwise, mark down vert
1009                                 else
1010                                 {
1011                                         arc->v2->degree = 3;
1012                                 }
1013                         }
1014                 }
1015         }
1016         
1017 }
1018
1019 int edgeEquals(ReebEdge *e1, ReebEdge *e2)
1020 {
1021         return (e1->v1 == e2->v1 && e1->v2 == e2->v2);
1022 }
1023
1024 ReebArc *nextArcMappedToEdge(ReebArc *arc, ReebEdge *e)
1025 {
1026         ReebEdge *nextEdge = NULL;
1027         ReebEdge *edge = NULL;
1028         ReebArc *result = NULL;
1029
1030         /* Find the ReebEdge in the edge list */
1031         for(edge = arc->edges.first; edge && !edgeEquals(edge, e); edge = edge->next)
1032         {       }
1033         
1034         nextEdge = edge->nextEdge;
1035         
1036         if (nextEdge != NULL)
1037         {
1038                 result = nextEdge->arc;
1039         }
1040
1041         return result;
1042 }
1043
1044 typedef enum {
1045         MERGE_LOWER,
1046         MERGE_HIGHER,
1047         MERGE_APPEND
1048 } MergeDirection;
1049
1050 void mergeArcEdges(ReebGraph *rg, ReebArc *aDst, ReebArc *aSrc, MergeDirection direction)
1051 {
1052         ReebEdge *e = NULL;
1053         
1054         if (direction == MERGE_APPEND)
1055         {
1056                 for(e = aSrc->edges.first; e; e = e->next)
1057                 {
1058                         e->arc = aDst; // Edge is stolen by new arc
1059                 }
1060                 
1061                 addlisttolist(&aDst->edges , &aSrc->edges);
1062         }
1063         else
1064         {
1065                 for(e = aSrc->edges.first; e; e = e->next)
1066                 {
1067                         ReebEdge *newEdge = copyEdge(e);
1068
1069                         newEdge->arc = aDst;
1070                         
1071                         BLI_addtail(&aDst->edges, newEdge);
1072                         
1073                         if (direction == MERGE_LOWER)
1074                         {
1075                                 void **p = BLI_edgehash_lookup_p(rg->emap, e->v1->index, e->v2->index);
1076                                 
1077                                 newEdge->nextEdge = e;
1078
1079                                 // if edge was the first in the list, point the edit edge to the new reeb edge instead.                                                 
1080                                 if (*p == e)
1081                                 {
1082                                         *p = (void*)newEdge;
1083                                 }
1084                                 // otherwise, advance in the list until the predecessor is found then insert it there
1085                                 else
1086                                 {
1087                                         ReebEdge *previous = (ReebEdge*)*p;
1088                                         
1089                                         while(previous->nextEdge != e)
1090                                         {
1091                                                 previous = previous->nextEdge;
1092                                         }
1093                                         
1094                                         previous->nextEdge = newEdge;
1095                                 }
1096                         }
1097                         else
1098                         {
1099                                 newEdge->nextEdge = e->nextEdge;
1100                                 e->nextEdge = newEdge;
1101                         }
1102                 }
1103         }
1104
1105
1106 // return 1 on full merge
1107 int mergeConnectedArcs(ReebGraph *rg, ReebArc *a0, ReebArc *a1)
1108 {
1109         int result = 0;
1110         ReebNode *removedNode = NULL;
1111         
1112         mergeArcEdges(rg, a0, a1, MERGE_APPEND);
1113         
1114         // Bring a0 to the combine length of both arcs
1115         if (a0->v2 == a1->v1)
1116         {
1117                 removedNode = a0->v2;
1118                 a0->v2 = a1->v2;
1119         }
1120         else if (a0->v1 == a1->v2)
1121         {
1122                 removedNode = a0->v1;
1123                 a0->v1 = a1->v1;
1124         }
1125         
1126         resizeArcBuckets(a0);
1127         // Merge a1 in a0
1128         mergeArcBuckets(a0, a1, a0->v1->weight, a0->v2->weight);
1129         
1130         // remove a1 from graph
1131         BLI_remlink(&rg->arcs, a1);
1132         freeArc(a1);
1133         
1134         BLI_freelinkN(&rg->nodes, removedNode);
1135         result = 1;
1136         
1137         return result;
1138 }
1139 // return 1 on full merge
1140 int mergeArcs(ReebGraph *rg, ReebArc *a0, ReebArc *a1)
1141 {
1142         int result = 0;
1143         // TRIANGLE POINTS DOWN
1144         if (a0->v1->weight == a1->v1->weight) // heads are the same
1145         {
1146                 if (a0->v2->weight == a1->v2->weight) // tails also the same, arcs can be totally merge together
1147                 {
1148                         mergeArcEdges(rg, a0, a1, MERGE_APPEND);
1149                         
1150                         mergeArcBuckets(a0, a1, a0->v1->weight, a0->v2->weight);
1151                         
1152                         // Adjust node degree
1153                         a1->v1->degree--;
1154                         a1->v2->degree--;
1155                         
1156                         // remove a1 from graph
1157                         BLI_remlink(&rg->arcs, a1);
1158                         
1159                         freeArc(a1);
1160                         result = 1;
1161                 }
1162                 else if (a0->v2->weight > a1->v2->weight) // a1->v2->weight is in the middle
1163                 {
1164                         mergeArcEdges(rg, a1, a0, MERGE_LOWER);
1165
1166                         // Adjust node degree
1167                         a0->v1->degree--;
1168                         a1->v2->degree++;
1169                         
1170                         mergeArcBuckets(a1, a0, a1->v1->weight, a1->v2->weight);
1171                         a0->v1 = a1->v2;
1172                         resizeArcBuckets(a0);
1173                 }
1174                 else // a0>n2 is in the middle
1175                 {
1176                         mergeArcEdges(rg, a0, a1, MERGE_LOWER);
1177                         
1178                         // Adjust node degree
1179                         a1->v1->degree--;
1180                         a0->v2->degree++;
1181                         
1182                         mergeArcBuckets(a0, a1, a0->v1->weight, a0->v2->weight);
1183                         a1->v1 = a0->v2;
1184                         resizeArcBuckets(a1);
1185                 }
1186         }
1187         // TRIANGLE POINTS UP
1188         else if (a0->v2->weight == a1->v2->weight) // tails are the same
1189         {
1190                 if (a0->v1->weight > a1->v1->weight) // a0->v1->weight is in the middle
1191                 {
1192                         mergeArcEdges(rg, a0, a1, MERGE_HIGHER);
1193                         
1194                         // Adjust node degree
1195                         a1->v2->degree--;
1196                         a0->v1->degree++;
1197                         
1198                         mergeArcBuckets(a0, a1, a0->v1->weight, a0->v2->weight);
1199                         a1->v2 = a0->v1;
1200                         resizeArcBuckets(a1);
1201                 }
1202                 else // a1->v1->weight is in the middle
1203                 {
1204                         mergeArcEdges(rg, a1, a0, MERGE_HIGHER);
1205
1206                         // Adjust node degree
1207                         a0->v2->degree--;
1208                         a1->v1->degree++;
1209
1210                         mergeArcBuckets(a1, a0, a1->v1->weight, a1->v2->weight);
1211                         a0->v2 = a1->v1;
1212                         resizeArcBuckets(a0);
1213                 }
1214         }
1215         else
1216         {
1217                 // Need something here (OR NOT)
1218         }
1219         
1220         return result;
1221 }
1222
1223 void glueByMergeSort(ReebGraph *rg, ReebArc *a0, ReebArc *a1, ReebEdge *e0, ReebEdge *e1)
1224 {
1225         int total = 0;
1226         while (total == 0 && a0 != a1 && a0 != NULL && a1 != NULL)
1227         {
1228                 total = mergeArcs(rg, a0, a1);
1229                 
1230                 if (total == 0) // if it wasn't a total merge, go forward
1231                 {
1232                         if (a0->v2->weight < a1->v2->weight)
1233                         {
1234                                 a0 = nextArcMappedToEdge(a0, e0);
1235                         }
1236                         else
1237                         {
1238                                 a1 = nextArcMappedToEdge(a1, e1);
1239                         }
1240                 }
1241         }
1242 }
1243
1244 void mergePaths(ReebGraph *rg, ReebEdge *e0, ReebEdge *e1, ReebEdge *e2)
1245 {
1246         ReebArc *a0, *a1, *a2;
1247         a0 = e0->arc;
1248         a1 = e1->arc;
1249         a2 = e2->arc;
1250         
1251         glueByMergeSort(rg, a0, a1, e0, e1);
1252         glueByMergeSort(rg, a0, a2, e0, e2);
1253
1254
1255 ReebNode * addNode(ReebGraph *rg, EditVert *eve, float weight)
1256 {
1257         ReebNode *node = NULL;
1258         
1259         node = MEM_callocN(sizeof(ReebNode), "reeb node");
1260         
1261         node->flags = 0; // clear flags on init
1262         node->arcs = NULL;
1263         node->degree = 0;
1264         node->weight = weight;
1265         node->index = rg->totnodes;
1266         VECCOPY(node->p, eve->co);      
1267         
1268         BLI_addtail(&rg->nodes, node);
1269         rg->totnodes++;
1270         
1271         return node;
1272 }
1273
1274 ReebEdge * createArc(ReebGraph *rg, ReebNode *node1, ReebNode *node2)
1275 {
1276         ReebEdge *edge;
1277         
1278         edge = BLI_edgehash_lookup(rg->emap, node1->index, node2->index);
1279         
1280         // Only add existing edges that haven't been added yet
1281         if (edge == NULL)
1282         {
1283                 ReebArc *arc;
1284                 ReebNode *v1, *v2;
1285                 float len, offset;
1286                 int i;
1287                 
1288                 arc = MEM_callocN(sizeof(ReebArc), "reeb arc");
1289                 edge = MEM_callocN(sizeof(ReebEdge), "reeb edge");
1290                 
1291                 arc->flags = 0; // clear flags on init
1292                 
1293                 if (node1->weight <= node2->weight)
1294                 {
1295                         v1 = node1;     
1296                         v2 = node2;     
1297                 }
1298                 else
1299                 {
1300                         v1 = node2;     
1301                         v2 = node1;     
1302                 }
1303                 
1304                 arc->v1 = v1;
1305                 arc->v2 = v2;
1306                 
1307                 // increase node degree
1308                 v1->degree++;
1309                 v2->degree++;
1310
1311                 BLI_edgehash_insert(rg->emap, node1->index, node2->index, edge);
1312                 
1313                 edge->arc = arc;
1314                 edge->nextEdge = NULL;
1315                 edge->v1 = v1;
1316                 edge->v2 = v2;
1317                 
1318                 BLI_addtail(&rg->arcs, arc);
1319                 BLI_addtail(&arc->edges, edge);
1320                 
1321                 /* adding buckets for embedding */
1322                 allocArcBuckets(arc);
1323                 
1324                 offset = arc->v1->weight;
1325                 len = arc->v2->weight - arc->v1->weight;
1326
1327 #if 0
1328                 /* This is the actual embedding filling described in the paper
1329                  * the problem is that it only works with really dense meshes
1330                  */
1331                 if (arc->bcount > 0)
1332                 {
1333                         addVertToBucket(&(arc->buckets[0]), arc->v1->co);
1334                         addVertToBucket(&(arc->buckets[arc->bcount - 1]), arc->v2->co);
1335                 }
1336 #else
1337                 for(i = 0; i < arc->bcount; i++)
1338                 {
1339                         float co[3];
1340                         float f = (arc->buckets[i].val - offset) / len;
1341                         
1342                         VecLerpf(co, v1->p, v2->p, f);
1343                         addVertToBucket(&(arc->buckets[i]), co);
1344                 }
1345 #endif
1346
1347         }
1348         
1349         return edge;
1350 }
1351
1352 void addTriangleToGraph(ReebGraph *rg, ReebNode * n1, ReebNode * n2, ReebNode * n3)
1353 {
1354         ReebEdge *re1, *re2, *re3;
1355         ReebEdge *e1, *e2, *e3;
1356         float len1, len2, len3;
1357         
1358         re1 = createArc(rg, n1, n2);
1359         re2 = createArc(rg, n2, n3);
1360         re3 = createArc(rg, n3, n1);
1361         
1362         len1 = (float)fabs(n1->weight - n2->weight);
1363         len2 = (float)fabs(n2->weight - n3->weight);
1364         len3 = (float)fabs(n3->weight - n1->weight);
1365         
1366         /* The rest of the algorithm assumes that e1 is the longest edge */
1367         
1368         if (len1 >= len2 && len1 >= len3)
1369         {
1370                 e1 = re1;
1371                 e2 = re2;
1372                 e3 = re3;
1373         }
1374         else if (len2 >= len1 && len2 >= len3)
1375         {
1376                 e1 = re2;
1377                 e2 = re1;
1378                 e3 = re3;
1379         }
1380         else
1381         {
1382                 e1 = re3;
1383                 e2 = re2;
1384                 e3 = re1;
1385         }
1386         
1387         /* And e2 is the lowest edge
1388          * If e3 is lower than e2, swap them
1389          */
1390         if (e3->v1->weight < e2->v1->weight)
1391         {
1392                 ReebEdge *etmp = e2;
1393                 e2 = e3;
1394                 e3 = etmp;
1395         }
1396         
1397         
1398         mergePaths(rg, e1, e2, e3);
1399 }
1400
1401 ReebGraph * generateReebGraph(EditMesh *em, int subdivisions)
1402 {
1403         ReebGraph *rg;
1404         struct DynamicList * dlist;
1405         EditVert *eve;
1406         EditFace *efa;
1407         int index;
1408         int totvert;
1409         int totfaces;
1410         
1411 #ifdef DEBUG_REEB
1412         int countfaces = 0;
1413 #endif
1414         
1415         rg = MEM_callocN(sizeof(ReebGraph), "reeb graph");
1416         
1417         rg->totnodes = 0;
1418         rg->emap = BLI_edgehash_new();
1419         
1420         totvert = BLI_countlist(&em->verts);
1421         totfaces = BLI_countlist(&em->faces);
1422         
1423         renormalizeWeight(em, 1.0f);
1424         
1425         /* Spread weight to minimize errors */
1426         spreadWeight(em);
1427
1428         renormalizeWeight(em, (float)subdivisions);
1429
1430         /* Adding vertice */
1431         for(index = 0, eve = em->verts.first; eve; index++, eve = eve->next)
1432         {
1433                 eve->hash = index;
1434                 eve->f2 = 0;
1435                 eve->tmp.p = addNode(rg, eve, eve->tmp.fp);
1436         }
1437         
1438         /* Temporarely convert node list to dynamic list, for indexed access */
1439         dlist = BLI_dlist_from_listbase(&rg->nodes);
1440         
1441         /* Adding face, edge per edge */
1442         for(efa = em->faces.first; efa; efa = efa->next)
1443         {
1444                 ReebNode *n1, *n2, *n3;
1445                 
1446                 n1 = (ReebNode*)BLI_dlist_find_link(dlist, efa->v1->hash);
1447                 n2 = (ReebNode*)BLI_dlist_find_link(dlist, efa->v2->hash);
1448                 n3 = (ReebNode*)BLI_dlist_find_link(dlist, efa->v3->hash);
1449                 
1450                 addTriangleToGraph(rg, n1, n2, n3);
1451                 
1452                 if (efa->v4)
1453                 {
1454                         ReebNode *n4 = (ReebNode*)efa->v4->tmp.p;
1455                         addTriangleToGraph(rg, n1, n3, n4);
1456                 }
1457
1458 #ifdef DEBUG_REEB
1459                 countfaces++;
1460                 if (countfaces % 100 == 0)
1461                 {
1462                         printf("face %i of %i\n", countfaces, totfaces);
1463                 }
1464 #endif
1465                 
1466                 
1467         }
1468         BLI_listbase_from_dlist(dlist, &rg->nodes);
1469         
1470         removeNormalNodes(rg);
1471         
1472         return rg;
1473 }
1474
1475 /***************************************** WEIGHT UTILS **********************************************/
1476
1477 void renormalizeWeight(EditMesh *em, float newmax)
1478 {
1479         EditVert *eve;
1480         float minimum, maximum, range;
1481         
1482         if (em == NULL || BLI_countlist(&em->verts) == 0)
1483                 return;
1484
1485         /* First pass, determine maximum and minimum */
1486         eve = em->verts.first;
1487         minimum = eve->tmp.fp;
1488         maximum = eve->tmp.fp;
1489         for(eve = em->verts.first; eve; eve = eve->next)
1490         {
1491                 maximum = MAX2(maximum, eve->tmp.fp);
1492                 minimum = MIN2(minimum, eve->tmp.fp);
1493         }
1494         
1495         range = maximum - minimum;
1496
1497         /* Normalize weights */
1498         for(eve = em->verts.first; eve; eve = eve->next)
1499         {
1500                 eve->tmp.fp = (eve->tmp.fp - minimum) / range * newmax;
1501         }
1502 }
1503
1504
1505 int weightFromLoc(EditMesh *em, int axis)
1506 {
1507         EditVert *eve;
1508         
1509         if (em == NULL || BLI_countlist(&em->verts) == 0 || axis < 0 || axis > 2)
1510                 return 0;
1511
1512         /* Copy coordinate in weight */
1513         for(eve = em->verts.first; eve; eve = eve->next)
1514         {
1515                 eve->tmp.fp = eve->co[axis];
1516         }
1517
1518         return 1;
1519 }
1520
1521 static float cotan_weight(float *v1, float *v2, float *v3)
1522 {
1523         float a[3], b[3], c[3], clen;
1524
1525         VecSubf(a, v2, v1);
1526         VecSubf(b, v3, v1);
1527         Crossf(c, a, b);
1528
1529         clen = VecLength(c);
1530
1531         if (clen == 0.0f)
1532                 return 0.0f;
1533         
1534         return Inpf(a, b)/clen;
1535 }
1536
1537 int weightToHarmonic(EditMesh *em)
1538 {
1539         NLboolean success;
1540         EditVert *eve;
1541         EditEdge *eed;
1542         EditFace *efa;
1543         int totvert = 0;
1544         int index;
1545         int rval;
1546         
1547         /* Find local extrema */
1548         for(eve = em->verts.first; eve; eve = eve->next)
1549         {
1550                 totvert++;
1551         }
1552
1553         /* Solve with openNL */
1554         
1555         nlNewContext();
1556
1557         nlSolverParameteri(NL_NB_VARIABLES, totvert);
1558
1559         nlBegin(NL_SYSTEM);
1560         
1561         /* Find local extrema */
1562         for(index = 0, eve = em->verts.first; eve; index++, eve = eve->next)
1563         {
1564                 EditEdge *eed;
1565                 int maximum = 1;
1566                 int minimum = 1;
1567                 
1568                 eve->hash = index; /* Assign index to vertex */
1569                 
1570                 NextEdgeForVert(NULL, NULL); /* Reset next edge */
1571                 for(eed = NextEdgeForVert(em, eve); eed && (maximum || minimum); eed = NextEdgeForVert(em, eve))
1572                 {
1573                         EditVert *eve2;
1574                         
1575                         if (eed->v1 == eve)
1576                         {
1577                                 eve2 = eed->v2;
1578                         }
1579                         else
1580                         {
1581                                 eve2 = eed->v1;
1582                         }
1583                         
1584                         /* Adjacent vertex is bigger, not a local maximum */
1585                         if (eve2->tmp.fp > eve->tmp.fp)
1586                         {
1587                                 maximum = 0;
1588                         }
1589                         /* Adjacent vertex is smaller, not a local minimum */
1590                         else if (eve2->tmp.fp < eve->tmp.fp)
1591                         {
1592                                 minimum = 0;
1593                         }
1594                 }
1595                 
1596                 if (maximum || minimum)
1597                 {
1598                         float w = eve->tmp.fp;
1599                         eve->f1 = 0;
1600                         nlSetVariable(0, index, w);
1601                         nlLockVariable(index);
1602                 }
1603                 else
1604                 {
1605                         eve->f1 = 1;
1606                 }
1607         }
1608         
1609         nlBegin(NL_MATRIX);
1610
1611         /* Zero edge weight */
1612         for(eed = em->edges.first; eed; eed = eed->next)
1613         {
1614                 eed->tmp.l = 0;
1615         }
1616         
1617         /* Add faces count to the edge weight */
1618         for(efa = em->faces.first; efa; efa = efa->next)
1619         {
1620                 efa->e1->tmp.l++;
1621                 efa->e2->tmp.l++;
1622                 efa->e3->tmp.l++;
1623         }
1624
1625         /* Add faces angle to the edge weight */
1626         for(efa = em->faces.first; efa; efa = efa->next)
1627         {
1628                 /* Angle opposite e1 */
1629                 float t1= cotan_weight(efa->v1->co, efa->v2->co, efa->v3->co) / efa->e2->tmp.l;
1630                 
1631                 /* Angle opposite e2 */
1632                 float t2 = cotan_weight(efa->v2->co, efa->v3->co, efa->v1->co) / efa->e3->tmp.l;
1633
1634                 /* Angle opposite e3 */
1635                 float t3 = cotan_weight(efa->v3->co, efa->v1->co, efa->v2->co) / efa->e1->tmp.l;
1636                 
1637                 int i1 = efa->v1->hash;
1638                 int i2 = efa->v2->hash;
1639                 int i3 = efa->v3->hash;
1640                 
1641                 nlMatrixAdd(i1, i1, t2+t3);
1642                 nlMatrixAdd(i2, i2, t1+t3);
1643                 nlMatrixAdd(i3, i3, t1+t2);
1644         
1645                 nlMatrixAdd(i1, i2, -t3);
1646                 nlMatrixAdd(i2, i1, -t3);
1647         
1648                 nlMatrixAdd(i2, i3, -t1);
1649                 nlMatrixAdd(i3, i2, -t1);
1650         
1651                 nlMatrixAdd(i3, i1, -t2);
1652                 nlMatrixAdd(i1, i3, -t2);
1653         }
1654         
1655         nlEnd(NL_MATRIX);
1656
1657         nlEnd(NL_SYSTEM);
1658
1659         success = nlSolveAdvanced(NULL, NL_TRUE);
1660
1661         if (success)
1662         {
1663                 rval = 1;
1664                 for(index = 0, eve = em->verts.first; eve; index++, eve = eve->next)
1665                 {
1666                         eve->tmp.fp = nlGetVariable(0, index);
1667                 }
1668         }
1669         else
1670         {
1671                 rval = 0;
1672         }
1673
1674         nlDeleteContext(nlGetCurrent());
1675
1676         return rval;
1677 }
1678
1679
1680 EditEdge * NextEdgeForVert(EditMesh *em, EditVert *v)
1681 {
1682         static EditEdge *e = NULL;
1683         
1684         /* Reset method, call with NULL mesh pointer */
1685         if (em == NULL)
1686         {
1687                 e = NULL;
1688                 return NULL;
1689         }
1690         
1691         /* first pass, start at the head of the list */
1692         if (e == NULL)
1693         {
1694                 e = em->edges.first;
1695         }
1696         /* subsequent passes, start on the next edge */
1697         else
1698         {
1699                 e = e->next;
1700         }
1701
1702         for( ; e ; e = e->next)
1703         {
1704                 if (e->v1 == v || e->v2 == v)
1705                 {
1706                         break;
1707                 }
1708         }       
1709         
1710         return e;
1711 }
1712
1713 int weightFromDistance(EditMesh *em)
1714 {
1715         EditVert *eve;
1716         int totedge = 0;
1717         int vCount = 0;
1718         
1719         if (em == NULL || BLI_countlist(&em->verts) == 0)
1720         {
1721                 return 0;
1722         }
1723         
1724         totedge = BLI_countlist(&em->edges);
1725         
1726         if (totedge == 0)
1727         {
1728                 return 0;
1729         }
1730
1731         /* Initialize vertice flags and find at least one selected vertex */
1732         for(eve = em->verts.first; eve && vCount == 0; eve = eve->next)
1733         {
1734                 eve->f1 = 0;
1735                 if (eve->f & SELECT)
1736                 {
1737                         vCount = 1;
1738                 }
1739         }
1740         
1741         if (vCount == 0)
1742         {
1743                 return 0; /* no selected vert, failure */
1744         }
1745         else
1746         {
1747                 EditVert *eve, *current_eve = NULL;
1748                 /* Apply dijkstra spf for each selected vert */
1749                 for(eve = em->verts.first; eve; eve = eve->next)
1750                 {
1751                         if (eve->f & SELECT)
1752                         {
1753                                 current_eve = eve;
1754                                 eve->f1 = 1;
1755                                 
1756                                 {
1757                                         EditEdge *eed = NULL;
1758                                         EditEdge *select_eed = NULL;
1759                                         EditEdge **edges = NULL;
1760                                         float    currentWeight = 0;
1761                                         int      eIndex = 0;
1762                                         
1763                                         edges = MEM_callocN(totedge * sizeof(EditEdge*), "Edges");
1764                                         
1765                                         /* Calculate edge weight and initialize edge flags */
1766                                         for(eed= em->edges.first; eed; eed= eed->next)
1767                                         {
1768                                                 eed->tmp.fp = VecLenf(eed->v1->co, eed->v2->co);
1769                                                 eed->f1 = 0;
1770                                         }
1771                                         
1772                                         do {
1773                                                 int i;
1774                                                 
1775                                                 current_eve->f1 = 1; /* mark vertex as selected */
1776                                                 
1777                                                 /* Add all new edges connected to current_eve to the list */
1778                                                 NextEdgeForVert(NULL, NULL); // Reset next edge
1779                                                 for(eed = NextEdgeForVert(em, current_eve); eed; eed = NextEdgeForVert(em, current_eve))
1780                                                 { 
1781                                                         if (eed->f1 == 0)
1782                                                         {
1783                                                                 edges[eIndex] = eed;
1784                                                                 eed->f1 = 1;
1785                                                                 eIndex++;
1786                                                         }
1787                                                 }
1788                                                 
1789                                                 /* Find next shortest edge */
1790                                                 select_eed = NULL;
1791                                                 for(i = 0; i < eIndex; i++)
1792                                                 {
1793                                                         eed = edges[i];
1794                                                         
1795                                                         if (eed->f1 != 2 && (eed->v1->f1 == 0 || eed->v2->f1 == 0)) /* eed is not selected yet and leads to a new node */
1796                                                         {
1797                                                                 float newWeight = 0;
1798                                                                 if (eed->v1->f1 == 1)
1799                                                                 {
1800                                                                         newWeight = eed->v1->tmp.fp + eed->tmp.fp;
1801                                                                 }
1802                                                                 else
1803                                                                 {
1804                                                                         newWeight = eed->v2->tmp.fp + eed->tmp.fp;
1805                                                                 }
1806                                                                 
1807                                                                 if (select_eed == NULL || newWeight < currentWeight) /* no selected edge or current smaller than selected */
1808                                                                 {
1809                                                                         currentWeight = newWeight;
1810                                                                         select_eed = eed;
1811                                                                 }
1812                                                         }
1813                                                 }
1814                                                 
1815                                                 if (select_eed != NULL)
1816                                                 {
1817                                                         select_eed->f1 = 2;
1818                                                         
1819                                                         if (select_eed->v1->f1 == 0) /* v1 is the new vertex */
1820                                                         {
1821                                                                 current_eve = select_eed->v1;
1822                                                         }
1823                                                         else /* otherwise, it's v2 */
1824                                                         {
1825                                                                 current_eve = select_eed->v2;
1826                                                         }                               
1827                                                         current_eve->tmp.fp = currentWeight;
1828                                                 }
1829                                         } while (select_eed != NULL);
1830                                         
1831                                         MEM_freeN(edges);
1832                                 }
1833                         }
1834                 }
1835         }
1836
1837         return 1;
1838 }
1839
1840 MCol MColFromWeight(EditVert *eve)
1841 {
1842         MCol col;
1843         col.a = 255;
1844         col.b = (char)(eve->tmp.fp * 255);
1845         col.g = 0;
1846         col.r = (char)((1.0f - eve->tmp.fp) * 255);
1847         return col;
1848 }
1849
1850 void weightToVCol(EditMesh *em)
1851 {
1852         EditFace *efa;
1853         MCol *mcol;
1854         if (!EM_vertColorCheck()) {
1855                 return;
1856         }
1857         
1858         for(efa=em->faces.first; efa; efa=efa->next) {
1859                 mcol = CustomData_em_get(&em->fdata, efa->data, CD_MCOL);
1860                                 
1861                 mcol[0] = MColFromWeight(efa->v1);
1862                 mcol[1] = MColFromWeight(efa->v2);
1863                 mcol[2] = MColFromWeight(efa->v3);
1864
1865                 if(efa->v4) {
1866                         mcol[3] = MColFromWeight(efa->v4);
1867                 }
1868         }
1869 }
1870
1871 /****************************************** BUCKET ITERATOR **************************************************/
1872
1873 void initArcIterator(ReebArcIterator *iter, ReebArc *arc, ReebNode *head)
1874 {
1875         iter->arc = arc;
1876         
1877         if (head == arc->v1)
1878         {
1879                 iter->start = 0;
1880                 iter->end = arc->bcount - 1;
1881                 iter->stride = 1;
1882         }
1883         else
1884         {
1885                 iter->start = arc->bcount - 1;
1886                 iter->end = 0;
1887                 iter->stride = -1;
1888         }
1889         
1890         iter->index = iter->start - iter->stride;
1891 }
1892
1893 void initArcIterator2(ReebArcIterator *iter, ReebArc *arc, int start, int end)
1894 {
1895         iter->arc = arc;
1896         
1897         iter->start = start;
1898         iter->end = end;
1899         
1900         if (end > start)
1901         {
1902                 iter->stride = 1;
1903         }
1904         else
1905         {
1906                 iter->stride = -1;
1907         }
1908
1909         iter->index = iter->start - iter->stride;
1910 }
1911
1912 EmbedBucket * nextBucket(ReebArcIterator *iter)
1913 {
1914         EmbedBucket *result = NULL;
1915         
1916         if (iter->index != iter->end)
1917         {
1918                 iter->index += iter->stride;
1919                 result = &(iter->arc->buckets[iter->index]);
1920         }
1921         
1922         return result;
1923 }