source/blender/editors/armature/reeb.c

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