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