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