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