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