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