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