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