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