made most variables which are only used in a single file and not defined in header...
[blender.git] / source / blender / editors / armature / 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, 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 #include <float.h>
33
34 #include "DNA_scene_types.h"
35 #include "DNA_object_types.h"
36
37 #include "MEM_guardedalloc.h"
38
39 #include "BKE_context.h"
40
41 #include "BLI_blenlib.h"
42 #include "BLI_math.h"
43 #include "BLI_utildefines.h"
44 #include "BLI_editVert.h"
45 #include "BLI_edgehash.h"
46 #include "BLI_ghash.h"
47 #include "BLI_heap.h"
48
49 //#include "BDR_editobject.h"
50
51 //#include "BIF_interface.h"
52 //#include "BIF_toolbox.h"
53 //#include "BIF_graphics.h"
54
55
56 //#include "blendef.h"
57
58 #include "ONL_opennl.h"
59
60 #include "reeb.h"
61
62
63 static ReebGraph *GLOBAL_RG = NULL;
64 static ReebGraph *FILTERED_RG = NULL;
65
66 /*
67  * Skeleton generation algorithm based on: 
68  * "Harmonic Skeleton for Realistic Character Animation"
69  * Gregoire Aujay, Franck Hetroy, Francis Lazarus and Christine Depraz
70  * SIGGRAPH 2007
71  * 
72  * Reeb graph generation algorithm based on: 
73  * "Robust On-line Computation of Reeb Graphs: Simplicity and Speed"
74  * Valerio Pascucci, Giorgio Scorzelli, Peer-Timo Bremer and Ajith Mascarenhas
75  * SIGGRAPH 2007
76  * 
77  * */
78  
79 #define DEBUG_REEB
80 #define DEBUG_REEB_NODE
81
82 typedef struct VertexData
83 {
84         float w; /* weight */
85         int i; /* index */
86         ReebNode *n;
87 } VertexData;
88
89 typedef struct EdgeIndex
90 {
91         EditEdge **edges;
92         int              *offset;
93 } EdgeIndex;
94
95 typedef enum {
96         MERGE_LOWER,
97         MERGE_HIGHER,
98         MERGE_APPEND
99 } MergeDirection;
100
101 int mergeArcs(ReebGraph *rg, ReebArc *a0, ReebArc *a1);
102 void mergeArcEdges(ReebGraph *rg, ReebArc *aDst, ReebArc *aSrc, MergeDirection direction);
103 int mergeConnectedArcs(ReebGraph *rg, ReebArc *a0, ReebArc *a1);
104 EditEdge * NextEdgeForVert(EdgeIndex *indexed_edges, int index);
105 void mergeArcFaces(ReebGraph *rg, ReebArc *aDst, ReebArc *aSrc);
106 void addFacetoArc(ReebArc *arc, EditFace *efa);
107
108 void REEB_RadialSymmetry(BNode* root_node, RadialArc* ring, int count);
109 void REEB_AxialSymmetry(BNode* root_node, BNode* node1, BNode* node2, struct BArc* barc1, BArc* barc2);
110
111 void flipArcBuckets(ReebArc *arc);
112
113
114 /***************************************** UTILS **********************************************/
115
116 static VertexData *allocVertexData(EditMesh *em)
117 {
118         VertexData *data;
119         EditVert *eve;
120         int totvert, index;
121         
122         totvert = BLI_countlist(&em->verts);
123         
124         data = MEM_callocN(sizeof(VertexData) * totvert, "VertexData");
125
126         for(index = 0, eve = em->verts.first; eve; index++, eve = eve->next)
127         {
128                 data[index].i = index;
129                 data[index].w = 0;
130                 eve->tmp.p = data + index;
131         }
132                 
133         return data;
134 }
135
136 static int indexData(EditVert *eve)
137 {
138         return ((VertexData*)eve->tmp.p)->i;
139 }
140
141 static float weightData(EditVert *eve)
142 {
143         return ((VertexData*)eve->tmp.p)->w;
144 }
145
146 static void weightSetData(EditVert *eve, float w)
147 {
148         ((VertexData*)eve->tmp.p)->w = w;
149 }
150
151 static ReebNode* nodeData(EditVert *eve)
152 {
153         return ((VertexData*)eve->tmp.p)->n;
154 }
155
156 static void nodeSetData(EditVert *eve, ReebNode *n)
157 {
158         ((VertexData*)eve->tmp.p)->n = n;
159 }
160
161 void REEB_freeArc(BArc *barc)
162 {
163         ReebArc *arc = (ReebArc*)barc;
164         BLI_freelistN(&arc->edges);
165         
166         if (arc->buckets)
167                 MEM_freeN(arc->buckets);
168                 
169         if (arc->faces)
170                 BLI_ghash_free(arc->faces, NULL, NULL);
171         
172         MEM_freeN(arc);
173 }
174
175 void REEB_freeGraph(ReebGraph *rg)
176 {
177         ReebArc *arc;
178         ReebNode *node;
179         
180         // free nodes
181         for( node = rg->nodes.first; node; node = node->next )
182         {
183                 BLI_freeNode((BGraph*)rg, (BNode*)node);
184         }
185         BLI_freelistN(&rg->nodes);
186         
187         // free arcs
188         arc = rg->arcs.first;
189         while( arc )
190         {
191                 ReebArc *next = arc->next;
192                 REEB_freeArc((BArc*)arc);
193                 arc = next;
194         }
195         
196         // free edge map
197         BLI_edgehash_free(rg->emap, NULL);
198         
199         /* free linked graph */
200         if (rg->link_up)
201         {
202                 REEB_freeGraph(rg->link_up);
203         }
204         
205         MEM_freeN(rg);
206 }
207
208 ReebGraph * newReebGraph(void)
209 {
210         ReebGraph *rg;
211         rg = MEM_callocN(sizeof(ReebGraph), "reeb graph");
212         
213         rg->totnodes = 0;
214         rg->emap = BLI_edgehash_new();
215         
216         
217         rg->free_arc = REEB_freeArc;
218         rg->free_node = NULL;
219         rg->radial_symmetry = REEB_RadialSymmetry;
220         rg->axial_symmetry = REEB_AxialSymmetry;
221         
222         return rg;
223 }
224
225 void BIF_flagMultiArcs(ReebGraph *rg, int flag)
226 {
227         for ( ; rg; rg = rg->link_up)
228         {
229                 BLI_flagArcs((BGraph*)rg, flag);
230         }
231 }
232
233 static ReebNode * addNode(ReebGraph *rg, EditVert *eve)
234 {
235         float weight;
236         ReebNode *node = NULL;
237         
238         weight = weightData(eve);
239         
240         node = MEM_callocN(sizeof(ReebNode), "reeb node");
241         
242         node->flag = 0; // clear flag on init
243         node->symmetry_level = 0;
244         node->arcs = NULL;
245         node->degree = 0;
246         node->weight = weight;
247         node->index = rg->totnodes;
248         VECCOPY(node->p, eve->co);      
249         
250         BLI_addtail(&rg->nodes, node);
251         rg->totnodes++;
252         
253         nodeSetData(eve, node);
254         
255         return node;
256 }
257
258 static ReebNode * copyNode(ReebGraph *rg, ReebNode *node)
259 {
260         ReebNode *cp_node = NULL;
261         
262         cp_node = MEM_callocN(sizeof(ReebNode), "reeb node copy");
263         
264         memcpy(cp_node, node, sizeof(ReebNode));
265         
266         cp_node->prev = NULL;
267         cp_node->next = NULL;
268         cp_node->arcs = NULL;
269         
270         cp_node->link_up = NULL;
271         cp_node->link_down = NULL;
272         
273         BLI_addtail(&rg->nodes, cp_node);
274         rg->totnodes++;
275         
276         return cp_node; 
277 }
278
279 static void relinkNodes(ReebGraph *low_rg, ReebGraph *high_rg)
280 {
281         ReebNode *low_node, *high_node;
282         
283         if (low_rg == NULL || high_rg == NULL)
284         {
285                 return;
286         }
287         
288         for (low_node = low_rg->nodes.first; low_node; low_node = low_node->next)
289         {
290                 for (high_node = high_rg->nodes.first; high_node; high_node = high_node->next)
291                 {
292                         if (low_node->index == high_node->index)
293                         {
294                                 high_node->link_down = low_node;
295                                 low_node->link_up = high_node;
296                                 break;
297                         }
298                 }
299         }
300 }
301
302 ReebNode *BIF_otherNodeFromIndex(ReebArc *arc, ReebNode *node)
303 {
304         return (arc->head->index == node->index) ? arc->tail : arc->head;
305 }
306
307 ReebNode *BIF_NodeFromIndex(ReebArc *arc, ReebNode *node)
308 {
309         return (arc->head->index == node->index) ? arc->head : arc->tail;
310 }
311
312 ReebNode *BIF_lowestLevelNode(ReebNode *node)
313 {
314         while (node->link_down)
315         {
316                 node = node->link_down;
317         }
318         
319         return node;
320 }
321
322 static ReebArc * copyArc(ReebGraph *rg, ReebArc *arc)
323 {
324         ReebArc *cp_arc;
325         ReebNode *node;
326         
327         cp_arc = MEM_callocN(sizeof(ReebArc), "reeb arc copy");
328
329         memcpy(cp_arc, arc, sizeof(ReebArc));
330         
331         cp_arc->link_up = arc;
332         
333         cp_arc->head = NULL;
334         cp_arc->tail = NULL;
335
336         cp_arc->prev = NULL;
337         cp_arc->next = NULL;
338
339         cp_arc->edges.first = NULL;
340         cp_arc->edges.last = NULL;
341
342         /* copy buckets */      
343         cp_arc->buckets = MEM_callocN(sizeof(EmbedBucket) * cp_arc->bcount, "embed bucket");
344         memcpy(cp_arc->buckets, arc->buckets, sizeof(EmbedBucket) * cp_arc->bcount);
345         
346         /* copy faces map */
347         cp_arc->faces = BLI_ghash_new(BLI_ghashutil_ptrhash, BLI_ghashutil_ptrcmp, "copyArc gh");
348         mergeArcFaces(rg, cp_arc, arc);
349         
350         /* find corresponding head and tail */
351         for (node = rg->nodes.first; node && (cp_arc->head == NULL || cp_arc->tail == NULL); node = node->next)
352         {
353                 if (node->index == arc->head->index)
354                 {
355                         cp_arc->head = node;
356                 }
357                 else if (node->index == arc->tail->index)
358                 {
359                         cp_arc->tail = node;
360                 }
361         }
362         
363         BLI_addtail(&rg->arcs, cp_arc);
364         
365         return cp_arc;
366 }
367
368 static ReebGraph * copyReebGraph(ReebGraph *rg, int level)
369 {
370         ReebNode *node;
371         ReebArc *arc;
372         ReebGraph *cp_rg = newReebGraph();
373         
374         cp_rg->resolution = rg->resolution;
375         cp_rg->length = rg->length;
376         cp_rg->link_up = rg;
377         cp_rg->multi_level = level;
378
379         /* Copy nodes */        
380         for (node = rg->nodes.first; node; node = node->next)
381         {
382                 ReebNode *cp_node = copyNode(cp_rg, node);
383                 cp_node->multi_level = level;
384         }
385         
386         /* Copy arcs */
387         for (arc = rg->arcs.first; arc; arc = arc->next)
388         {
389                 copyArc(cp_rg, arc);
390         }
391         
392         BLI_buildAdjacencyList((BGraph*)cp_rg);
393         
394         return cp_rg;
395 }
396
397 ReebGraph *BIF_graphForMultiNode(ReebGraph *rg, ReebNode *node)
398 {
399         ReebGraph *multi_rg = rg;
400         
401         while(multi_rg && multi_rg->multi_level != node->multi_level)
402         {
403                 multi_rg = multi_rg->link_up;
404         }
405         
406         return multi_rg;
407 }
408
409 static ReebEdge * copyEdge(ReebEdge *edge)
410 {
411         ReebEdge *newEdge = NULL;
412         
413         newEdge = MEM_callocN(sizeof(ReebEdge), "reeb edge");
414         memcpy(newEdge, edge, sizeof(ReebEdge));
415         
416         newEdge->next = NULL;
417         newEdge->prev = NULL;
418         
419         return newEdge;
420 }
421
422 static void printArc(ReebArc *arc)
423 {
424         ReebEdge *edge;
425         ReebNode *head = (ReebNode*)arc->head;
426         ReebNode *tail = (ReebNode*)arc->tail;
427         printf("arc: (%i) %f -> (%i) %f\n", head->index, head->weight, tail->index, tail->weight);
428         
429         for(edge = arc->edges.first; edge ; edge = edge->next)
430         {
431                 printf("\tedge (%i, %i)\n", edge->v1->index, edge->v2->index);
432         }
433 }
434
435 static void flipArc(ReebArc *arc)
436 {
437         ReebNode *tmp;
438         tmp = arc->head;
439         arc->head = arc->tail;
440         arc->tail = tmp;
441         
442         flipArcBuckets(arc);
443 }
444
445 #ifdef DEBUG_REEB_NODE
446 static void NodeDegreeDecrement(ReebGraph *UNUSED(rg), ReebNode *node)
447 {
448         node->degree--;
449
450 //      if (node->degree == 0)
451 //      {
452 //              printf("would remove node %i\n", node->index);
453 //      }
454 }
455
456 static void NodeDegreeIncrement(ReebGraph *UNUSED(rg), ReebNode *node)
457 {
458 //      if (node->degree == 0)
459 //      {
460 //              printf("first connect node %i\n", node->index);
461 //      }
462
463         node->degree++;
464 }
465
466 #else
467 #define NodeDegreeDecrement(rg, node) {node->degree--;}
468 #define NodeDegreeIncrement(rg, node) {node->degree++;}
469 #endif
470
471 void repositionNodes(ReebGraph *rg)
472 {
473         BArc *arc = NULL;
474         BNode *node = NULL;
475         
476         // Reset node positions
477         for(node = rg->nodes.first; node; node = node->next)
478         {
479                 node->p[0] = node->p[1] = node->p[2] = 0;
480         }
481         
482         for(arc = rg->arcs.first; arc; arc = arc->next)
483         {
484                 if (((ReebArc*)arc)->bcount > 0)
485                 {
486                         float p[3];
487                         
488                         VECCOPY(p, ((ReebArc*)arc)->buckets[0].p);
489                         mul_v3_fl(p, 1.0f / arc->head->degree);
490                         add_v3_v3(arc->head->p, p);
491                         
492                         VECCOPY(p, ((ReebArc*)arc)->buckets[((ReebArc*)arc)->bcount - 1].p);
493                         mul_v3_fl(p, 1.0f / arc->tail->degree);
494                         add_v3_v3(arc->tail->p, p);
495                 }
496         }
497 }
498
499 static void verifyNodeDegree(ReebGraph *rg)
500 {
501 #ifdef DEBUG_REEB
502         ReebNode *node = NULL;
503         ReebArc *arc = NULL;
504
505         for(node = rg->nodes.first; node; node = node->next)
506         {
507                 int count = 0;
508                 for(arc = rg->arcs.first; arc; arc = arc->next)
509                 {
510                         if (arc->head == node || arc->tail == node)
511                         {
512                                 count++;
513                         }
514                 }
515                 if (count != node->degree)
516                 {
517                         printf("degree error in node %i: expected %i got %i\n", node->index, count, node->degree);
518                 }
519                 if (node->degree == 0)
520                 {
521                         printf("zero degree node %i with weight %f\n", node->index, node->weight);
522                 }
523         }
524 #endif
525 }
526
527 static void verifyBucketsArc(ReebGraph *UNUSED(rg), ReebArc *arc)
528 {
529         ReebNode *head = (ReebNode*)arc->head;
530         ReebNode *tail = (ReebNode*)arc->tail;
531
532         if (arc->bcount > 0)
533         {
534                 int i;
535                 for(i = 0; i < arc->bcount; i++)
536                 {
537                         if (arc->buckets[i].nv == 0)
538                         {
539                                 printArc(arc);
540                                 printf("count error in bucket %i/%i\n", i+1, arc->bcount);
541                         }
542                 }
543                 
544                 if (ceil(head->weight) != arc->buckets[0].val)
545                 {
546                         printArc(arc);
547                         printf("alloc error in first bucket: %f should be %f \n", arc->buckets[0].val, ceil(head->weight));
548                 }
549                 if (floor(tail->weight) != arc->buckets[arc->bcount - 1].val)
550                 {
551                         printArc(arc);
552                         printf("alloc error in last bucket: %f should be %f \n", arc->buckets[arc->bcount - 1].val, floor(tail->weight));
553                 }
554         }
555 }
556
557 void verifyBuckets(ReebGraph *rg)
558 {
559 #ifdef DEBUG_REEB
560         ReebArc *arc = NULL;
561         for(arc = rg->arcs.first; arc; arc = arc->next)
562         {
563                 verifyBucketsArc(rg, arc);
564         }
565 #endif
566 }
567
568 void verifyFaces(ReebGraph *rg)
569 {
570 #ifdef DEBUG_REEB
571         int total = 0;
572         ReebArc *arc = NULL;
573         for(arc = rg->arcs.first; arc; arc = arc->next)
574         {
575                 total += BLI_ghash_size(arc->faces);
576         }
577         
578 #endif
579 }
580
581 static void verifyArcs(ReebGraph *rg)
582 {
583         ReebArc *arc;
584         
585         for (arc = rg->arcs.first; arc; arc = arc->next)
586         {
587                 if (arc->head->weight > arc->tail->weight)
588                 {
589                         printf("FLIPPED ARC!\n");
590                 }
591         }
592 }
593
594 static void verifyMultiResolutionLinks(ReebGraph *rg, int level)
595 {
596 #ifdef DEBUG_REEB
597         ReebGraph *lower_rg = rg->link_up;
598         
599         if (lower_rg)
600         {
601                 ReebArc *arc;
602                 
603                 for (arc = rg->arcs.first; arc; arc = arc->next)
604                 {
605                         if (BLI_findindex(&lower_rg->arcs, arc->link_up) == -1)
606                         {
607                                 printf("missing arc %p for level %i\n", (void *)arc->link_up, level);
608                                 printf("Source arc was ---\n");
609                                 printArc(arc);
610
611                                 arc->link_up = NULL;
612                         }
613                 }
614                 
615                 
616                 verifyMultiResolutionLinks(lower_rg, level + 1);
617         }
618 #endif
619 }
620 /***************************************** BUCKET UTILS **********************************************/
621
622 static void addVertToBucket(EmbedBucket *b, float co[3])
623 {
624         b->nv++;
625         interp_v3_v3v3(b->p, b->p, co, 1.0f / b->nv);
626 }
627
628 static void removeVertFromBucket(EmbedBucket *b, float co[3])
629 {
630         mul_v3_fl(b->p, (float)b->nv);
631         sub_v3_v3(b->p, co);
632         b->nv--;
633         mul_v3_fl(b->p, 1.0f / (float)b->nv);
634 }
635
636 static void mergeBuckets(EmbedBucket *bDst, EmbedBucket *bSrc)
637 {
638         if (bDst->nv > 0 && bSrc->nv > 0)
639         {
640                 bDst->nv += bSrc->nv;
641                 interp_v3_v3v3(bDst->p, bDst->p, bSrc->p, (float)bSrc->nv / (float)(bDst->nv));
642         }
643         else if (bSrc->nv > 0)
644         {
645                 bDst->nv = bSrc->nv;
646                 VECCOPY(bDst->p, bSrc->p);
647         }
648 }
649
650 static void mergeArcBuckets(ReebArc *aDst, ReebArc *aSrc, float start, float end)
651 {
652         if (aDst->bcount > 0 && aSrc->bcount > 0)
653         {
654                 int indexDst = 0, indexSrc = 0;
655                 
656                 start = MAX3(start, aDst->buckets[0].val, aSrc->buckets[0].val);
657                 
658                 while(indexDst < aDst->bcount && aDst->buckets[indexDst].val < start)
659                 {
660                         indexDst++;
661                 }
662
663                 while(indexSrc < aSrc->bcount && aSrc->buckets[indexSrc].val < start)
664                 {
665                         indexSrc++;
666                 }
667                 
668                 for( ;  indexDst < aDst->bcount &&
669                                 indexSrc < aSrc->bcount &&
670                                 aDst->buckets[indexDst].val <= end &&
671                                 aSrc->buckets[indexSrc].val <= end
672                                 
673                          ;      indexDst++, indexSrc++)
674                 {
675                         mergeBuckets(aDst->buckets + indexDst, aSrc->buckets + indexSrc);
676                 }
677         }
678 }
679
680 void flipArcBuckets(ReebArc *arc)
681 {
682         int i, j;
683         
684         for (i = 0, j = arc->bcount - 1; i < j; i++, j--)
685         {
686                 EmbedBucket tmp;
687                 
688                 tmp = arc->buckets[i];
689                 arc->buckets[i] = arc->buckets[j];
690                 arc->buckets[j] = tmp;
691         }
692 }
693
694 static int countArcBuckets(ReebArc *arc)
695 {
696         return (int)(floor(arc->tail->weight) - ceil(arc->head->weight)) + 1;
697 }
698
699 static void allocArcBuckets(ReebArc *arc)
700 {
701         int i;
702         float start = ceil(arc->head->weight);
703         arc->bcount = countArcBuckets(arc);
704         
705         if (arc->bcount > 0)
706         {
707                 arc->buckets = MEM_callocN(sizeof(EmbedBucket) * arc->bcount, "embed bucket");
708                 
709                 for(i = 0; i < arc->bcount; i++)
710                 {
711                         arc->buckets[i].val = start + i;
712                 }
713         }
714         else
715         {
716                 arc->buckets = NULL;
717         }
718         
719 }
720
721 static void resizeArcBuckets(ReebArc *arc)
722 {
723         EmbedBucket *oldBuckets = arc->buckets;
724         int oldBCount = arc->bcount;
725         
726         if (countArcBuckets(arc) == oldBCount)
727         {
728                 return;
729         }
730         
731         allocArcBuckets(arc);
732         
733         if (oldBCount != 0 && arc->bcount != 0)
734         {
735                 int oldStart = (int)oldBuckets[0].val;
736                 int oldEnd = (int)oldBuckets[oldBCount - 1].val;
737                 int newStart = (int)arc->buckets[0].val;
738                 int newEnd = (int)arc->buckets[arc->bcount - 1].val;
739                 int oldOffset = 0;
740                 int newOffset = 0;
741                 int len;
742                 
743                 if (oldStart < newStart)
744                 {
745                         oldOffset = newStart - oldStart;
746                 }
747                 else
748                 {
749                         newOffset = oldStart - newStart;
750                 }
751                 
752                 len = MIN2(oldEnd - (oldStart + oldOffset) + 1, newEnd - (newStart - newOffset) + 1);
753                 
754                 memcpy(arc->buckets + newOffset, oldBuckets + oldOffset, len * sizeof(EmbedBucket)); 
755         }
756
757         if (oldBuckets != NULL)
758         {
759                 MEM_freeN(oldBuckets);
760         }
761 }
762
763 static void reweightBuckets(ReebArc *arc)
764 {
765         int i;
766         float start = ceil((arc->head)->weight);
767         
768         if (arc->bcount > 0)
769         {
770                 for(i = 0; i < arc->bcount; i++)
771                 {
772                         arc->buckets[i].val = start + i;
773                 }
774         }
775 }
776
777 static void interpolateBuckets(ReebArc *arc, float *start_p, float *end_p, int start_index, int end_index)
778 {
779         int total;
780         int j;
781         
782         total = end_index - start_index + 2;
783         
784         for (j = start_index; j <= end_index; j++)
785         {
786                 EmbedBucket *empty = arc->buckets + j;
787                 empty->nv = 1;
788                 interp_v3_v3v3(empty->p, start_p, end_p, (float)(j - start_index + 1) / total);
789         }
790 }
791
792 static void fillArcEmptyBuckets(ReebArc *arc)
793 {
794         float *start_p, *end_p;
795         int start_index = 0, end_index = 0;
796         int missing = 0;
797         int i;
798         
799         start_p = arc->head->p;
800         
801         for(i = 0; i < arc->bcount; i++)
802         {
803                 EmbedBucket *bucket = arc->buckets + i;
804                 
805                 if (missing)
806                 {
807                         if (bucket->nv > 0)
808                         {
809                                 missing = 0;
810                                 
811                                 end_p = bucket->p;
812                                 end_index = i - 1;
813                                 
814                                 interpolateBuckets(arc, start_p, end_p, start_index, end_index);
815                         }
816                 }
817                 else
818                 {
819                         if (bucket->nv == 0)
820                         {
821                                 missing = 1;
822                                 
823                                 if (i > 0)
824                                 {
825                                         start_p = arc->buckets[i - 1].p;
826                                 }
827                                 start_index = i;
828                         }
829                 }
830         }
831         
832         if (missing)
833         {
834                 end_p = arc->tail->p;
835                 end_index = arc->bcount - 1;
836                 
837                 interpolateBuckets(arc, start_p, end_p, start_index, end_index);
838         }
839 }
840
841 static void ExtendArcBuckets(ReebArc *arc)
842 {
843         ReebArcIterator arc_iter;
844         BArcIterator *iter = (BArcIterator*)&arc_iter;
845         EmbedBucket *last_bucket, *first_bucket;
846         float *previous = NULL;
847         float average_length = 0, length;
848         int padding_head = 0, padding_tail = 0;
849         
850         if (arc->bcount == 0)
851         {
852                 return; /* failsafe, shouldn't happen */
853         }
854         
855         initArcIterator(iter, arc, arc->head);
856         IT_next(iter);
857         previous = iter->p;
858         
859         for (   IT_next(iter);
860                         IT_stopped(iter) == 0;
861                         previous = iter->p, IT_next(iter)
862                 )
863         {
864                 average_length += len_v3v3(previous, iter->p);
865         }
866         average_length /= (arc->bcount - 1);
867         
868         first_bucket = arc->buckets;
869         last_bucket = arc->buckets + (arc->bcount - 1);
870         
871         length = len_v3v3(first_bucket->p, arc->head->p);
872         if (length > 2 * average_length)
873         {
874                 padding_head = (int)floor(length / average_length);
875         }
876
877         length = len_v3v3(last_bucket->p, arc->tail->p);
878         if (length > 2 * average_length)
879         {
880                 padding_tail = (int)floor(length / average_length);
881         }
882         
883         if (padding_head + padding_tail > 0)
884         {
885                 EmbedBucket *old_buckets = arc->buckets;
886                 
887                 arc->buckets = MEM_callocN(sizeof(EmbedBucket) * (padding_head + arc->bcount + padding_tail), "embed bucket");
888                 memcpy(arc->buckets + padding_head, old_buckets, arc->bcount * sizeof(EmbedBucket));
889                 
890                 arc->bcount = padding_head + arc->bcount + padding_tail;
891                 
892                 MEM_freeN(old_buckets);
893         }
894         
895         if (padding_head > 0)
896         {
897                 interpolateBuckets(arc, arc->head->p, first_bucket->p, 0, padding_head);
898         }
899         
900         if (padding_tail > 0)
901         {
902                 interpolateBuckets(arc, last_bucket->p, arc->tail->p, arc->bcount - padding_tail, arc->bcount - 1);
903         }
904 }
905
906 /* CALL THIS ONLY AFTER FILTERING, SINCE IT MESSES UP WEIGHT DISTRIBUTION */
907 static void extendGraphBuckets(ReebGraph *rg)
908 {
909         ReebArc *arc;
910         
911         for (arc = rg->arcs.first; arc; arc = arc->next)
912         {
913                 ExtendArcBuckets(arc);
914         }
915 }
916
917 /**************************************** LENGTH CALCULATIONS ****************************************/
918
919 static void calculateArcLength(ReebArc *arc)
920 {
921         ReebArcIterator arc_iter;
922         BArcIterator *iter = (BArcIterator*)&arc_iter;
923         float *vec0, *vec1;
924
925         arc->length = 0;
926         
927         initArcIterator(iter, arc, arc->head);
928
929         vec0 = arc->head->p;
930         vec1 = arc->head->p; /* in case there's no embedding */
931
932         while (IT_next(iter))   
933         {
934                 vec1 = iter->p;
935                 
936                 arc->length += len_v3v3(vec0, vec1);
937                 
938                 vec0 = vec1;
939         }
940         
941         arc->length += len_v3v3(arc->tail->p, vec1);    
942 }
943
944 static void calculateGraphLength(ReebGraph *rg)
945 {
946         ReebArc *arc;
947         
948         for (arc = rg->arcs.first; arc; arc = arc->next)
949         {
950                 calculateArcLength(arc);
951         }
952 }
953
954 /**************************************** SYMMETRY HANDLING ******************************************/
955
956 void REEB_RadialSymmetry(BNode* root_node, RadialArc* ring, int count)
957 {
958         ReebNode *node = (ReebNode*)root_node;
959         float axis[3];
960         int i;
961         
962         VECCOPY(axis, root_node->symmetry_axis);
963         
964         /* first pass, merge incrementally */
965         for (i = 0; i < count - 1; i++)
966         {
967                 ReebNode *node1, *node2;
968                 ReebArc *arc1, *arc2;
969                 float tangent[3];
970                 float normal[3];
971                 int j = i + 1;
972
973                 add_v3_v3v3(tangent, ring[i].n, ring[j].n);
974                 cross_v3_v3v3(normal, tangent, axis);
975                 
976                 node1 = (ReebNode*)BLI_otherNode(ring[i].arc, root_node);
977                 node2 = (ReebNode*)BLI_otherNode(ring[j].arc, root_node);
978                 
979                 arc1 = (ReebArc*)ring[i].arc;
980                 arc2 = (ReebArc*)ring[j].arc;
981
982                 /* mirror first node and mix with the second */
983                 BLI_mirrorAlongAxis(node1->p, root_node->p, normal);
984                 interp_v3_v3v3(node2->p, node2->p, node1->p, 1.0f / (j + 1));
985                 
986                 /* Merge buckets
987                  * there shouldn't be any null arcs here, but just to be safe 
988                  * */
989                 if (arc1->bcount > 0 && arc2->bcount > 0)
990                 {
991                         ReebArcIterator arc_iter1, arc_iter2;
992                         BArcIterator *iter1 = (BArcIterator*)&arc_iter1;
993                         BArcIterator *iter2 = (BArcIterator*)&arc_iter2;
994                         EmbedBucket *bucket1 = NULL, *bucket2 = NULL;
995                         
996                         initArcIterator(iter1, arc1, (ReebNode*)root_node);
997                         initArcIterator(iter2, arc2, (ReebNode*)root_node);
998                         
999                         bucket1 = IT_next(iter1);
1000                         bucket2 = IT_next(iter2);
1001                 
1002                         /* Make sure they both start at the same value */       
1003                         while(bucket1 && bucket2 && bucket1->val < bucket2->val)
1004                         {
1005                                 bucket1 = IT_next(iter1);
1006                         }
1007                         
1008                         while(bucket1 && bucket2 && bucket2->val < bucket1->val)
1009                         {
1010                                 bucket2 = IT_next(iter2);
1011                         }
1012         
1013         
1014                         for ( ;bucket1 && bucket2; bucket1 = IT_next(iter1), bucket2 = IT_next(iter2))
1015                         {
1016                                 bucket2->nv += bucket1->nv; /* add counts */
1017                                 
1018                                 /* mirror on axis */
1019                                 BLI_mirrorAlongAxis(bucket1->p, root_node->p, normal);
1020                                 /* add bucket2 in bucket1 */
1021                                 interp_v3_v3v3(bucket2->p, bucket2->p, bucket1->p, (float)bucket1->nv / (float)(bucket2->nv));
1022                         }
1023                 }
1024         }
1025         
1026         /* second pass, mirror back on previous arcs */
1027         for (i = count - 1; i > 0; i--)
1028         {
1029                 ReebNode *node1, *node2;
1030                 ReebArc *arc1, *arc2;
1031                 float tangent[3];
1032                 float normal[3];
1033                 int j = i - 1;
1034
1035                 add_v3_v3v3(tangent, ring[i].n, ring[j].n);
1036                 cross_v3_v3v3(normal, tangent, axis);
1037                 
1038                 node1 = (ReebNode*)BLI_otherNode(ring[i].arc, root_node);
1039                 node2 = (ReebNode*)BLI_otherNode(ring[j].arc, root_node);
1040                 
1041                 arc1 = (ReebArc*)ring[i].arc;
1042                 arc2 = (ReebArc*)ring[j].arc;
1043
1044                 /* copy first node than mirror */
1045                 VECCOPY(node2->p, node1->p);
1046                 BLI_mirrorAlongAxis(node2->p, root_node->p, normal);
1047                 
1048                 /* Copy buckets
1049                  * there shouldn't be any null arcs here, but just to be safe 
1050                  * */
1051                 if (arc1->bcount > 0 && arc2->bcount > 0)
1052                 {
1053                         ReebArcIterator arc_iter1, arc_iter2;
1054                         BArcIterator *iter1 = (BArcIterator*)&arc_iter1;
1055                         BArcIterator *iter2 = (BArcIterator*)&arc_iter2;
1056                         EmbedBucket *bucket1 = NULL, *bucket2 = NULL;
1057                         
1058                         initArcIterator(iter1, arc1, node);
1059                         initArcIterator(iter2, arc2, node);
1060                         
1061                         bucket1 = IT_next(iter1);
1062                         bucket2 = IT_next(iter2);
1063                 
1064                         /* Make sure they both start at the same value */       
1065                         while(bucket1 && bucket1->val < bucket2->val)
1066                         {
1067                                 bucket1 = IT_next(iter1);
1068                         }
1069                         
1070                         while(bucket2 && bucket2->val < bucket1->val)
1071                         {
1072                                 bucket2 = IT_next(iter2);
1073                         }
1074         
1075         
1076                         for ( ;bucket1 && bucket2; bucket1 = IT_next(iter1), bucket2 = IT_next(iter2))
1077                         {
1078                                 /* copy and mirror back to bucket2 */                   
1079                                 bucket2->nv = bucket1->nv;
1080                                 VECCOPY(bucket2->p, bucket1->p);
1081                                 BLI_mirrorAlongAxis(bucket2->p, node->p, normal);
1082                         }
1083                 }
1084         }
1085 }
1086
1087 void REEB_AxialSymmetry(BNode* root_node, BNode* node1, BNode* node2, struct BArc* barc1, BArc* barc2)
1088 {
1089         ReebArc *arc1, *arc2;
1090         float nor[3], p[3];
1091
1092         arc1 = (ReebArc*)barc1;
1093         arc2 = (ReebArc*)barc2;
1094
1095         VECCOPY(nor, root_node->symmetry_axis);
1096         
1097         /* mirror node2 along axis */
1098         VECCOPY(p, node2->p);
1099         BLI_mirrorAlongAxis(p, root_node->p, nor);
1100
1101         /* average with node1 */
1102         add_v3_v3(node1->p, p);
1103         mul_v3_fl(node1->p, 0.5f);
1104         
1105         /* mirror back on node2 */
1106         VECCOPY(node2->p, node1->p);
1107         BLI_mirrorAlongAxis(node2->p, root_node->p, nor);
1108         
1109         /* Merge buckets
1110          * there shouldn't be any null arcs here, but just to be safe 
1111          * */
1112         if (arc1->bcount > 0 && arc2->bcount > 0)
1113         {
1114                 ReebArcIterator arc_iter1, arc_iter2;
1115                 BArcIterator *iter1 = (BArcIterator*)&arc_iter1;
1116                 BArcIterator *iter2 = (BArcIterator*)&arc_iter2;
1117                 EmbedBucket *bucket1 = NULL, *bucket2 = NULL;
1118                 
1119                 initArcIterator(iter1, arc1, (ReebNode*)root_node);
1120                 initArcIterator(iter2, arc2, (ReebNode*)root_node);
1121                 
1122                 bucket1 = IT_next(iter1);
1123                 bucket2 = IT_next(iter2);
1124         
1125                 /* Make sure they both start at the same value */       
1126                 while(bucket1 && bucket1->val < bucket2->val)
1127                 {
1128                         bucket1 = IT_next(iter1);
1129                 }
1130                 
1131                 while(bucket2 && bucket2->val < bucket1->val)
1132                 {
1133                         bucket2 = IT_next(iter2);
1134                 }
1135
1136
1137                 for ( ;bucket1 && bucket2; bucket1 = IT_next(iter1), bucket2 = IT_next(iter2))
1138                 {
1139                         bucket1->nv += bucket2->nv; /* add counts */
1140                         
1141                         /* mirror on axis */
1142                         BLI_mirrorAlongAxis(bucket2->p, root_node->p, nor);
1143                         /* add bucket2 in bucket1 */
1144                         interp_v3_v3v3(bucket1->p, bucket1->p, bucket2->p, (float)bucket2->nv / (float)(bucket1->nv));
1145
1146                         /* copy and mirror back to bucket2 */                   
1147                         bucket2->nv = bucket1->nv;
1148                         VECCOPY(bucket2->p, bucket1->p);
1149                         BLI_mirrorAlongAxis(bucket2->p, root_node->p, nor);
1150                 }
1151         }
1152 }
1153
1154 /************************************** ADJACENCY LIST *************************************************/
1155
1156
1157 /****************************************** SMOOTHING **************************************************/
1158
1159 void postprocessGraph(ReebGraph *rg, char mode)
1160 {
1161         ReebArc *arc;
1162         float fac1 = 0, fac2 = 1, fac3 = 0;
1163
1164         switch(mode)
1165         {
1166         case SKGEN_AVERAGE:
1167                 fac1 = fac2 = fac3 = 1.0f / 3.0f;
1168                 break;
1169         case SKGEN_SMOOTH:
1170                 fac1 = fac3 = 0.25f;
1171                 fac2 = 0.5f;
1172                 break;
1173         case SKGEN_SHARPEN:
1174                 fac1 = fac2 = -0.25f;
1175                 fac2 = 1.5f;
1176                 break;
1177         default:
1178 //              XXX
1179 //              error("Unknown post processing mode");
1180                 return;
1181         }
1182         
1183         for(arc = rg->arcs.first; arc; arc = arc->next)
1184         {
1185                 EmbedBucket *buckets = arc->buckets;
1186                 int bcount = arc->bcount;
1187                 int index;
1188
1189                 for(index = 1; index < bcount - 1; index++)
1190                 {
1191                         interp_v3_v3v3(buckets[index].p, buckets[index].p, buckets[index - 1].p, fac1 / (fac1 + fac2));
1192                         interp_v3_v3v3(buckets[index].p, buckets[index].p, buckets[index + 1].p, fac3 / (fac1 + fac2 + fac3));
1193                 }
1194         }
1195 }
1196
1197 /********************************************SORTING****************************************************/
1198
1199 static int compareNodesWeight(void *vnode1, void *vnode2)
1200 {
1201         ReebNode *node1 = (ReebNode*)vnode1;
1202         ReebNode *node2 = (ReebNode*)vnode2;
1203         
1204         if (node1->weight < node2->weight)
1205         {
1206                 return -1;
1207         }
1208         if (node1->weight > node2->weight)
1209         {
1210                 return 1;
1211         }
1212         else
1213         {
1214                 return 0;
1215         }
1216 }
1217
1218 void sortNodes(ReebGraph *rg)
1219 {
1220         BLI_sortlist(&rg->nodes, compareNodesWeight);
1221 }
1222
1223 static int compareArcsWeight(void *varc1, void *varc2)
1224 {
1225         ReebArc *arc1 = (ReebArc*)varc1;
1226         ReebArc *arc2 = (ReebArc*)varc2;
1227         ReebNode *node1 = (ReebNode*)arc1->head; 
1228         ReebNode *node2 = (ReebNode*)arc2->head; 
1229         
1230         if (node1->weight < node2->weight)
1231         {
1232                 return -1;
1233         }
1234         if (node1->weight > node2->weight)
1235         {
1236                 return 1;
1237         }
1238         else
1239         {
1240                 return 0;
1241         }
1242 }
1243
1244 void sortArcs(ReebGraph *rg)
1245 {
1246         BLI_sortlist(&rg->arcs, compareArcsWeight);
1247 }
1248 /******************************************* JOINING ***************************************************/
1249
1250 static void reweightArc(ReebGraph *rg, ReebArc *arc, ReebNode *start_node, float start_weight)
1251 {
1252         ReebNode *node;
1253         float old_weight;
1254         float end_weight = start_weight + ABS(arc->tail->weight - arc->head->weight);
1255         int i;
1256         
1257         node = (ReebNode*)BLI_otherNode((BArc*)arc, (BNode*)start_node);
1258         
1259         /* prevent backtracking */
1260         if (node->flag == 1)
1261         {
1262                 return;
1263         }
1264
1265         if (arc->tail == start_node)
1266         {
1267                 flipArc(arc);
1268         }
1269         
1270         start_node->flag = 1;
1271         
1272         for (i = 0; i < node->degree; i++)
1273         {
1274                 ReebArc *next_arc = node->arcs[i];
1275                 
1276                 reweightArc(rg, next_arc, node, end_weight);
1277         }
1278
1279         /* update only if needed */     
1280         if (arc->head->weight != start_weight || arc->tail->weight != end_weight)
1281         {
1282                 old_weight = arc->head->weight; /* backup head weight, other arcs need it intact, it will be fixed by the source arc */
1283                 
1284                 arc->head->weight = start_weight;
1285                 arc->tail->weight = end_weight;
1286                 
1287                 reweightBuckets(arc);
1288                 resizeArcBuckets(arc);
1289                 fillArcEmptyBuckets(arc);
1290                 
1291                 arc->head->weight = old_weight;
1292         }
1293
1294
1295 static void reweightSubgraph(ReebGraph *rg, ReebNode *start_node, float start_weight)
1296 {
1297         int i;
1298                 
1299         BLI_flagNodes((BGraph*)rg, 0);
1300
1301         for (i = 0; i < start_node->degree; i++)
1302         {
1303                 ReebArc *next_arc = start_node->arcs[i];
1304                 
1305                 reweightArc(rg, next_arc, start_node, start_weight);
1306         }
1307         start_node->weight = start_weight;
1308 }
1309
1310 static int joinSubgraphsEnds(ReebGraph *rg, float threshold, int nb_subgraphs)
1311 {
1312         int joined = 0;
1313         int subgraph;
1314         
1315         for (subgraph = 1; subgraph <= nb_subgraphs; subgraph++)
1316         {
1317                 ReebNode *start_node, *end_node;
1318                 ReebNode *min_node_start = NULL, *min_node_end = NULL;
1319                 float min_distance = FLT_MAX;
1320                 
1321                 for (start_node = rg->nodes.first; start_node; start_node = start_node->next)
1322                 {
1323                         if (start_node->subgraph_index == subgraph && start_node->degree == 1)
1324                         {
1325                                 
1326                                 for (end_node = rg->nodes.first; end_node; end_node = end_node->next)
1327                                 {
1328                                         if (end_node->subgraph_index != subgraph)
1329                                         {
1330                                                 float distance = len_v3v3(start_node->p, end_node->p);
1331                                                 
1332                                                 if (distance < threshold && distance < min_distance)
1333                                                 {
1334                                                         min_distance = distance;
1335                                                         min_node_end = end_node;
1336                                                         min_node_start = start_node;
1337                                                 }
1338                                         }
1339                                 }
1340                         }
1341                 }
1342                 
1343                 end_node = min_node_end;
1344                 start_node = min_node_start;
1345                 
1346                 if (end_node && start_node)
1347                 {
1348                         ReebArc *start_arc, *end_arc;
1349                         int merging = 0;
1350                         
1351                         start_arc = start_node->arcs[0];
1352                         end_arc = end_node->arcs[0];
1353                         
1354                         if (start_arc->tail == start_node)
1355                         {
1356                                 reweightSubgraph(rg, end_node, start_node->weight);
1357                                 
1358                                 start_arc->tail = end_node;
1359                                 
1360                                 merging = 1;
1361                         }
1362                         else if (start_arc->head == start_node)
1363                         {
1364                                 reweightSubgraph(rg, start_node, end_node->weight);
1365
1366                                 start_arc->head = end_node;
1367
1368                                 merging = 2;
1369                         }
1370                         
1371                         if (merging)
1372                         {
1373                                 BLI_ReflagSubgraph((BGraph*)rg, end_node->flag, subgraph);
1374                                                                         
1375                                 resizeArcBuckets(start_arc);
1376                                 fillArcEmptyBuckets(start_arc);
1377                                 
1378                                 NodeDegreeIncrement(rg, end_node);
1379                                 BLI_rebuildAdjacencyListForNode((BGraph*)rg, (BNode*)end_node);
1380                                 
1381                                 BLI_removeNode((BGraph*)rg, (BNode*)start_node);
1382                         }
1383                         
1384                         joined = 1;
1385                 }               
1386         }
1387         
1388         return joined;
1389 }
1390
1391 /* Reweight graph from smallest node, fix fliped arcs */
1392 static void fixSubgraphsOrientation(ReebGraph *rg, int nb_subgraphs)
1393 {
1394         int subgraph;
1395         
1396         for (subgraph = 1; subgraph <= nb_subgraphs; subgraph++)
1397         {
1398                 ReebNode *node;
1399                 ReebNode *start_node = NULL;
1400                 
1401                 for (node = rg->nodes.first; node; node = node->next)
1402                 {
1403                         if (node->subgraph_index == subgraph)
1404                         {
1405                                 if (start_node == NULL || node->weight < start_node->weight)
1406                                 {
1407                                         start_node = node;
1408                                 }
1409                         }
1410                 }
1411                 
1412                 if (start_node)
1413                 {
1414                         reweightSubgraph(rg, start_node, start_node->weight);
1415                 }
1416         }
1417 }
1418
1419 static int joinSubgraphs(ReebGraph *rg, float threshold)
1420 {
1421         int nb_subgraphs;
1422         int joined = 0;
1423         
1424         BLI_buildAdjacencyList((BGraph*)rg);
1425         
1426         if (BLI_isGraphCyclic((BGraph*)rg))
1427         {
1428                 /* don't deal with cyclic graphs YET */
1429                 return 0;
1430         }
1431         
1432         /* sort nodes before flagging subgraphs to make sure root node is subgraph 0 */
1433         sortNodes(rg);
1434         
1435         nb_subgraphs = BLI_FlagSubgraphs((BGraph*)rg);
1436         
1437         /* Harmonic function can create flipped arcs, take the occasion to fix them */
1438 //      XXX
1439 //      if (G.scene->toolsettings->skgen_options & SKGEN_HARMONIC)
1440 //      {
1441                 fixSubgraphsOrientation(rg, nb_subgraphs);
1442 //      }
1443
1444         if (nb_subgraphs > 1)
1445         {
1446                 joined |= joinSubgraphsEnds(rg, threshold, nb_subgraphs);
1447                 
1448                 if (joined)
1449                 {
1450                         removeNormalNodes(rg);
1451                         BLI_buildAdjacencyList((BGraph*)rg);
1452                 }
1453         }
1454         
1455         return joined;
1456 }
1457
1458 /****************************************** FILTERING **************************************************/
1459
1460 static float lengthArc(ReebArc *arc)
1461 {
1462 #if 0
1463         ReebNode *head = (ReebNode*)arc->head;
1464         ReebNode *tail = (ReebNode*)arc->tail;
1465         
1466         return tail->weight - head->weight;
1467 #else
1468         return arc->length;
1469 #endif
1470 }
1471
1472 static int compareArcs(void *varc1, void *varc2)
1473 {
1474         ReebArc *arc1 = (ReebArc*)varc1;
1475         ReebArc *arc2 = (ReebArc*)varc2;
1476         float len1 = lengthArc(arc1);
1477         float len2 = lengthArc(arc2);
1478         
1479         if (len1 < len2)
1480         {
1481                 return -1;
1482         }
1483         if (len1 > len2)
1484         {
1485                 return 1;
1486         }
1487         else
1488         {
1489                 return 0;
1490         }
1491 }
1492
1493 static void filterArc(ReebGraph *rg, ReebNode *newNode, ReebNode *removedNode, ReebArc * srcArc, int merging)
1494 {
1495         ReebArc *arc = NULL, *nextArc = NULL;
1496
1497         if (merging)
1498         {
1499                 /* first pass, merge buckets for arcs that spawned the two nodes into the source arc*/
1500                 for(arc = rg->arcs.first; arc; arc = arc->next)
1501                 {
1502                         if (arc->head == srcArc->head && arc->tail == srcArc->tail && arc != srcArc)
1503                         {
1504                                 ReebNode *head = srcArc->head;
1505                                 ReebNode *tail = srcArc->tail;
1506                                 mergeArcBuckets(srcArc, arc, head->weight, tail->weight);
1507                         }
1508                 }
1509         }
1510
1511         /* second pass, replace removedNode by newNode, remove arcs that are collapsed in a loop */
1512         arc = rg->arcs.first;
1513         while(arc)
1514         {
1515                 nextArc = arc->next;
1516                 
1517                 if (arc->head == removedNode || arc->tail == removedNode)
1518                 {
1519                         if (arc->head == removedNode)
1520                         {
1521                                 arc->head = newNode;
1522                         }
1523                         else
1524                         {
1525                                 arc->tail = newNode;
1526                         }
1527
1528                         // Remove looped arcs                   
1529                         if (arc->head == arc->tail)
1530                         {
1531                                 // v1 or v2 was already newNode, since we're removing an arc, decrement degree
1532                                 NodeDegreeDecrement(rg, newNode);
1533                                 
1534                                 // If it's srcArc, it'll be removed later, so keep it for now
1535                                 if (arc != srcArc)
1536                                 {
1537                                         BLI_remlink(&rg->arcs, arc);
1538                                         REEB_freeArc((BArc*)arc);
1539                                 }
1540                         }
1541                         else
1542                         {
1543                                 /* flip arcs that flipped, can happen on diamond shapes, mostly on null arcs */
1544                                 if (arc->head->weight > arc->tail->weight)
1545                                 {
1546                                         flipArc(arc);
1547                                 }
1548                                 //newNode->degree++; // incrementing degree since we're adding an arc
1549                                 NodeDegreeIncrement(rg, newNode);
1550                                 mergeArcFaces(rg, arc, srcArc);
1551
1552                                 if (merging)
1553                                 {
1554                                         ReebNode *head = arc->head;
1555                                         ReebNode *tail = arc->tail;
1556
1557                                         // resize bucket list
1558                                         resizeArcBuckets(arc);
1559                                         mergeArcBuckets(arc, srcArc, head->weight, tail->weight);
1560                                         
1561                                         /* update length */
1562                                         arc->length += srcArc->length;
1563                                 }
1564                         }
1565                 }
1566                 
1567                 arc = nextArc;
1568         }
1569 }
1570
1571 void filterNullReebGraph(ReebGraph *rg)
1572 {
1573         ReebArc *arc = NULL, *nextArc = NULL;
1574         
1575         arc = rg->arcs.first;
1576         while(arc)
1577         {
1578                 nextArc = arc->next;
1579                 // Only collapse arcs too short to have any embed bucket
1580                 if (arc->bcount == 0)
1581                 {
1582                         ReebNode *newNode = (ReebNode*)arc->head;
1583                         ReebNode *removedNode = (ReebNode*)arc->tail;
1584                         float blend;
1585                         
1586                         blend = (float)newNode->degree / (float)(newNode->degree + removedNode->degree); // blending factors
1587                         
1588                         interp_v3_v3v3(newNode->p, removedNode->p, newNode->p, blend);
1589                         
1590                         filterArc(rg, newNode, removedNode, arc, 0);
1591
1592                         // Reset nextArc, it might have changed
1593                         nextArc = arc->next;
1594                         
1595                         BLI_remlink(&rg->arcs, arc);
1596                         REEB_freeArc((BArc*)arc);
1597                         
1598                         BLI_removeNode((BGraph*)rg, (BNode*)removedNode);
1599                 }
1600                 
1601                 arc = nextArc;
1602         }
1603 }
1604
1605 static int filterInternalExternalReebGraph(ReebGraph *rg, float threshold_internal, float threshold_external)
1606 {
1607         ReebArc *arc = NULL, *nextArc = NULL;
1608         int value = 0;
1609         
1610         BLI_sortlist(&rg->arcs, compareArcs);
1611         
1612         for (arc = rg->arcs.first; arc; arc = nextArc)
1613         {
1614                 nextArc = arc->next;
1615
1616                 // Only collapse non-terminal arcs that are shorter than threshold
1617                 if (threshold_internal > 0 && arc->head->degree > 1 && arc->tail->degree > 1 && (lengthArc(arc) < threshold_internal))
1618                 {
1619                         ReebNode *newNode = NULL;
1620                         ReebNode *removedNode = NULL;
1621                         
1622                         /* Always remove lower node, so arcs don't flip */
1623                         newNode = arc->head;
1624                         removedNode = arc->tail;
1625
1626                         filterArc(rg, newNode, removedNode, arc, 1);
1627
1628                         // Reset nextArc, it might have changed
1629                         nextArc = arc->next;
1630                         
1631                         BLI_remlink(&rg->arcs, arc);
1632                         REEB_freeArc((BArc*)arc);
1633                         
1634                         BLI_removeNode((BGraph*)rg, (BNode*)removedNode);
1635                         value = 1;
1636                 }
1637                 
1638                 // Only collapse terminal arcs that are shorter than threshold
1639                 else if (threshold_external > 0 && (arc->head->degree == 1 || arc->tail->degree == 1) && (lengthArc(arc) < threshold_external))
1640                 {
1641                         ReebNode *terminalNode = NULL;
1642                         ReebNode *middleNode = NULL;
1643                         ReebNode *removedNode = NULL;
1644                         
1645                         // Assign terminal and middle nodes
1646                         if (arc->head->degree == 1)
1647                         {
1648                                 terminalNode = arc->head;
1649                                 middleNode = arc->tail;
1650                         }
1651                         else
1652                         {
1653                                 terminalNode = arc->tail;
1654                                 middleNode = arc->head;
1655                         }
1656                         
1657                         if (middleNode->degree == 2 && middleNode != rg->nodes.first)
1658                         {
1659 #if 1
1660                                 // If middle node is a normal node, it will be removed later
1661                                 // Only if middle node is not the root node
1662                                 /* USE THIS IF YOU WANT TO PROLONG ARCS TO THEIR TERMINAL NODES
1663                                  * FOR HANDS, THIS IS NOT THE BEST RESULT 
1664                                  * */
1665                                 continue;
1666 #else
1667                                 removedNode = terminalNode;
1668
1669                                 // removing arc, so we need to decrease the degree of the remaining node
1670                                 NodeDegreeDecrement(rg, middleNode);
1671 #endif
1672                         }
1673                         // Otherwise, just plain remove of the arc
1674                         else
1675                         {
1676                                 removedNode = terminalNode;
1677
1678                                 // removing arc, so we need to decrease the degree of the remaining node
1679                                 NodeDegreeDecrement(rg, middleNode);
1680                         }
1681
1682                         // Reset nextArc, it might have changed
1683                         nextArc = arc->next;
1684                         
1685                         BLI_remlink(&rg->arcs, arc);
1686                         REEB_freeArc((BArc*)arc);
1687                         
1688                         BLI_removeNode((BGraph*)rg, (BNode*)removedNode);
1689                         value = 1;
1690                 }
1691         }
1692         
1693         return value;
1694 }
1695
1696 static int filterCyclesReebGraph(ReebGraph *rg, float UNUSED(distance_threshold))
1697 {
1698         ReebArc *arc1, *arc2;
1699         ReebArc *next2;
1700         int filtered = 0;
1701         
1702         for (arc1 = rg->arcs.first; arc1; arc1 = arc1->next)
1703         {
1704                 for (arc2 = arc1->next; arc2; arc2 = next2)
1705                 {
1706                         next2 = arc2->next;
1707                         if (arc1 != arc2 && arc1->head == arc2->head && arc1->tail == arc2->tail)
1708                         {
1709                                 mergeArcEdges(rg, arc1, arc2, MERGE_APPEND);
1710                                 mergeArcFaces(rg, arc1, arc2);
1711                                 mergeArcBuckets(arc1, arc2, arc1->head->weight, arc1->tail->weight);
1712
1713                                 NodeDegreeDecrement(rg, arc1->head);
1714                                 NodeDegreeDecrement(rg, arc1->tail);
1715
1716                                 BLI_remlink(&rg->arcs, arc2);
1717                                 REEB_freeArc((BArc*)arc2);
1718                                 
1719                                 filtered = 1;
1720                         }
1721                 }
1722         }
1723         
1724         return filtered;
1725 }
1726
1727 int filterSmartReebGraph(ReebGraph *UNUSED(rg), float UNUSED(threshold))
1728 {
1729         int value = 0;
1730 #if 0 //XXX
1731         ReebArc *arc = NULL, *nextArc = NULL;
1732         
1733         BLI_sortlist(&rg->arcs, compareArcs);
1734
1735 #ifdef DEBUG_REEB
1736         {       
1737                 EditFace *efa;
1738                 for(efa=G.editMesh->faces.first; efa; efa=efa->next) {
1739                         efa->tmp.fp = -1;
1740                 }
1741         }
1742 #endif
1743
1744         arc = rg->arcs.first;
1745         while(arc)
1746         {
1747                 nextArc = arc->next;
1748                 
1749                 /* need correct normals and center */
1750                 recalc_editnormals();
1751
1752                 // Only test terminal arcs
1753                 if (arc->head->degree == 1 || arc->tail->degree == 1)
1754                 {
1755                         GHashIterator ghi;
1756                         int merging = 0;
1757                         int total = BLI_ghash_size(arc->faces);
1758                         float avg_angle = 0;
1759                         float avg_vec[3] = {0,0,0};
1760                         
1761                         for(BLI_ghashIterator_init(&ghi, arc->faces);
1762                                 !BLI_ghashIterator_isDone(&ghi);
1763                                 BLI_ghashIterator_step(&ghi))
1764                         {
1765                                 EditFace *efa = BLI_ghashIterator_getValue(&ghi);
1766
1767 #if 0
1768                                 ReebArcIterator arc_iter;
1769                                 BArcIterator *iter = (BArcIterator*)&arc_iter;
1770                                 EmbedBucket *bucket = NULL;
1771                                 EmbedBucket *previous = NULL;
1772                                 float min_distance = -1;
1773                                 float angle = 0;
1774                 
1775                                 initArcIterator(iter, arc, arc->head);
1776                 
1777                                 bucket = nextBucket(iter);
1778                                 
1779                                 while (bucket != NULL)
1780                                 {
1781                                         float *vec0 = NULL;
1782                                         float *vec1 = bucket->p;
1783                                         float midpoint[3], tangent[3];
1784                                         float distance;
1785                 
1786                                         /* first bucket. Previous is head */
1787                                         if (previous == NULL)
1788                                         {
1789                                                 vec0 = arc->head->p;
1790                                         }
1791                                         /* Previous is a valid bucket */
1792                                         else
1793                                         {
1794                                                 vec0 = previous->p;
1795                                         }
1796                                         
1797                                         VECCOPY(midpoint, vec1);
1798                                         
1799                                         distance = len_v3v3(midpoint, efa->cent);
1800                                         
1801                                         if (min_distance == -1 || distance < min_distance)
1802                                         {
1803                                                 min_distance = distance;
1804                                         
1805                                                 sub_v3_v3v3(tangent, vec1, vec0);
1806                                                 normalize_v3(tangent);
1807                                                 
1808                                                 angle = dot_v3v3(tangent, efa->n);
1809                                         }
1810                                         
1811                                         previous = bucket;
1812                                         bucket = nextBucket(iter);
1813                                 }
1814                                 
1815                                 avg_angle += saacos(fabs(angle));
1816 #ifdef DEBUG_REEB
1817                                 efa->tmp.fp = saacos(fabs(angle));
1818 #endif
1819 #else
1820                                 add_v3_v3(avg_vec, efa->n);             
1821 #endif
1822                         }
1823
1824
1825 #if 0                   
1826                         avg_angle /= total;
1827 #else
1828                         mul_v3_fl(avg_vec, 1.0 / total);
1829                         avg_angle = dot_v3v3(avg_vec, avg_vec);
1830 #endif
1831                         
1832                         arc->angle = avg_angle;
1833                         
1834                         if (avg_angle > threshold)
1835                                 merging = 1;
1836                         
1837                         if (merging)
1838                         {
1839                                 ReebNode *terminalNode = NULL;
1840                                 ReebNode *middleNode = NULL;
1841                                 ReebNode *newNode = NULL;
1842                                 ReebNode *removedNode = NULL;
1843                                 int merging = 0;
1844                                 
1845                                 // Assign terminal and middle nodes
1846                                 if (arc->head->degree == 1)
1847                                 {
1848                                         terminalNode = arc->head;
1849                                         middleNode = arc->tail;
1850                                 }
1851                                 else
1852                                 {
1853                                         terminalNode = arc->tail;
1854                                         middleNode = arc->head;
1855                                 }
1856                                 
1857                                 // If middle node is a normal node, merge to terminal node
1858                                 if (middleNode->degree == 2)
1859                                 {
1860                                         merging = 1;
1861                                         newNode = terminalNode;
1862                                         removedNode = middleNode;
1863                                 }
1864                                 // Otherwise, just plain remove of the arc
1865                                 else
1866                                 {
1867                                         merging = 0;
1868                                         newNode = middleNode;
1869                                         removedNode = terminalNode;
1870                                 }
1871                                 
1872                                 // Merging arc
1873                                 if (merging)
1874                                 {
1875                                         filterArc(rg, newNode, removedNode, arc, 1);
1876                                 }
1877                                 else
1878                                 {
1879                                         // removing arc, so we need to decrease the degree of the remaining node
1880                                         //newNode->degree--;
1881                                         NodeDegreeDecrement(rg, newNode);
1882                                 }
1883         
1884                                 // Reset nextArc, it might have changed
1885                                 nextArc = arc->next;
1886                                 
1887                                 BLI_remlink(&rg->arcs, arc);
1888                                 REEB_freeArc((BArc*)arc);
1889                                 
1890                                 BLI_freelinkN(&rg->nodes, removedNode);
1891                                 value = 1;
1892                         }
1893                 }
1894                 
1895                 arc = nextArc;
1896         }
1897         
1898         #endif
1899         
1900         return value;
1901 }
1902
1903 static void filterGraph(ReebGraph *rg, short options, float threshold_internal, float threshold_external)
1904 {
1905         int done = 1;
1906         
1907         calculateGraphLength(rg);
1908
1909         if ((options & SKGEN_FILTER_EXTERNAL) == 0)
1910         {
1911                 threshold_external = 0;
1912         }
1913
1914         if ((options & SKGEN_FILTER_INTERNAL) == 0)
1915         {
1916                 threshold_internal = 0;
1917         }
1918
1919         if (threshold_internal > 0 || threshold_external > 0)
1920         { 
1921                 /* filter until there's nothing more to do */
1922                 while (done == 1)
1923                 {
1924                         done = 0; /* no work done yet */
1925                         
1926                         done = filterInternalExternalReebGraph(rg, threshold_internal, threshold_external);
1927                 }
1928         }
1929
1930         if (options & SKGEN_FILTER_SMART)
1931         {
1932                 filterSmartReebGraph(rg, 0.5);
1933                 filterCyclesReebGraph(rg, 0.5);
1934         }
1935
1936         repositionNodes(rg);
1937
1938         /* Filtering might have created degree 2 nodes, so remove them */
1939         removeNormalNodes(rg);
1940 }
1941
1942 static void finalizeGraph(ReebGraph *rg, char passes, char method)
1943 {
1944         int i;
1945         
1946         BLI_buildAdjacencyList((BGraph*)rg);
1947
1948         sortNodes(rg);
1949         
1950         sortArcs(rg);
1951         
1952         for(i = 0; i <  passes; i++)
1953         {
1954                 postprocessGraph(rg, method);
1955         }
1956         
1957         extendGraphBuckets(rg);
1958 }
1959
1960 /************************************** WEIGHT SPREADING ***********************************************/
1961
1962 static int compareVerts( const void* a, const void* b )
1963 {
1964         EditVert *va = *(EditVert**)a;
1965         EditVert *vb = *(EditVert**)b;
1966         int value = 0;
1967         
1968         if (weightData(va) < weightData(vb))
1969         {
1970                 value = -1;
1971         }
1972         else if (weightData(va) > weightData(vb))
1973         {
1974                 value = 1;
1975         }
1976
1977         return value;           
1978 }
1979
1980 static void spreadWeight(EditMesh *em)
1981 {
1982         EditVert **verts, *eve;
1983         float lastWeight = 0.0f;
1984         int totvert = BLI_countlist(&em->verts);
1985         int i;
1986         int work_needed = 1;
1987         
1988         verts = MEM_callocN(sizeof(EditVert*) * totvert, "verts array");
1989         
1990         for(eve = em->verts.first, i = 0; eve; eve = eve->next, i++)
1991         {
1992                 verts[i] = eve;
1993         }
1994         
1995         while(work_needed == 1)
1996         {
1997                 work_needed = 0;
1998                 qsort(verts, totvert, sizeof(EditVert*), compareVerts);
1999                 
2000                 for(i = 0; i < totvert; i++)
2001                 {
2002                         eve = verts[i];
2003                         
2004                         if (i == 0 || (weightData(eve) - lastWeight) > FLT_EPSILON)
2005                         {
2006                                 lastWeight = weightData(eve);
2007                         }
2008                         else
2009                         {
2010                                 work_needed = 1;
2011                                 weightSetData(eve, lastWeight + FLT_EPSILON * 2);
2012                                 lastWeight = weightData(eve);
2013                         }
2014                 }
2015         }
2016         
2017         MEM_freeN(verts);
2018 }
2019
2020 /******************************************** EXPORT ***************************************************/
2021
2022 static void exportNode(FILE *f, const char *text, ReebNode *node)
2023 {
2024         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]);
2025 }
2026
2027 void REEB_exportGraph(ReebGraph *rg, int count)
2028 {
2029         ReebArc *arc;
2030         char filename[128];
2031         FILE *f;
2032         
2033         if (count == -1)
2034         {
2035                 sprintf(filename, "test.txt");
2036         }
2037         else
2038         {
2039                 sprintf(filename, "test%05i.txt", count);
2040         }
2041         f = fopen(filename, "w");
2042
2043         for(arc = rg->arcs.first; arc; arc = arc->next)
2044         {
2045                 int i;
2046                 float p[3];
2047                 
2048                 exportNode(f, "v1", arc->head);
2049                 
2050                 for(i = 0; i < arc->bcount; i++)
2051                 {
2052                         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]);
2053                 }
2054                 
2055                 add_v3_v3v3(p, arc->tail->p, arc->head->p);
2056                 mul_v3_fl(p, 0.5f);
2057                 
2058                 fprintf(f, "angle %0.3f %0.3f %0.3f %0.3f %i\n", p[0], p[1], p[2], arc->angle, BLI_ghash_size(arc->faces));
2059                 exportNode(f, "v2", arc->tail);
2060         }       
2061         
2062         fclose(f);
2063 }
2064
2065 /***************************************** MAIN ALGORITHM **********************************************/
2066
2067 /* edges alone will create zero degree nodes, use this function to remove them */
2068 static void removeZeroNodes(ReebGraph *rg)
2069 {
2070         ReebNode *node, *next_node;
2071         
2072         for (node = rg->nodes.first; node; node = next_node)
2073         {
2074                 next_node = node->next;
2075                 
2076                 if (node->degree == 0)
2077                 {
2078                         BLI_removeNode((BGraph*)rg, (BNode*)node);
2079                 }
2080         }
2081 }
2082
2083 void removeNormalNodes(ReebGraph *rg)
2084 {
2085         ReebArc *arc, *nextArc;
2086         
2087         // Merge degree 2 nodes
2088         for(arc = rg->arcs.first; arc; arc = nextArc)
2089         {
2090                 nextArc = arc->next;
2091                 
2092                 while (arc->head->degree == 2 || arc->tail->degree == 2)
2093                 {
2094                         // merge at v1
2095                         if (arc->head->degree == 2)
2096                         {
2097                                 ReebArc *connectedArc = (ReebArc*)BLI_findConnectedArc((BGraph*)rg, (BArc*)arc, (BNode*)arc->head);
2098
2099                                 /* If arcs are one after the other */
2100                                 if (arc->head == connectedArc->tail)
2101                                 {               
2102                                         /* remove furthest arc */               
2103                                         if (arc->tail->weight < connectedArc->head->weight)
2104                                         {
2105                                                 mergeConnectedArcs(rg, arc, connectedArc);
2106                                                 nextArc = arc->next;
2107                                         }
2108                                         else
2109                                         {
2110                                                 mergeConnectedArcs(rg, connectedArc, arc);
2111                                                 break; /* arc was removed, move to next */
2112                                         }
2113                                 }
2114                                 /* Otherwise, arcs are side by side */
2115                                 else
2116                                 {
2117                                         /* Don't do anything, we need to keep the lowest node, even if degree 2 */
2118                                         break;
2119                                 }
2120                         }
2121                         
2122                         // merge at v2
2123                         if (arc->tail->degree == 2)
2124                         {
2125                                 ReebArc *connectedArc = (ReebArc*)BLI_findConnectedArc((BGraph*)rg, (BArc*)arc, (BNode*)arc->tail);
2126                                 
2127                                 /* If arcs are one after the other */
2128                                 if (arc->tail == connectedArc->head)
2129                                 {                               
2130                                         /* remove furthest arc */               
2131                                         if (arc->head->weight < connectedArc->tail->weight)
2132                                         {
2133                                                 mergeConnectedArcs(rg, arc, connectedArc);
2134                                                 nextArc = arc->next;
2135                                         }
2136                                         else
2137                                         {
2138                                                 mergeConnectedArcs(rg, connectedArc, arc);
2139                                                 break; /* arc was removed, move to next */
2140                                         }
2141                                 }
2142                                 /* Otherwise, arcs are side by side */
2143                                 else
2144                                 {
2145                                         /* Don't do anything, we need to keep the lowest node, even if degree 2 */
2146                                         break;
2147                                 }
2148                         }
2149                 }
2150         }
2151         
2152 }
2153
2154 static int edgeEquals(ReebEdge *e1, ReebEdge *e2)
2155 {
2156         return (e1->v1 == e2->v1 && e1->v2 == e2->v2);
2157 }
2158
2159 static ReebArc *nextArcMappedToEdge(ReebArc *arc, ReebEdge *e)
2160 {
2161         ReebEdge *nextEdge = NULL;
2162         ReebEdge *edge = NULL;
2163         ReebArc *result = NULL;
2164
2165         /* Find the ReebEdge in the edge list */
2166         for(edge = arc->edges.first; edge && !edgeEquals(edge, e); edge = edge->next)
2167         {       }
2168         
2169         nextEdge = edge->nextEdge;
2170         
2171         if (nextEdge != NULL)
2172         {
2173                 result = nextEdge->arc;
2174         }
2175
2176         return result;
2177 }
2178
2179 void addFacetoArc(ReebArc *arc, EditFace *efa)
2180 {
2181         BLI_ghash_insert(arc->faces, efa, efa);
2182 }
2183
2184 void mergeArcFaces(ReebGraph *UNUSED(rg), ReebArc *aDst, ReebArc *aSrc)
2185 {
2186         GHashIterator ghi;
2187         
2188         for(BLI_ghashIterator_init(&ghi, aSrc->faces);
2189                 !BLI_ghashIterator_isDone(&ghi);
2190                 BLI_ghashIterator_step(&ghi))
2191         {
2192                 EditFace *efa = BLI_ghashIterator_getValue(&ghi);
2193                 BLI_ghash_insert(aDst->faces, efa, efa);
2194         }
2195
2196
2197 void mergeArcEdges(ReebGraph *rg, ReebArc *aDst, ReebArc *aSrc, MergeDirection direction)
2198 {
2199         ReebEdge *e = NULL;
2200         
2201         if (direction == MERGE_APPEND)
2202         {
2203                 for(e = aSrc->edges.first; e; e = e->next)
2204                 {
2205                         e->arc = aDst; // Edge is stolen by new arc
2206                 }
2207                 
2208                 BLI_movelisttolist(&aDst->edges , &aSrc->edges);
2209         }
2210         else
2211         {
2212                 for(e = aSrc->edges.first; e; e = e->next)
2213                 {
2214                         ReebEdge *newEdge = copyEdge(e);
2215
2216                         newEdge->arc = aDst;
2217                         
2218                         BLI_addtail(&aDst->edges, newEdge);
2219                         
2220                         if (direction == MERGE_LOWER)
2221                         {
2222                                 void **p = BLI_edgehash_lookup_p(rg->emap, e->v1->index, e->v2->index);
2223                                 
2224                                 newEdge->nextEdge = e;
2225
2226                                 // if edge was the first in the list, point the edit edge to the new reeb edge instead.                                                 
2227                                 if (*p == e)
2228                                 {
2229                                         *p = (void*)newEdge;
2230                                 }
2231                                 // otherwise, advance in the list until the predecessor is found then insert it there
2232                                 else
2233                                 {
2234                                         ReebEdge *previous = (ReebEdge*)*p;
2235                                         
2236                                         while(previous->nextEdge != e)
2237                                         {
2238                                                 previous = previous->nextEdge;
2239                                         }
2240                                         
2241                                         previous->nextEdge = newEdge;
2242                                 }
2243                         }
2244                         else
2245                         {
2246                                 newEdge->nextEdge = e->nextEdge;
2247                                 e->nextEdge = newEdge;
2248                         }
2249                 }
2250         }
2251
2252
2253 // return 1 on full merge
2254 int mergeConnectedArcs(ReebGraph *rg, ReebArc *a0, ReebArc *a1)
2255 {
2256         int result = 0;
2257         ReebNode *removedNode = NULL;
2258         
2259         a0->length += a1->length;
2260         
2261         mergeArcEdges(rg, a0, a1, MERGE_APPEND);
2262         mergeArcFaces(rg, a0, a1);
2263         
2264         // Bring a0 to the combine length of both arcs
2265         if (a0->tail == a1->head)
2266         {
2267                 removedNode = a0->tail;
2268                 a0->tail = a1->tail;
2269         }
2270         else if (a0->head == a1->tail)
2271         {
2272                 removedNode = a0->head;
2273                 a0->head = a1->head;
2274         }
2275         
2276         resizeArcBuckets(a0);
2277         // Merge a1 in a0
2278         mergeArcBuckets(a0, a1, a0->head->weight, a0->tail->weight);
2279         
2280         // remove a1 from graph
2281         BLI_remlink(&rg->arcs, a1);
2282         REEB_freeArc((BArc*)a1);
2283         
2284         BLI_removeNode((BGraph*)rg, (BNode*)removedNode);
2285         result = 1;
2286         
2287         return result;
2288 }
2289 // return 1 on full merge
2290 int mergeArcs(ReebGraph *rg, ReebArc *a0, ReebArc *a1)
2291 {
2292         int result = 0;
2293         // TRIANGLE POINTS DOWN
2294         if (a0->head->weight == a1->head->weight) // heads are the same
2295         {
2296                 if (a0->tail->weight == a1->tail->weight) // tails also the same, arcs can be totally merge together
2297                 {
2298                         mergeArcEdges(rg, a0, a1, MERGE_APPEND);
2299                         mergeArcFaces(rg, a0, a1);
2300                         
2301                         mergeArcBuckets(a0, a1, a0->head->weight, a0->tail->weight);
2302                         
2303                         // Adjust node degree
2304                         //a1->head->degree--;
2305                         NodeDegreeDecrement(rg, a1->head);
2306                         //a1->tail->degree--;
2307                         NodeDegreeDecrement(rg, a1->tail);
2308                         
2309                         // remove a1 from graph
2310                         BLI_remlink(&rg->arcs, a1);
2311                         
2312                         REEB_freeArc((BArc*)a1);
2313                         result = 1;
2314                 }
2315                 else if (a0->tail->weight > a1->tail->weight) // a1->tail->weight is in the middle
2316                 {
2317                         mergeArcEdges(rg, a1, a0, MERGE_LOWER);
2318                         mergeArcFaces(rg, a1, a0);
2319
2320                         // Adjust node degree
2321                         //a0->head->degree--;
2322                         NodeDegreeDecrement(rg, a0->head);
2323                         //a1->tail->degree++;
2324                         NodeDegreeIncrement(rg, a1->tail);
2325                         
2326                         mergeArcBuckets(a1, a0, a1->head->weight, a1->tail->weight);
2327                         a0->head = a1->tail;
2328                         resizeArcBuckets(a0);
2329                 }
2330                 else // a0>n2 is in the middle
2331                 {
2332                         mergeArcEdges(rg, a0, a1, MERGE_LOWER);
2333                         mergeArcFaces(rg, a0, a1);
2334                         
2335                         // Adjust node degree
2336                         //a1->head->degree--;
2337                         NodeDegreeDecrement(rg, a1->head);
2338                         //a0->tail->degree++;
2339                         NodeDegreeIncrement(rg, a0->tail);
2340                         
2341                         mergeArcBuckets(a0, a1, a0->head->weight, a0->tail->weight);
2342                         a1->head = a0->tail;
2343                         resizeArcBuckets(a1);
2344                 }
2345         }
2346         // TRIANGLE POINTS UP
2347         else if (a0->tail->weight == a1->tail->weight) // tails are the same
2348         {
2349                 if (a0->head->weight > a1->head->weight) // a0->head->weight is in the middle
2350                 {
2351                         mergeArcEdges(rg, a0, a1, MERGE_HIGHER);
2352                         mergeArcFaces(rg, a0, a1);
2353                         
2354                         // Adjust node degree
2355                         //a1->tail->degree--;
2356                         NodeDegreeDecrement(rg, a1->tail);
2357                         //a0->head->degree++;
2358                         NodeDegreeIncrement(rg, a0->head);
2359                         
2360                         mergeArcBuckets(a0, a1, a0->head->weight, a0->tail->weight);
2361                         a1->tail = a0->head;
2362                         resizeArcBuckets(a1);
2363                 }
2364                 else // a1->head->weight is in the middle
2365                 {
2366                         mergeArcEdges(rg, a1, a0, MERGE_HIGHER);
2367                         mergeArcFaces(rg, a1, a0);
2368
2369                         // Adjust node degree
2370                         //a0->tail->degree--;
2371                         NodeDegreeDecrement(rg, a0->tail);
2372                         //a1->head->degree++;
2373                         NodeDegreeIncrement(rg, a1->head);
2374
2375                         mergeArcBuckets(a1, a0, a1->head->weight, a1->tail->weight);
2376                         a0->tail = a1->head;
2377                         resizeArcBuckets(a0);
2378                 }
2379         }
2380         else
2381         {
2382                 // Need something here (OR NOT)
2383         }
2384         
2385         return result;
2386 }
2387
2388 static void glueByMergeSort(ReebGraph *rg, ReebArc *a0, ReebArc *a1, ReebEdge *e0, ReebEdge *e1)
2389 {
2390         int total = 0;
2391         while (total == 0 && a0 != a1 && a0 != NULL && a1 != NULL)
2392         {
2393                 total = mergeArcs(rg, a0, a1);
2394                 
2395                 if (total == 0) // if it wasn't a total merge, go forward
2396                 {
2397                         if (a0->tail->weight < a1->tail->weight)
2398                         {
2399                                 a0 = nextArcMappedToEdge(a0, e0);
2400                         }
2401                         else
2402                         {
2403                                 a1 = nextArcMappedToEdge(a1, e1);
2404                         }
2405                 }
2406         }
2407 }
2408
2409 static void mergePaths(ReebGraph *rg, ReebEdge *e0, ReebEdge *e1, ReebEdge *e2)
2410 {
2411         ReebArc *a0, *a1, *a2;
2412         a0 = e0->arc;
2413         a1 = e1->arc;
2414         a2 = e2->arc;
2415         
2416         glueByMergeSort(rg, a0, a1, e0, e1);
2417         glueByMergeSort(rg, a0, a2, e0, e2);
2418
2419
2420 static ReebEdge * createArc(ReebGraph *rg, ReebNode *node1, ReebNode *node2)
2421 {
2422         ReebEdge *edge;
2423         
2424         edge = BLI_edgehash_lookup(rg->emap, node1->index, node2->index);
2425         
2426         // Only add existing edges that haven't been added yet
2427         if (edge == NULL)
2428         {
2429                 ReebArc *arc;
2430                 ReebNode *v1, *v2;
2431                 float len, offset;
2432                 int i;
2433                 
2434                 arc = MEM_callocN(sizeof(ReebArc), "reeb arc");
2435                 edge = MEM_callocN(sizeof(ReebEdge), "reeb edge");
2436                 
2437                 arc->flag = 0; // clear flag on init
2438                 arc->symmetry_level = 0;
2439                 arc->faces = BLI_ghash_new(BLI_ghashutil_ptrhash, BLI_ghashutil_ptrcmp, "createArc gh");
2440                 
2441                 if (node1->weight <= node2->weight)
2442                 {
2443                         v1 = node1;     
2444                         v2 = node2;     
2445                 }
2446                 else
2447                 {
2448                         v1 = node2;     
2449                         v2 = node1;     
2450                 }
2451                 
2452                 arc->head = v1;
2453                 arc->tail = v2;
2454                 
2455                 // increase node degree
2456                 //v1->degree++;
2457                 NodeDegreeIncrement(rg, v1);
2458                 //v2->degree++;
2459                 NodeDegreeIncrement(rg, v2);
2460
2461                 BLI_edgehash_insert(rg->emap, node1->index, node2->index, edge);
2462                 
2463                 edge->arc = arc;
2464                 edge->nextEdge = NULL;
2465                 edge->v1 = v1;
2466                 edge->v2 = v2;
2467                 
2468                 BLI_addtail(&rg->arcs, arc);
2469                 BLI_addtail(&arc->edges, edge);
2470                 
2471                 /* adding buckets for embedding */
2472                 allocArcBuckets(arc);
2473                 
2474                 offset = arc->head->weight;
2475                 len = arc->tail->weight - arc->head->weight;
2476
2477 #if 0
2478                 /* This is the actual embedding filling described in the paper
2479                  * the problem is that it only works with really dense meshes
2480                  */
2481                 if (arc->bcount > 0)
2482                 {
2483                         addVertToBucket(&(arc->buckets[0]), arc->head->co);
2484                         addVertToBucket(&(arc->buckets[arc->bcount - 1]), arc->tail->co);
2485                 }
2486 #else
2487                 for(i = 0; i < arc->bcount; i++)
2488                 {
2489                         float co[3];
2490                         float f = (arc->buckets[i].val - offset) / len;
2491                         
2492                         interp_v3_v3v3(co, v1->p, v2->p, f);
2493                         addVertToBucket(&(arc->buckets[i]), co);
2494                 }
2495 #endif
2496
2497         }
2498         
2499         return edge;
2500 }
2501
2502 static void addTriangleToGraph(ReebGraph *rg, ReebNode * n1, ReebNode * n2, ReebNode * n3, EditFace *efa)
2503 {
2504         ReebEdge *re1, *re2, *re3;
2505         ReebEdge *e1, *e2, *e3;
2506         float len1, len2, len3;
2507         
2508         re1 = createArc(rg, n1, n2);
2509         re2 = createArc(rg, n2, n3);
2510         re3 = createArc(rg, n3, n1);
2511         
2512         addFacetoArc(re1->arc, efa);
2513         addFacetoArc(re2->arc, efa);
2514         addFacetoArc(re3->arc, efa);
2515         
2516         len1 = (float)fabs(n1->weight - n2->weight);
2517         len2 = (float)fabs(n2->weight - n3->weight);
2518         len3 = (float)fabs(n3->weight - n1->weight);
2519         
2520         /* The rest of the algorithm assumes that e1 is the longest edge */
2521         
2522         if (len1 >= len2 && len1 >= len3)
2523         {
2524                 e1 = re1;
2525                 e2 = re2;
2526                 e3 = re3;
2527         }
2528         else if (len2 >= len1 && len2 >= len3)
2529         {
2530                 e1 = re2;
2531                 e2 = re1;
2532                 e3 = re3;
2533         }
2534         else
2535         {
2536                 e1 = re3;
2537                 e2 = re2;
2538                 e3 = re1;
2539         }
2540         
2541         /* And e2 is the lowest edge
2542          * If e3 is lower than e2, swap them
2543          */
2544         if (e3->v1->weight < e2->v1->weight)
2545         {
2546                 ReebEdge *etmp = e2;
2547                 e2 = e3;
2548                 e3 = etmp;
2549         }
2550         
2551         
2552         mergePaths(rg, e1, e2, e3);
2553 }
2554
2555 ReebGraph * generateReebGraph(EditMesh *em, int subdivisions)
2556 {
2557         ReebGraph *rg;
2558         EditVert *eve;
2559         EditFace *efa;
2560         int index;
2561         /*int totvert;*/
2562         
2563 #ifdef DEBUG_REEB
2564         int totfaces;
2565         int countfaces = 0;
2566 #endif
2567         
2568         rg = newReebGraph();
2569         
2570         rg->resolution = subdivisions;
2571         
2572         /*totvert = BLI_countlist(&em->verts);*/ /*UNUSED*/
2573 #ifdef DEBUG_REEB
2574         totfaces = BLI_countlist(&em->faces);
2575 #endif
2576         
2577         renormalizeWeight(em, 1.0f);
2578         
2579         /* Spread weight to minimize errors */
2580         spreadWeight(em);
2581
2582         renormalizeWeight(em, (float)rg->resolution);
2583
2584         /* Adding vertice */
2585         for(index = 0, eve = em->verts.first; eve; eve = eve->next)
2586         {
2587                 if (eve->h == 0)
2588                 {
2589                         addNode(rg, eve);
2590                         eve->f2 = 0;
2591                         index++;
2592                 }
2593         }
2594         
2595         /* Adding face, edge per edge */
2596         for(efa = em->faces.first; efa; efa = efa->next)
2597         {
2598                 if (efa->h == 0)
2599                 {
2600                         ReebNode *n1, *n2, *n3;
2601                         
2602                         n1 = nodeData(efa->v1);
2603                         n2 = nodeData(efa->v2);
2604                         n3 = nodeData(efa->v3);
2605                         
2606                         addTriangleToGraph(rg, n1, n2, n3, efa);
2607                         
2608                         if (efa->v4)
2609                         {
2610                                 ReebNode *n4 = nodeData(efa->v4);
2611                                 addTriangleToGraph(rg, n1, n3, n4, efa);
2612                         }
2613 #ifdef DEBUG_REEB
2614                         countfaces++;
2615                         if (countfaces % 100 == 0)
2616                         {
2617                                 printf("\rface %i of %i", countfaces, totfaces);
2618                         }
2619 #endif
2620                 }
2621         }
2622         
2623         printf("\n");
2624         
2625         removeZeroNodes(rg);
2626         
2627         removeNormalNodes(rg);
2628         
2629         return rg;
2630 }
2631
2632 /***************************************** WEIGHT UTILS **********************************************/
2633
2634 void renormalizeWeight(EditMesh *em, float newmax)
2635 {
2636         EditVert *eve;
2637         float minimum, maximum, range;
2638         
2639         if (em == NULL || BLI_countlist(&em->verts) == 0)
2640                 return;
2641
2642         /* First pass, determine maximum and minimum */
2643         eve = em->verts.first;
2644         minimum = weightData(eve);
2645         maximum = minimum;
2646         for(; eve; eve = eve->next)
2647         {
2648                 maximum = MAX2(maximum, weightData(eve));
2649                 minimum = MIN2(minimum, weightData(eve));
2650         }
2651         
2652         range = maximum - minimum;
2653
2654         /* Normalize weights */
2655         for(eve = em->verts.first; eve; eve = eve->next)
2656         {
2657                 float weight = (weightData(eve) - minimum) / range * newmax;
2658                 weightSetData(eve, weight);
2659         }
2660 }
2661
2662
2663 int weightFromLoc(EditMesh *em, int axis)
2664 {
2665         EditVert *eve;
2666         
2667         if (em == NULL || BLI_countlist(&em->verts) == 0 || axis < 0 || axis > 2)
2668                 return 0;
2669
2670         /* Copy coordinate in weight */
2671         for(eve = em->verts.first; eve; eve = eve->next)
2672         {
2673                 weightSetData(eve, eve->co[axis]);
2674         }
2675
2676         return 1;
2677 }
2678
2679 static float cotan_weight(float *v1, float *v2, float *v3)
2680 {
2681         float a[3], b[3], c[3], clen;
2682
2683         sub_v3_v3v3(a, v2, v1);
2684         sub_v3_v3v3(b, v3, v1);
2685         cross_v3_v3v3(c, a, b);
2686
2687         clen = len_v3(c);
2688
2689         if (clen == 0.0f)
2690                 return 0.0f;
2691         
2692         return dot_v3v3(a, b)/clen;
2693 }
2694
2695 static void addTriangle(EditVert *v1, EditVert *v2, EditVert *v3, int e1, int e2, int e3)
2696 {
2697         /* Angle opposite e1 */
2698         float t1= cotan_weight(v1->co, v2->co, v3->co) / e2;
2699         
2700         /* Angle opposite e2 */
2701         float t2 = cotan_weight(v2->co, v3->co, v1->co) / e3;
2702
2703         /* Angle opposite e3 */
2704         float t3 = cotan_weight(v3->co, v1->co, v2->co) / e1;
2705         
2706         int i1 = indexData(v1);
2707         int i2 = indexData(v2);
2708         int i3 = indexData(v3);
2709         
2710         nlMatrixAdd(i1, i1, t2+t3);
2711         nlMatrixAdd(i2, i2, t1+t3);
2712         nlMatrixAdd(i3, i3, t1+t2);
2713
2714         nlMatrixAdd(i1, i2, -t3);
2715         nlMatrixAdd(i2, i1, -t3);
2716
2717         nlMatrixAdd(i2, i3, -t1);
2718         nlMatrixAdd(i3, i2, -t1);
2719
2720         nlMatrixAdd(i3, i1, -t2);
2721         nlMatrixAdd(i1, i3, -t2);
2722 }
2723
2724 int weightToHarmonic(EditMesh *em, EdgeIndex *indexed_edges)
2725 {
2726         NLboolean success;
2727         EditVert *eve;
2728         EditEdge *eed;
2729         EditFace *efa;
2730         int totvert = 0;
2731         int index;
2732         int rval;
2733         
2734         /* Find local extrema */
2735         for(eve = em->verts.first; eve; eve = eve->next)
2736         {
2737                 totvert++;
2738         }
2739
2740         /* Solve with openNL */
2741         
2742         nlNewContext();
2743
2744         nlSolverParameteri(NL_NB_VARIABLES, totvert);
2745
2746         nlBegin(NL_SYSTEM);
2747         
2748         /* Find local extrema */
2749         for(index = 0, eve = em->verts.first; eve; index++, eve = eve->next)
2750         {
2751                 if (eve->h == 0)
2752                 {
2753                         EditEdge *eed;
2754                         int maximum = 1;
2755                         int minimum = 1;
2756                         
2757                         NextEdgeForVert(indexed_edges, -1); /* Reset next edge */
2758                         for(eed = NextEdgeForVert(indexed_edges, index); eed && (maximum || minimum); eed = NextEdgeForVert(indexed_edges, index))
2759                         {
2760                                 EditVert *eve2;
2761                                 
2762                                 if (eed->v1 == eve)
2763                                 {
2764                                         eve2 = eed->v2;
2765                                 }
2766                                 else
2767                                 {
2768                                         eve2 = eed->v1;
2769                                 }
2770                                 
2771                                 if (eve2->h == 0)
2772                                 {
2773                                         /* Adjacent vertex is bigger, not a local maximum */
2774                                         if (weightData(eve2) > weightData(eve))
2775                                         {
2776                                                 maximum = 0;
2777                                         }
2778                                         /* Adjacent vertex is smaller, not a local minimum */
2779                                         else if (weightData(eve2) < weightData(eve))
2780                                         {
2781                                                 minimum = 0;
2782                                         }
2783                                 }
2784                         }
2785                         
2786                         if (maximum || minimum)
2787                         {
2788                                 float w = weightData(eve);
2789                                 eve->f1 = 0;
2790                                 nlSetVariable(0, index, w);
2791                                 nlLockVariable(index);
2792                         }
2793                         else
2794                         {
2795                                 eve->f1 = 1;
2796                         }
2797                 }
2798         }
2799         
2800         nlBegin(NL_MATRIX);
2801
2802         /* Zero edge weight */
2803         for(eed = em->edges.first; eed; eed = eed->next)
2804         {
2805                 eed->tmp.l = 0;
2806         }
2807         
2808         /* Add faces count to the edge weight */
2809         for(efa = em->faces.first; efa; efa = efa->next)
2810         {
2811                 if (efa->h == 0)
2812                 {
2813                         efa->e1->tmp.l++;
2814                         efa->e2->tmp.l++;
2815                         efa->e3->tmp.l++;
2816                         
2817                         if (efa->e4)
2818                         {
2819                                 efa->e4->tmp.l++;
2820                         }
2821                 }
2822         }
2823
2824         /* Add faces angle to the edge weight */
2825         for(efa = em->faces.first; efa; efa = efa->next)
2826         {
2827                 if (efa->h == 0)
2828                 {
2829                         if (efa->v4 == NULL)
2830                         {
2831                                 addTriangle(efa->v1, efa->v2, efa->v3, efa->e1->tmp.l, efa->e2->tmp.l, efa->e3->tmp.l);
2832                         }
2833                         else
2834                         {
2835                                 addTriangle(efa->v1, efa->v2, efa->v3, efa->e1->tmp.l, efa->e2->tmp.l, 2);
2836                                 addTriangle(efa->v3, efa->v4, efa->v1, efa->e3->tmp.l, efa->e4->tmp.l, 2);
2837                         }
2838                 }
2839         }
2840         
2841         nlEnd(NL_MATRIX);
2842
2843         nlEnd(NL_SYSTEM);
2844
2845         success = nlSolveAdvanced(NULL, NL_TRUE);
2846
2847         if (success)
2848         {
2849                 rval = 1;
2850                 for(index = 0, eve = em->verts.first; eve; index++, eve = eve->next)
2851                 {
2852                         weightSetData(eve, nlGetVariable(0, index));
2853                 }
2854         }
2855         else
2856         {
2857                 rval = 0;
2858         }
2859
2860         nlDeleteContext(nlGetCurrent());
2861
2862         return rval;
2863 }
2864
2865
2866 EditEdge * NextEdgeForVert(EdgeIndex *indexed_edges, int index)
2867 {
2868         static int offset = -1;
2869         
2870         /* Reset method, call with NULL mesh pointer */
2871         if (index == -1)
2872         {
2873                 offset = -1;
2874                 return NULL;
2875         }
2876         
2877         /* first pass, start at the head of the list */
2878         if (offset == -1)
2879         {
2880                 offset = indexed_edges->offset[index];
2881         }
2882         /* subsequent passes, start on the next edge */
2883         else
2884         {
2885                 offset++;
2886         }
2887         
2888         return indexed_edges->edges[offset];
2889 }
2890
2891 static void shortestPathsFromVert(EditMesh *em, EditVert *starting_vert, EdgeIndex *indexed_edges)
2892 {
2893         Heap     *edge_heap;
2894         EditVert *current_eve = NULL;
2895         EditEdge *eed = NULL;
2896         EditEdge *select_eed = NULL;
2897         
2898         edge_heap = BLI_heap_new();
2899         
2900         current_eve = starting_vert;
2901         
2902         /* insert guard in heap, when that is returned, no more edges */
2903         BLI_heap_insert(edge_heap, FLT_MAX, NULL);
2904
2905         /* Initialize edge flag */
2906         for(eed= em->edges.first; eed; eed= eed->next)
2907         {
2908                 eed->f1 = 0;
2909         }
2910         
2911         while (BLI_heap_size(edge_heap) > 0)
2912         {
2913                 float current_weight;
2914                 
2915                 current_eve->f1 = 1; /* mark vertex as selected */
2916                 
2917                 /* Add all new edges connected to current_eve to the list */
2918                 NextEdgeForVert(indexed_edges, -1); // Reset next edge
2919                 for(eed = NextEdgeForVert(indexed_edges, indexData(current_eve)); eed; eed = NextEdgeForVert(indexed_edges, indexData(current_eve)))
2920                 { 
2921                         if (eed->f1 == 0)
2922                         {
2923                                 BLI_heap_insert(edge_heap, weightData(current_eve) + eed->tmp.fp, eed);
2924                                 eed->f1 = 1;
2925                         }
2926                 }
2927                 
2928                 /* Find next shortest edge with unselected verts */
2929                 do
2930                 {
2931                         current_weight = BLI_heap_node_value(BLI_heap_top(edge_heap));
2932                         select_eed = BLI_heap_popmin(edge_heap);
2933                 } while (select_eed != NULL && select_eed->v1->f1 != 0 && select_eed->v2->f1);
2934                 
2935                 if (select_eed != NULL)
2936                 {
2937                         select_eed->f1 = 2;
2938                         
2939                         if (select_eed->v1->f1 == 0) /* v1 is the new vertex */
2940                         {
2941                                 current_eve = select_eed->v1;
2942                         }
2943                         else /* otherwise, it's v2 */
2944                         {
2945                                 current_eve = select_eed->v2;
2946                         }
2947                         
2948                         weightSetData(current_eve, current_weight);
2949                 }
2950         }
2951         
2952         BLI_heap_free(edge_heap, NULL);
2953 }
2954
2955 static void freeEdgeIndex(EdgeIndex *indexed_edges)
2956 {
2957         MEM_freeN(indexed_edges->offset);
2958         MEM_freeN(indexed_edges->edges);
2959 }
2960
2961 static void buildIndexedEdges(EditMesh *em, EdgeIndex *indexed_edges)
2962 {
2963         EditVert *eve;
2964         EditEdge *eed;
2965         int totvert = 0;
2966         int tot_indexed = 0;
2967         int offset = 0;
2968
2969         totvert = BLI_countlist(&em->verts);
2970
2971         indexed_edges->offset = MEM_callocN(totvert * sizeof(int), "EdgeIndex offset");
2972
2973         for(eed = em->edges.first; eed; eed = eed->next)
2974         {
2975                 if (eed->v1->h == 0 && eed->v2->h == 0)
2976                 {
2977                         tot_indexed += 2;
2978                         indexed_edges->offset[indexData(eed->v1)]++;
2979                         indexed_edges->offset[indexData(eed->v2)]++;
2980                 }
2981         }
2982         
2983         tot_indexed += totvert;
2984
2985         indexed_edges->edges = MEM_callocN(tot_indexed * sizeof(EditEdge*), "EdgeIndex edges");
2986
2987         /* setting vert offsets */
2988         for(eve = em->verts.first; eve; eve = eve->next)
2989         {
2990                 if (eve->h == 0)
2991                 {
2992                         int d = indexed_edges->offset[indexData(eve)];
2993                         indexed_edges->offset[indexData(eve)] = offset;
2994                         offset += d + 1;
2995                 }
2996         }
2997
2998         /* adding edges in array */
2999         for(eed = em->edges.first; eed; eed= eed->next)
3000         {
3001                 if (eed->v1->h == 0 && eed->v2->h == 0)
3002                 {
3003                         int i;
3004                         for (i = indexed_edges->offset[indexData(eed->v1)]; i < tot_indexed; i++)
3005                         {
3006                                 if (indexed_edges->edges[i] == NULL)
3007                                 {
3008                                         indexed_edges->edges[i] = eed;
3009                                         break;
3010                                 }
3011                         }
3012                         
3013                         for (i = indexed_edges->offset[indexData(eed->v2)]; i < tot_indexed; i++)
3014                         {
3015                                 if (indexed_edges->edges[i] == NULL)
3016                                 {
3017                                         indexed_edges->edges[i] = eed;
3018                                         break;
3019                                 }
3020                         }
3021                 }
3022         }
3023 }
3024
3025 int weightFromDistance(EditMesh *em, EdgeIndex *indexed_edges)
3026 {
3027         EditVert *eve;
3028         int totedge = 0;
3029         int totvert = 0;
3030         int vCount = 0;
3031         
3032         totvert = BLI_countlist(&em->verts);
3033         
3034         if (em == NULL || totvert == 0)
3035         {
3036                 return 0;
3037         }
3038         
3039         totedge = BLI_countlist(&em->edges);
3040         
3041         if (totedge == 0)
3042         {
3043                 return 0;
3044         }
3045         
3046         /* Initialize vertice flag and find at least one selected vertex */
3047         for(eve = em->verts.first; eve; eve = eve->next)
3048         {
3049                 eve->f1 = 0;
3050                 if (eve->f & SELECT)
3051                 {
3052                         vCount = 1;
3053                 }
3054         }
3055         
3056         if (vCount == 0)
3057         {
3058                 return 0; /* no selected vert, failure */
3059         }
3060         else
3061         {
3062                 EditEdge *eed;
3063                 int allDone = 0;
3064
3065                 /* Calculate edge weight */
3066                 for(eed = em->edges.first; eed; eed= eed->next)
3067                 {
3068                         if (eed->v1->h == 0 && eed->v2->h == 0)
3069                         {
3070                                 eed->tmp.fp = len_v3v3(eed->v1->co, eed->v2->co);
3071                         }
3072                 }
3073
3074                 /* Apply dijkstra spf for each selected vert */
3075                 for(eve = em->verts.first; eve; eve = eve->next)
3076                 {
3077                         if (eve->f & SELECT)
3078                         {
3079                                 shortestPathsFromVert(em, eve, indexed_edges);                          
3080                         }
3081                 }
3082                 
3083                 /* connect unselected islands */
3084                 while (allDone == 0)
3085                 {
3086                         EditVert *selected_eve = NULL;
3087                         float selected_weight = 0;
3088                         float min_distance = FLT_MAX;
3089                         
3090                         allDone = 1;
3091                         
3092                         for (eve = em->verts.first; eve; eve = eve->next)
3093                         {
3094                                 /* for every vertex visible that hasn't been processed yet */
3095                                 if (eve->h == 0 && eve->f1 != 1)
3096                                 {
3097                                         EditVert *closest_eve;
3098                                         
3099                                         /* find the closest processed vertex */
3100                                         for (closest_eve = em->verts.first; closest_eve; closest_eve = closest_eve->next)
3101                                         {
3102                                                 /* vertex is already processed and distance is smaller than current minimum */
3103                                                 if (closest_eve->f1 == 1)
3104                                                 {
3105                                                         float distance = len_v3v3(closest_eve->co, eve->co);
3106                                                         if (distance < min_distance)
3107                                                         {
3108                                                                 min_distance = distance;
3109                                                                 selected_eve = eve;
3110                                                                 selected_weight = weightData(closest_eve);
3111                                                         }
3112                                                 }
3113                                         }
3114                                 }
3115                         }
3116                         
3117                         if (selected_eve)
3118                         {
3119                                 allDone = 0;
3120
3121                                 weightSetData(selected_eve, selected_weight + min_distance);
3122                                 shortestPathsFromVert(em, selected_eve, indexed_edges);
3123                         }
3124                 }
3125         }
3126
3127         for(eve = em->verts.first; eve && vCount == 0; eve = eve->next)
3128         {
3129                 if (eve->f1 == 0)
3130                 {
3131                         printf("vertex not reached\n");
3132                         break;
3133                 }
3134         }
3135         
3136         return 1;
3137 }
3138
3139 /****************************************** BUCKET ITERATOR **************************************************/
3140
3141 static void* headNode(void *arg);
3142 static void* tailNode(void *arg);
3143 static void* nextBucket(void *arg);
3144 static void* nextNBucket(void *arg, int n);
3145 static void* peekBucket(void *arg, int n);
3146 static void* previousBucket(void *arg);
3147 static int   iteratorStopped(void *arg);
3148
3149 static void initIteratorFct(ReebArcIterator *iter)
3150 {
3151         iter->head = headNode;
3152         iter->tail = tailNode;
3153         iter->peek = peekBucket;
3154         iter->next = nextBucket;
3155         iter->nextN = nextNBucket;
3156         iter->previous = previousBucket;
3157         iter->stopped = iteratorStopped;        
3158 }
3159
3160 static void setIteratorValues(ReebArcIterator *iter, EmbedBucket *bucket)
3161 {
3162         if (bucket)
3163         {
3164                 iter->p = bucket->p;
3165                 iter->no = bucket->no;
3166         }
3167         else
3168         {
3169                 iter->p = NULL;
3170                 iter->no = NULL;
3171         }
3172         iter->size = 0;
3173 }
3174
3175 void initArcIterator(BArcIterator *arg, ReebArc *arc, ReebNode *head)
3176 {
3177         ReebArcIterator *iter = (ReebArcIterator*)arg;
3178
3179         initIteratorFct(iter);
3180         iter->arc = arc;
3181         
3182         if (head == arc->head)
3183         {
3184                 iter->start = 0;
3185                 iter->end = arc->bcount - 1;
3186                 iter->stride = 1;
3187         }
3188         else
3189         {
3190                 iter->start = arc->bcount - 1;
3191                 iter->end = 0;
3192                 iter->stride = -1;
3193         }
3194         
3195         iter->length = arc->bcount;
3196         
3197         iter->index = -1;
3198 }
3199
3200 void initArcIteratorStart(BArcIterator *arg, struct ReebArc *arc, struct ReebNode *head, int start)
3201 {
3202         ReebArcIterator *iter = (ReebArcIterator*)arg;
3203
3204         initIteratorFct(iter);
3205         iter->arc = arc;
3206         
3207         if (head == arc->head)
3208         {
3209                 iter->start = start;
3210                 iter->end = arc->bcount - 1;
3211                 iter->stride = 1;
3212         }
3213         else
3214         {
3215                 iter->start = arc->bcount - 1 - start;
3216                 iter->end = 0;
3217                 iter->stride = -1;
3218         }
3219         
3220         iter->index = -1;
3221