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