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