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