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