merge with trunk (15330 -> 15566)
authorMartin Poirier <theeth@yahoo.com>
Mon, 14 Jul 2008 14:09:36 +0000 (14:09 +0000)
committerMartin Poirier <theeth@yahoo.com>
Mon, 14 Jul 2008 14:09:36 +0000 (14:09 +0000)
18 files changed:
release/scripts/reeb.py [new file with mode: 0644]
source/blender/blenlib/BLI_ghash.h
source/blender/blenlib/BLI_graph.h [new file with mode: 0644]
source/blender/blenlib/intern/BLI_ghash.c
source/blender/blenlib/intern/graph.c [new file with mode: 0644]
source/blender/include/BIF_editarmature.h
source/blender/include/butspace.h
source/blender/include/reeb.h
source/blender/makesdna/DNA_scene_types.h
source/blender/src/autoarmature.c [new file with mode: 0644]
source/blender/src/buttons_editing.c
source/blender/src/drawview.c
source/blender/src/editarmature.c
source/blender/src/reeb.c
source/blender/src/usiblender.c
source/gameengine/Converter/BL_ShapeActionActuator.cpp
source/gameengine/Converter/BL_ShapeDeformer.cpp
source/gameengine/Converter/BL_ShapeDeformer.h

diff --git a/release/scripts/reeb.py b/release/scripts/reeb.py
new file mode 100644 (file)
index 0000000..63ab120
--- /dev/null
@@ -0,0 +1,110 @@
+#!BPY
+"""
+Name: 'Reeb graph import'
+Blender: 245
+Group: 'Import'
+Tooltip: 'Imports a reeb graph saved after skeleton generation'
+"""
+import Blender
+
+def name(count):
+       if count == -1:
+               return ""
+       else:
+               return "%05" % count
+
+def importGraph(count):
+       bNode = Blender.Draw.Create(1)
+       bSize = Blender.Draw.Create(0.01)
+
+       Block = []
+       
+       Block.append(("Size: ", bSize, 0.01, 10.0, "Size of the nodes"))
+       Block.append(("Nodes", bNode, "Import nodes as tetras"))
+       
+       retval = Blender.Draw.PupBlock("Reeb Graph Import", Block)
+       
+       if not retval:
+               return
+
+
+       me = Blender.Mesh.New("graph%s" % name(count))
+       scn = Blender.Scene.GetCurrent()
+       
+       f = open("test%s.txt" % name(count), "r")
+       
+       verts = []
+       edges = []
+       faces = []
+       
+       i = 0
+       first = False
+       
+       SIZE = float(bSize.val)
+       WITH_NODE = bool(bNode.val)
+       
+       def addNode(v, s, verts, faces):
+               if WITH_NODE:
+                       v1 = [v[0], v[1], v[2] + s]
+                       i1 = len(verts)
+                       verts.append(v1)
+                       v2 = [v[0], v[1] + 0.959 * s, v[2] - 0.283 * s]
+                       i2 = len(verts)
+                       verts.append(v2)
+                       v3 = [v[0] - 0.830 * s, v[1] - 0.479 * s, v[2] - 0.283 * s]
+                       i3 = len(verts)
+                       verts.append(v3)
+                       v4 = [v[0] + 0.830 * s, v[1] - 0.479 * s, v[2] - 0.283 * s]
+                       i4 = len(verts)
+                       verts.append(v4)
+                       
+                       faces.append([i1,i2,i3])
+                       faces.append([i1,i3,i4])
+                       faces.append([i2,i3,i4])
+                       faces.append([i1,i2,i4])
+                       
+                       return 4
+               else:
+                       return 0
+               
+       for line in f:
+               data = line.strip().split(" ")
+               if data[0] == "v1":
+                       v = [float(x) for x in data[-3:]]
+                       i += addNode(v, SIZE, verts, faces)
+                       verts.append(v)
+                       i += 1
+               elif data[0] == "v2":
+                       pass
+                       v = [float(x) for x in data[-3:]]
+                       verts.append(v)
+                       edges.append((i-1, i))
+                       i += 1
+                       i += addNode(v, SIZE, verts, faces)
+               elif data[0] == "b":
+                       verts.append([float(x) for x in data[-3:]])
+                       edges.append((i-1, i))
+                       i += 1
+#              elif data[0] == "angle":
+#                      obj = scn.objects.new('Empty')
+#                      obj.loc = (float(data[1]), float(data[2]), float(data[3]))
+#                      obj.properties["angle"] = data[4]
+#                      del obj
+                        
+                        
+       me.verts.extend(verts)
+       me.edges.extend(edges)
+       me.faces.extend(faces)
+       
+       
+       ob = scn.objects.new(me, "graph%s" % name(count))
+       del ob
+       del scn
+       
+
+#for i in range(16):
+#      importGraph(i)
+
+if __name__=='__main__':
+       importGraph(-1)
index aec77f5f3859ddcd684cc31c3af321cc18b5e25a..c77e82f0a2bbcb3ea6fc8830bcadf660cbd3dde1 100644 (file)
 
 struct GHash;
 typedef struct GHash GHash;
-typedef struct GHashIterator GHashIterator;
+
+typedef struct GHashIterator {
+       GHash *gh;
+       int curBucket;
+       struct Entry *curEntry;
+} GHashIterator;
 
 typedef unsigned int   (*GHashHashFP)          (void *key);
 typedef int                            (*GHashCmpFP)           (void *a, void *b);
@@ -62,6 +67,15 @@ int          BLI_ghash_size          (GHash *gh);
         * @return Pointer to a new DynStr.
         */
 GHashIterator* BLI_ghashIterator_new           (GHash *gh);
+       /**
+        * Init an already allocated GHashIterator. The hash table must not
+        * be mutated while the iterator is in use, and the iterator will
+        * step exactly BLI_ghash_size(gh) times before becoming done.
+        * 
+        * @param ghi The GHashIterator to initialize.
+        * @param gh The GHash to iterate over.
+        */
+void BLI_ghashIterator_init(GHashIterator *ghi, GHash *gh);
        /**
         * Free a GHashIterator.
         *
diff --git a/source/blender/blenlib/BLI_graph.h b/source/blender/blenlib/BLI_graph.h
new file mode 100644 (file)
index 0000000..d309d73
--- /dev/null
@@ -0,0 +1,114 @@
+#ifndef BLI_GRAPH_H_
+#define BLI_GRAPH_H_
+
+#include "DNA_listBase.h"
+
+struct BGraph;
+struct BNode;
+struct BArc;
+
+struct RadialArc;
+
+typedef void (*FreeArc)(struct BArc*);
+typedef void (*FreeNode)(struct BNode*);
+typedef void (*RadialSymmetry)(struct BNode* root_node, struct RadialArc* ring, int total);
+typedef void (*AxialSymmetry)(struct BNode* root_node, struct BNode* node1, struct BNode* node2, struct BArc* arc1, struct BArc* arc2);
+
+/* IF YOU MODIFY THOSE TYPES, YOU NEED TO UPDATE ALL THOSE THAT "INHERIT" FROM THEM
+ * 
+ * RigGraph, ReebGraph
+ * 
+ * */
+
+typedef struct BGraph {
+       ListBase        arcs;
+       ListBase        nodes;
+       
+       float length;
+       
+       /* function pointer to deal with custom fonctionnality */
+       FreeArc                 free_arc;
+       FreeNode                free_node;
+       RadialSymmetry  radial_symmetry;
+       AxialSymmetry   axial_symmetry;
+} BGraph;
+
+typedef struct BNode {
+       void *next, *prev;
+       float p[3];
+       int flag;
+
+       int degree;
+       struct BArc **arcs;
+
+       int symmetry_level;
+       int symmetry_flag;
+       float symmetry_axis[3];
+} BNode;
+
+typedef struct BArc {
+       void *next, *prev;
+       struct BNode *head, *tail;
+       int flag;
+
+       float length;
+
+       int symmetry_level;
+       int symmetry_group;
+       int symmetry_flag;
+} BArc;
+
+/* Helper structure for radial symmetry */
+typedef struct RadialArc
+{
+       struct BArc *arc; 
+       float n[3]; /* normalized vector joining the nodes of the arc */
+} RadialArc;
+
+BNode *BLI_otherNode(BArc *arc, BNode *node);
+
+void BLI_freeNode(BGraph *graph, BNode *node);
+void BLI_removeNode(BGraph *graph, BNode *node);
+
+void BLI_flagNodes(BGraph *graph, int flag);
+void BLI_flagArcs(BGraph *graph, int flag);
+
+int BLI_hasAdjacencyList(BGraph *rg);
+void BLI_buildAdjacencyList(BGraph *rg);
+void BLI_rebuildAdjacencyList(BGraph* rg);
+void BLI_freeAdjacencyList(BGraph *rg);
+
+int BLI_FlagSubgraphs(BGraph *graph);
+
+int BLI_subtreeShape(BNode *node, BArc *rootArc, int include_root);
+float BLI_subtreeLength(BNode *node, BArc *rootArc);
+void BLI_calcGraphLength(BGraph *graph);
+
+void BLI_replaceNode(BGraph *graph, BNode *node_src, BNode *node_replaced);
+void BLI_removeDoubleNodes(BGraph *graph, float limit);
+
+BArc * BLI_findConnectedArc(BGraph *graph, BArc *arc, BNode *v);
+
+int    BLI_isGraphCyclic(BGraph *graph);
+
+/*------------ Symmetry handling ------------*/
+void BLI_markdownSymmetry(BGraph *graph, BNode *root_node, float limit);
+
+void BLI_mirrorAlongAxis(float v[3], float center[3], float axis[3]);
+
+/* BNode symmetry flags */
+#define SYM_TOPOLOGICAL        1
+#define SYM_PHYSICAL   2
+
+/* the following two are exclusive */
+#define SYM_AXIAL              4
+#define SYM_RADIAL             8
+
+/* BArc symmetry flags
+ * 
+ * axial symetry sides */
+#define SYM_SIDE_POSITIVE              1
+#define SYM_SIDE_NEGATIVE              2
+/* Anything higher is the order in radial symmetry */
+
+#endif /*BLI_GRAPH_H_*/
index 227cb8f5e9a5e23c95f2c446f01bbe8020691d12..bae6ad428ea39c0f6cde431c2e161d7b309a4832 100644 (file)
@@ -198,12 +198,6 @@ void BLI_ghash_free(GHash *gh, GHashKeyFreeFP keyfreefp, GHashValFreeFP valfreef
 
 /***/
 
-struct GHashIterator {
-       GHash *gh;
-       int curBucket;
-       Entry *curEntry;
-};
-
 GHashIterator *BLI_ghashIterator_new(GHash *gh) {
        GHashIterator *ghi= malloc(sizeof(*ghi));
        ghi->gh= gh;
@@ -217,6 +211,17 @@ GHashIterator *BLI_ghashIterator_new(GHash *gh) {
        }
        return ghi;
 }
+void BLI_ghashIterator_init(GHashIterator *ghi, GHash *gh) {
+       ghi->gh= gh;
+       ghi->curEntry= NULL;
+       ghi->curBucket= -1;
+       while (!ghi->curEntry) {
+               ghi->curBucket++;
+               if (ghi->curBucket==ghi->gh->nbuckets)
+                       break;
+               ghi->curEntry= ghi->gh->buckets[ghi->curBucket];
+       }
+}
 void BLI_ghashIterator_free(GHashIterator *ghi) {
        free(ghi);
 }
diff --git a/source/blender/blenlib/intern/graph.c b/source/blender/blenlib/intern/graph.c
new file mode 100644 (file)
index 0000000..a308936
--- /dev/null
@@ -0,0 +1,964 @@
+/**
+ * $Id:
+ *
+ * ***** BEGIN GPL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+ *
+ * Contributor(s): Martin Poirier
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ * graph.c: Common graph interface and methods
+ */
+
+#include <float.h> 
+#include <math.h>
+
+#include "MEM_guardedalloc.h"
+
+#include "BLI_graph.h"
+#include "BLI_blenlib.h"
+#include "BLI_arithb.h"
+
+#include "BKE_utildefines.h"
+
+static void testRadialSymmetry(BGraph *graph, BNode* root_node, RadialArc* ring, int total, float axis[3], float limit, int group);
+
+static void handleAxialSymmetry(BGraph *graph, BNode *root_node, int depth, float axis[3], float limit);
+static void testAxialSymmetry(BGraph *graph, BNode* root_node, BNode* node1, BNode* node2, BArc* arc1, BArc* arc2, float axis[3], float limit, int group);
+static void flagAxialSymmetry(BNode *root_node, BNode *end_node, BArc *arc, int group);
+
+void BLI_freeNode(BGraph *graph, BNode *node)
+{
+       if (node->arcs)
+       {
+               MEM_freeN(node->arcs);
+       }
+       
+       if (graph->free_node)
+       {
+               graph->free_node(node);
+       }
+}
+
+void BLI_removeNode(BGraph *graph, BNode *node)
+{
+       BLI_freeNode(graph, node);
+       BLI_freelinkN(&graph->nodes, node);
+}
+
+BNode *BLI_otherNode(BArc *arc, BNode *node)
+{
+       return (arc->head == node) ? arc->tail : arc->head;
+}
+
+void BLI_flagNodes(BGraph *graph, int flag)
+{
+       BNode *node;
+       
+       for(node = graph->nodes.first; node; node = node->next)
+       {
+               node->flag = flag;
+       }
+}
+
+void BLI_flagArcs(BGraph *graph, int flag)
+{
+       BArc *arc;
+       
+       for(arc = graph->arcs.first; arc; arc = arc->next)
+       {
+               arc->flag = flag;
+       }
+}
+
+static void addArcToNodeAdjacencyList(BNode *node, BArc *arc)
+{
+       node->arcs[node->flag] = arc;
+       node->flag++;
+}
+
+void BLI_rebuildAdjacencyList(BGraph *rg)
+{
+       BLI_freeAdjacencyList(rg);
+       BLI_buildAdjacencyList(rg);
+}
+
+void BLI_buildAdjacencyList(BGraph *rg)
+{
+       BNode *node;
+       BArc *arc;
+
+       for(node = rg->nodes.first; node; node = node->next)
+       {
+               if (node->arcs != NULL)
+               {
+                       MEM_freeN(node->arcs);
+               }
+               
+               node->arcs = MEM_callocN((node->degree) * sizeof(BArc*), "adjacency list");
+               
+               /* temporary use to indicate the first index available in the lists */
+               node->flag = 0;
+       }
+
+       for(arc = rg->arcs.first; arc; arc= arc->next)
+       {
+               addArcToNodeAdjacencyList(arc->head, arc);
+               addArcToNodeAdjacencyList(arc->tail, arc);
+       }
+
+       for(node = rg->nodes.first; node; node = node->next)
+       {
+               if (node->degree != node->flag)
+               {
+                       printf("error in node [%p]. Added only %i arcs out of %i\n", node, node->flag, node->degree);
+               }
+       }
+}
+
+void BLI_freeAdjacencyList(BGraph *rg)
+{
+       BNode *node;
+
+       for(node = rg->nodes.first; node; node = node->next)
+       {
+               if (node->arcs != NULL)
+               {
+                       MEM_freeN(node->arcs);
+                       node->arcs = NULL;
+               }
+       }
+}
+
+int BLI_hasAdjacencyList(BGraph *rg)
+{
+       BNode *node;
+       
+       for(node = rg->nodes.first; node; node = node->next)
+       {
+               if (node->arcs == NULL)
+               {
+                       return 0;
+               }
+       }
+       
+       return 1;
+}
+
+void BLI_replaceNode(BGraph *graph, BNode *node_src, BNode *node_replaced)
+{
+       BArc *arc, *next_arc;
+       
+       for (arc = graph->arcs.first; arc; arc = next_arc)
+       {
+               next_arc = arc->next;
+               
+               if (arc->head == node_replaced)
+               {
+                       arc->head = node_src;
+                       node_src->degree++;
+               }
+
+               if (arc->tail == node_replaced)
+               {
+                       arc->tail = node_src;
+                       node_src->degree++;
+               }
+               
+               if (arc->head == arc->tail)
+               {
+                       node_src->degree -= 2;
+                       
+                       graph->free_arc(arc);
+                       BLI_freelinkN(&graph->arcs, arc);
+               }
+       }
+}
+
+void BLI_removeDoubleNodes(BGraph *graph, float limit)
+{
+       BNode *node_src, *node_replaced;
+       
+       for(node_src = graph->nodes.first; node_src; node_src = node_src->next)
+       {
+               for(node_replaced = graph->nodes.first; node_replaced; node_replaced = node_replaced->next)
+               {
+                       if (node_replaced != node_src && VecLenf(node_replaced->p, node_src->p) <= limit)
+                       {
+                               BLI_replaceNode(graph, node_src, node_replaced);
+                               
+                               BLI_removeNode(graph, node_replaced);
+                       }
+               }
+       }
+       
+}
+/************************************* SUBGRAPH DETECTION **********************************************/
+
+void flagSubgraph(BNode *node, int subgraph)
+{
+       if (node->flag == 0)
+       {
+               BArc *arc;
+               int i;
+               
+               node->flag = subgraph;
+               
+               for(i = 0; i < node->degree; i++)
+               {
+                       arc = node->arcs[i];
+                       flagSubgraph(BLI_otherNode(arc, node), subgraph);
+               }
+       }
+} 
+
+int BLI_FlagSubgraphs(BGraph *graph)
+{
+       BNode *node;
+       int subgraph = 0;
+
+       if (BLI_hasAdjacencyList(graph) == 0)
+       {
+               BLI_buildAdjacencyList(graph);
+       }
+       
+       BLI_flagNodes(graph, 0);
+       
+       for (node = graph->nodes.first; node; node = node->next)
+       {
+               if (node->flag == 0)
+               {
+                       subgraph++;
+                       flagSubgraph(node, subgraph);
+               }
+       }
+       
+       return subgraph;
+}
+
+/*************************************** CYCLE DETECTION ***********************************************/
+
+int detectCycle(BNode *node, BArc *src_arc)
+{
+       int value = 0;
+       
+       if (node->flag == 0)
+       {
+               int i;
+
+               /* mark node as visited */
+               node->flag = 1;
+
+               for(i = 0; i < node->degree && value == 0; i++)
+               {
+                       BArc *arc = node->arcs[i];
+                       
+                       /* don't go back on the source arc */
+                       if (arc != src_arc)
+                       {
+                               value = detectCycle(BLI_otherNode(arc, node), arc);
+                       }
+               }
+       }
+       else
+       {
+               value = 1;
+       }
+       
+       return value;
+}
+
+int    BLI_isGraphCyclic(BGraph *graph)
+{
+       BNode *node;
+       int value = 0;
+       
+       /* NEED TO CHECK IF ADJACENCY LIST EXIST */
+       
+       /* Mark all nodes as not visited */
+       BLI_flagNodes(graph, 0);
+
+       /* detectCycles in subgraphs */ 
+       for(node = graph->nodes.first; node && value == 0; node = node->next)
+       {
+               /* only for nodes in subgraphs that haven't been visited yet */
+               if (node->flag == 0)
+               {
+                       value = value || detectCycle(node, NULL);
+               }               
+       }
+       
+       return value;
+}
+
+BArc * BLI_findConnectedArc(BGraph *graph, BArc *arc, BNode *v)
+{
+       BArc *nextArc = arc->next;
+       
+       for(nextArc = graph->arcs.first; nextArc; nextArc = nextArc->next)
+       {
+               if (arc != nextArc && (nextArc->head == v || nextArc->tail == v))
+               {
+                       break;
+               }
+       }
+       
+       return nextArc;
+}
+
+/*********************************** GRAPH AS TREE FUNCTIONS *******************************************/
+
+int BLI_subtreeShape(BNode *node, BArc *rootArc, int include_root)
+{
+       int depth = 0;
+       
+       if (include_root)
+       {
+               BNode *newNode = BLI_otherNode(rootArc, node);
+               return BLI_subtreeShape(newNode, rootArc, 0);
+       }
+       else
+       {
+               /* Base case, no arcs leading away */
+               if (node->arcs == NULL || *(node->arcs) == NULL)
+               {
+                       return 0;
+               }
+               else
+               {
+                       int i;
+       
+                       for(i = 0; i < node->degree; i++)
+                       {
+                               BArc *arc = node->arcs[i];
+                               
+                               /* only arcs that go down the tree */
+                               if (arc != rootArc)
+                               {
+                                       BNode *newNode = BLI_otherNode(arc, node);
+                                       depth += BLI_subtreeShape(newNode, arc, 0);
+                               }
+                       }
+               }
+               
+               return 10 * depth + 1;
+       }
+}
+
+float BLI_subtreeLength(BNode *node, BArc *rootArc)
+{
+       float length = 0;
+       int i;
+
+       for(i = 0; i < node->degree; i++)
+       {
+               BArc *arc = node->arcs[i];
+               
+               /* don't go back on the root arc */
+               if (arc != rootArc)
+               {
+                       length = MAX2(length, BLI_subtreeLength(BLI_otherNode(arc, node), arc));
+               }
+       }
+       
+       if (rootArc)
+       {
+               length += rootArc->length;
+       }
+       
+       return length;
+}
+
+void BLI_calcGraphLength(BGraph *graph)
+{
+       if (BLI_isGraphCyclic(graph) == 0)
+       {
+               float length = 0;
+               int nb_subgraphs;
+               int i;
+               
+               nb_subgraphs = BLI_FlagSubgraphs(graph);
+               
+               for (i = 1; i <= nb_subgraphs; i++)
+               {
+                       BNode *node;
+                       
+                       for (node = graph->nodes.first; node; node = node->next)
+                       {
+                               /* start on an external node  of the subgraph */
+                               if (node->flag == i && node->degree == 1)
+                               {
+                                       length = MAX2(length, BLI_subtreeLength(node, NULL));
+                                       break;
+                               }
+                       }
+               }
+               
+               graph->length = length;
+       }
+}
+
+/********************************* SYMMETRY DETECTION **************************************************/
+
+void markdownSymmetryArc(BGraph *graph, BArc *arc, BNode *node, int level, float limit);
+
+void BLI_mirrorAlongAxis(float v[3], float center[3], float axis[3])
+{
+       float dv[3], pv[3];
+       
+       VecSubf(dv, v, center);
+       Projf(pv, dv, axis);
+       VecMulf(pv, -2);
+       VecAddf(v, v, pv);
+}
+
+static void testRadialSymmetry(BGraph *graph, BNode* root_node, RadialArc* ring, int total, float axis[3], float limit, int group)
+{
+       int symmetric = 1;
+       int i;
+       
+       /* sort ring by angle */
+       for (i = 0; i < total - 1; i++)
+       {
+               float minAngle = FLT_MAX;
+               int minIndex = -1;
+               int j;
+
+               for (j = i + 1; j < total; j++)
+               {
+                       float angle = Inpf(ring[i].n, ring[j].n);
+
+                       /* map negative values to 1..2 */
+                       if (angle < 0)
+                       {
+                               angle = 1 - angle;
+                       }
+
+                       if (angle < minAngle)
+                       {
+                               minIndex = j;
+                               minAngle = angle;
+                       }
+               }
+
+               /* swap if needed */
+               if (minIndex != i + 1)
+               {
+                       RadialArc tmp;
+                       tmp = ring[i + 1];
+                       ring[i + 1] = ring[minIndex];
+                       ring[minIndex] = tmp;
+               }
+       }
+
+       for (i = 0; i < total && symmetric; i++)
+       {
+               BNode *node1, *node2;
+               float tangent[3];
+               float normal[3];
+               float p[3];
+               int j = (i + 1) % total; /* next arc in the circular list */
+
+               VecAddf(tangent, ring[i].n, ring[j].n);
+               Crossf(normal, tangent, axis);
+               
+               node1 = BLI_otherNode(ring[i].arc, root_node);
+               node2 = BLI_otherNode(ring[j].arc, root_node);
+
+               VECCOPY(p, node2->p);
+               BLI_mirrorAlongAxis(p, root_node->p, normal);
+               
+               /* check if it's within limit before continuing */
+               if (VecLenf(node1->p, p) > limit)
+               {
+                       symmetric = 0;
+               }
+
+       }
+
+       if (symmetric)
+       {
+               /* mark node as symmetric physically */
+               VECCOPY(root_node->symmetry_axis, axis);
+               root_node->symmetry_flag |= SYM_PHYSICAL;
+               root_node->symmetry_flag |= SYM_RADIAL;
+               
+               /* FLAG SYMMETRY GROUP */
+               for (i = 0; i < total; i++)
+               {
+                       ring[i].arc->symmetry_group = group;
+                       ring[i].arc->symmetry_flag = i;
+               }
+
+               if (graph->radial_symmetry)
+               {
+                       graph->radial_symmetry(root_node, ring, total);
+               }
+       }
+}
+
+static void handleRadialSymmetry(BGraph *graph, BNode *root_node, int depth, float axis[3], float limit)
+{
+       RadialArc *ring = NULL;
+       RadialArc *unit;
+       int total = 0;
+       int group;
+       int first;
+       int i;
+
+       /* mark topological symmetry */
+       root_node->symmetry_flag |= SYM_TOPOLOGICAL;
+
+       /* total the number of arcs in the symmetry ring */
+       for (i = 0; i < root_node->degree; i++)
+       {
+               BArc *connectedArc = root_node->arcs[i];
+               
+               /* depth is store as a negative in flag. symmetry level is positive */
+               if (connectedArc->symmetry_level == -depth)
+               {
+                       total++;
+               }
+       }
+
+       ring = MEM_callocN(sizeof(RadialArc) * total, "radial symmetry ring");
+       unit = ring;
+
+       /* fill in the ring */
+       for (unit = ring, i = 0; i < root_node->degree; i++)
+       {
+               BArc *connectedArc = root_node->arcs[i];
+               
+               /* depth is store as a negative in flag. symmetry level is positive */
+               if (connectedArc->symmetry_level == -depth)
+               {
+                       BNode *otherNode = BLI_otherNode(connectedArc, root_node);
+                       float vec[3];
+
+                       unit->arc = connectedArc;
+
+                       /* project the node to node vector on the symmetry plane */
+                       VecSubf(unit->n, otherNode->p, root_node->p);
+                       Projf(vec, unit->n, axis);
+                       VecSubf(unit->n, unit->n, vec);
+
+                       Normalize(unit->n);
+
+                       unit++;
+               }
+       }
+
+       /* sort ring by arc length
+        * using a rather bogus insertion sort
+        * butrings will never get too big to matter
+        * */
+       for (i = 0; i < total; i++)
+       {
+               int j;
+
+               for (j = i - 1; j >= 0; j--)
+               {
+                       BArc *arc1, *arc2;
+                       
+                       arc1 = ring[j].arc;
+                       arc2 = ring[j + 1].arc;
+                       
+                       if (arc1->length > arc2->length)
+                       {
+                               /* swap with smaller */
+                               RadialArc tmp;
+                               
+                               tmp = ring[j + 1];
+                               ring[j + 1] = ring[j];
+                               ring[j] = tmp;
+                       }
+                       else
+                       {
+                               break;
+                       }
+               }
+       }
+
+       /* Dispatch to specific symmetry tests */
+       first = 0;
+       group = 0;
+       
+       for (i = 1; i < total; i++)
+       {
+               int dispatch = 0;
+               int last = i - 1;
+               
+               if (fabs(ring[first].arc->length - ring[i].arc->length) > limit)
+               {
+                       dispatch = 1;
+               }
+
+               /* if not dispatching already and on last arc
+                * Dispatch using current arc as last
+                * */           
+               if (dispatch == 0 && i == total - 1)
+               {
+                       last = i;
+                       dispatch = 1;
+               } 
+               
+               if (dispatch)
+               {
+                       int sub_total = last - first + 1; 
+
+                       group += 1;
+
+                       if (sub_total == 1)
+                       {
+                               /* NOTHING TO DO */
+                       }
+                       else if (sub_total == 2)
+                       {
+                               BArc *arc1, *arc2;
+                               BNode *node1, *node2;
+                               
+                               arc1 = ring[first].arc;
+                               arc2 = ring[last].arc;
+                               
+                               node1 = BLI_otherNode(arc1, root_node);
+                               node2 = BLI_otherNode(arc2, root_node);
+                               
+                               testAxialSymmetry(graph, root_node, node1, node2, arc1, arc2, axis, limit, group);
+                       }
+                       else if (sub_total != total) /* allocate a new sub ring if needed */
+                       {
+                               RadialArc *sub_ring = MEM_callocN(sizeof(RadialArc) * sub_total, "radial symmetry ring");
+                               int sub_i;
+                               
+                               /* fill in the sub ring */
+                               for (sub_i = 0; sub_i < sub_total; sub_i++)
+                               {
+                                       sub_ring[sub_i] = ring[first + sub_i];
+                               }
+                               
+                               testRadialSymmetry(graph, root_node, sub_ring, sub_total, axis, limit, group);
+                       
+                               MEM_freeN(sub_ring);
+                       }
+                       else if (sub_total == total)
+                       {
+                               testRadialSymmetry(graph, root_node, ring, total, axis, limit, group);
+                       }
+                       
+                       first = i;
+               }
+       }
+
+
+       MEM_freeN(ring);
+}
+
+static void flagAxialSymmetry(BNode *root_node, BNode *end_node, BArc *arc, int group)
+{
+       float vec[3];
+       
+       arc->symmetry_group = group;
+       
+       VecSubf(vec, end_node->p, root_node->p);
+       
+       if (Inpf(vec, root_node->symmetry_axis) < 0)
+       {
+               arc->symmetry_flag |= SYM_SIDE_NEGATIVE;
+       }
+       else
+       {
+               arc->symmetry_flag |= SYM_SIDE_POSITIVE;
+       }
+}
+
+static void testAxialSymmetry(BGraph *graph, BNode* root_node, BNode* node1, BNode* node2, BArc* arc1, BArc* arc2, float axis[3], float limit, int group)
+{
+       float nor[3], vec[3], p[3];
+
+       VecSubf(vec, node1->p, root_node->p);
+       Normalize(vec);
+       VecSubf(p, root_node->p, node2->p);
+       Normalize(p);
+       VecAddf(p, p, vec);
+
+       Crossf(vec, p, axis);
+       Crossf(nor, vec, axis);
+       
+       /* mirror node2 along axis */
+       VECCOPY(p, node2->p);
+       BLI_mirrorAlongAxis(p, root_node->p, nor);
+       
+       /* check if it's within limit before continuing */
+       if (VecLenf(node1->p, p) <= limit)
+       {
+               /* mark node as symmetric physically */
+               VECCOPY(root_node->symmetry_axis, nor);
+               root_node->symmetry_flag |= SYM_PHYSICAL;
+               root_node->symmetry_flag |= SYM_AXIAL;
+
+               /* flag side on arcs */
+               flagAxialSymmetry(root_node, node1, arc1, group);
+               flagAxialSymmetry(root_node, node2, arc2, group);
+               
+               if (graph->axial_symmetry)
+               {
+                       graph->axial_symmetry(root_node, node1, node2, arc1, arc2);
+               }
+       }
+       else
+       {
+               /* NOT SYMMETRIC */
+       }
+}
+
+static void handleAxialSymmetry(BGraph *graph, BNode *root_node, int depth, float axis[3], float limit)
+{
+       BArc *arc1 = NULL, *arc2 = NULL;
+       BNode *node1 = NULL, *node2 = NULL;
+       int i;
+       
+       /* mark topological symmetry */
+       root_node->symmetry_flag |= SYM_TOPOLOGICAL;
+
+       for (i = 0; i < root_node->degree; i++)
+       {
+               BArc *connectedArc = root_node->arcs[i];
+               
+               /* depth is store as a negative in flag. symmetry level is positive */
+               if (connectedArc->symmetry_level == -depth)
+               {
+                       if (arc1 == NULL)
+                       {
+                               arc1 = connectedArc;
+                               node1 = BLI_otherNode(arc1, root_node);
+                       }
+                       else
+                       {
+                               arc2 = connectedArc;
+                               node2 = BLI_otherNode(arc2, root_node);
+                               break; /* Can stop now, the two arcs have been found */
+                       }
+               }
+       }
+       
+       /* shouldn't happen, but just to be sure */
+       if (node1 == NULL || node2 == NULL)
+       {
+               return;
+       }
+       
+       testAxialSymmetry(graph, root_node, node1, node2, arc1, arc2, axis, limit, 1);
+}
+
+static void markdownSecondarySymmetry(BGraph *graph, BNode *node, int depth, int level, float limit)
+{
+       float axis[3] = {0, 0, 0};
+       int count = 0;
+       int i;
+
+       /* count the number of branches in this symmetry group
+        * and determinte the axis of symmetry
+        *  */  
+       for (i = 0; i < node->degree; i++)
+       {
+               BArc *connectedArc = node->arcs[i];
+               
+               /* depth is store as a negative in flag. symmetry level is positive */
+               if (connectedArc->symmetry_level == -depth)
+               {
+                       count++;
+               }
+               /* If arc is on the axis */
+               else if (connectedArc->symmetry_level == level)
+               {
+                       VecAddf(axis, axis, connectedArc->head->p);
+                       VecSubf(axis, axis, connectedArc->tail->p);
+               }
+       }
+
+       Normalize(axis);
+
+       /* Split between axial and radial symmetry */
+       if (count == 2)
+       {
+               handleAxialSymmetry(graph, node, depth, axis, limit);
+       }
+       else
+       {
+               handleRadialSymmetry(graph, node, depth, axis, limit);
+       }
+       
+       /* markdown secondary symetries */      
+       for (i = 0; i < node->degree; i++)
+       {
+               BArc *connectedArc = node->arcs[i];
+               
+               if (connectedArc->symmetry_level == -depth)
+               {
+                       /* markdown symmetry for branches corresponding to the depth */
+                       markdownSymmetryArc(graph, connectedArc, node, level + 1, limit);
+               }
+       }
+}
+
+void markdownSymmetryArc(BGraph *graph, BArc *arc, BNode *node, int level, float limit)
+{
+       int i;
+
+       /* if arc is null, we start straight from a node */     
+       if (arc)
+       {
+               arc->symmetry_level = level;
+               
+               node = BLI_otherNode(arc, node);
+       }
+       
+       for (i = 0; i < node->degree; i++)
+       {
+               BArc *connectedArc = node->arcs[i];
+               
+               if (connectedArc != arc)
+               {
+                       BNode *connectedNode = BLI_otherNode(connectedArc, node);
+                       
+                       /* symmetry level is positive value, negative values is subtree depth */
+                       connectedArc->symmetry_level = -BLI_subtreeShape(connectedNode, connectedArc, 0);
+               }
+       }
+
+       arc = NULL;
+
+       for (i = 0; i < node->degree; i++)
+       {
+               int issymmetryAxis = 0;
+               BArc *connectedArc = node->arcs[i];
+               
+               /* only arcs not already marked as symetric */
+               if (connectedArc->symmetry_level < 0)
+               {
+                       int j;
+                       
+                       /* true by default */
+                       issymmetryAxis = 1;
+                       
+                       for (j = 0; j < node->degree; j++)
+                       {
+                               BArc *otherArc = node->arcs[j];
+                               
+                               /* different arc, same depth */
+                               if (otherArc != connectedArc && otherArc->symmetry_level == connectedArc->symmetry_level)
+                               {
+                                       /* not on the symmetry axis */
+                                       issymmetryAxis = 0;
+                                       break;
+                               } 
+                       }
+               }
+               
+               /* arc could be on the symmetry axis */
+               if (issymmetryAxis == 1)
+               {
+                       /* no arc as been marked previously, keep this one */
+                       if (arc == NULL)
+                       {
+                               arc = connectedArc;
+                       }
+                       else if (connectedArc->symmetry_level < arc->symmetry_level)
+                       {
+                               /* go with more complex subtree as symmetry arc */
+                               arc = connectedArc;
+                       }
+               }
+       }
+       
+       /* go down the arc continuing the symmetry axis */
+       if (arc)
+       {
+               markdownSymmetryArc(graph, arc, node, level, limit);
+       }
+
+       
+       /* secondary symmetry */
+       for (i = 0; i < node->degree; i++)
+       {
+               BArc *connectedArc = node->arcs[i];
+               
+               /* only arcs not already marked as symetric and is not the next arc on the symmetry axis */
+               if (connectedArc->symmetry_level < 0)
+               {
+                       /* subtree depth is store as a negative value in the symmetry */
+                       markdownSecondarySymmetry(graph, node, -connectedArc->symmetry_level, level, limit);
+               }
+       }
+}
+
+void BLI_markdownSymmetry(BGraph *graph, BNode *root_node, float limit)
+{
+       BNode *node;
+       BArc *arc;
+       
+       if (BLI_isGraphCyclic(graph))
+       {
+               return;
+       }
+       
+       /* mark down all arcs as non-symetric */
+       BLI_flagArcs(graph, 0);
+       
+       /* mark down all nodes as not on the symmetry axis */
+       BLI_flagNodes(graph, 0);
+
+       node = root_node;
+       
+       /* sanity check REMOVE ME */
+       if (node->degree > 0)
+       {
+               arc = node->arcs[0];
+               
+               if (node->degree == 1)
+               {
+                       markdownSymmetryArc(graph, arc, node, 1, limit);
+               }
+               else
+               {
+                       markdownSymmetryArc(graph, NULL, node, 1, limit);
+               }
+               
+
+
+               /* mark down non-symetric arcs */
+               for (arc = graph->arcs.first; arc; arc = arc->next)
+               {
+                       if (arc->symmetry_level < 0)
+                       {
+                               arc->symmetry_level = 0;
+                       }
+                       else
+                       {
+                               /* mark down nodes with the lowest level symmetry axis */
+                               if (arc->head->symmetry_level == 0 || arc->head->symmetry_level > arc->symmetry_level)
+                               {
+                                       arc->head->symmetry_level = arc->symmetry_level;
+                               }
+                               if (arc->tail->symmetry_level == 0 || arc->tail->symmetry_level > arc->symmetry_level)
+                               {
+                                       arc->tail->symmetry_level = arc->symmetry_level;
+                               }
+                       }
+               }
+       }
+}
+
index 24112c7f11afd2664d69804428f2748f19aaa5e5..c3660f700b3a9d62ca748d6a349f7a1e6b5ce61b 100644 (file)
@@ -68,6 +68,8 @@ typedef struct EditBone
 
 } EditBone;
 
+void   make_boneList(struct ListBase *list, struct ListBase *bones, EditBone *parent);
+void   editbones_to_armature (struct ListBase *list, struct Object *ob);
 
 void   adduplicate_armature(void);
 void   addvert_armature(void);
@@ -142,6 +144,9 @@ void        show_all_armature_bones(void);
 
 #define BONESEL_NOSEL  0x80000000      /* Indicates a negative number */
 
+/* from autoarmature */
+void BIF_retargetArmature();
+
 #endif
 
 
index c0542e3f34c17ca11a4bffaa1b89eabdf3258dd6..4207047d24e831a2307db88f06606ccae472cd9e 100644 (file)
@@ -446,7 +446,8 @@ void curvemap_buttons(struct uiBlock *block, struct CurveMapping *cumap, char la
 #define B_SETMCOL_RND          2083
 #define B_DRAWBWEIGHTS         2084
 
-#define B_GEN_SKELETON         2090
+#define B_GEN_SKELETON         2085
+#define B_RETARGET_SKELETON    2086
 
 /* *********************** */
 #define B_VGROUPBUTS           2100
index c8352aedec54e63710c5fba10cff01ea4e8b0982..6baaaf1ca963ff99f57734d96656c2d1ad019be9 100644 (file)
 
 #include "DNA_listBase.h"
 
+#include "BLI_graph.h"
+
+struct GHash;
 struct EdgeHash;
 struct ReebArc;
 struct ReebEdge;
 struct ReebNode;
 
 typedef struct ReebGraph {
-       ListBase arcs;
-       ListBase nodes;
+       ListBase        arcs;
+       ListBase        nodes;
+       
+       float length;
+       
+       FreeArc                 free_arc;
+       FreeNode                free_node;
+       RadialSymmetry  radial_symmetry;
+       AxialSymmetry   axial_symmetry;
+       /*********************************/
+       
+       int resolution;
        int totnodes;
        struct EdgeHash *emap;
+       struct ReebGraph *link_up; /* for multi resolution filtering, points to higher levels */
 } ReebGraph;
 
 typedef struct EmbedBucket {
@@ -49,13 +63,22 @@ typedef struct EmbedBucket {
 } EmbedBucket;
 
 typedef struct ReebNode {
-       struct ReebNode *next, *prev;
+       void *next, *prev;
+       float p[3];
+       int flag;
+
+       int degree;
        struct ReebArc **arcs;
+
+       int symmetry_level;
+       int symmetry_flag;
+       float symmetry_axis[3];
+       /*********************************/
+
        int index;
-       int degree;
        float weight;
-       float p[3];
-       int flags;
+       struct ReebNode *link_down; /* for multi resolution filtering, points to lower levels, if present */
+       struct ReebNode *link_up;
 } ReebNode;
 
 typedef struct ReebEdge {
@@ -63,15 +86,28 @@ typedef struct ReebEdge {
        struct ReebArc  *arc;
        struct ReebNode *v1, *v2;
        struct ReebEdge *nextEdge;
+       int flag;
 } ReebEdge;
 
 typedef struct ReebArc {
-       struct ReebArc *next, *prev;
+       void *next, *prev;
+       struct ReebNode *head, *tail;
+       int flag;
+
+       float length;
+
+       int symmetry_level;
+       int symmetry_group;
+       int symmetry_flag;
+       /*********************************/
+
        ListBase edges;
-       struct ReebNode *v1, *v2;
+       int bcount;
        struct EmbedBucket *buckets;
-       int     bcount;
-       int flags;
+
+       struct GHash *faces;    
+       float angle;
+       struct ReebArc *link_up; /* for multi resolution filtering, points to higher levels */
 } ReebArc;
 
 typedef struct ReebArcIterator {
@@ -80,6 +116,7 @@ typedef struct ReebArcIterator {
        int start;
        int end;
        int stride;
+       int length;
 } ReebArcIterator;
 
 struct EditMesh;
@@ -87,21 +124,27 @@ struct EditMesh;
 int weightToHarmonic(struct EditMesh *em);
 int weightFromDistance(struct EditMesh *em);
 int weightFromLoc(struct EditMesh *me, int axis);
-void weightToVCol(struct EditMesh *em);
+void weightToVCol(struct EditMesh *em, int index);
+void arcToVCol(struct ReebGraph *rg, struct EditMesh *em, int index);
+void angleToVCol(struct EditMesh *em, int index);
 void renormalizeWeight(struct EditMesh *em, float newmax);
 
 ReebGraph * generateReebGraph(struct EditMesh *me, int subdivisions);
-void freeGraph(ReebGraph *rg);
-void exportGraph(ReebGraph *rg, int count);
-
-#define OTHER_NODE(arc, node) ((arc->v1 == node) ? arc->v2 : arc->v1)
+ReebGraph * newReebGraph();
 
 void initArcIterator(struct ReebArcIterator *iter, struct ReebArc *arc, struct ReebNode *head);
 void initArcIterator2(struct ReebArcIterator *iter, struct ReebArc *arc, int start, int end);
+void initArcIteratorStart(struct ReebArcIterator *iter, struct ReebArc *arc, struct ReebNode *head, int start);
 struct EmbedBucket * nextBucket(struct ReebArcIterator *iter);
+struct EmbedBucket * nextNBucket(ReebArcIterator *iter, int n);
+struct EmbedBucket * peekBucket(ReebArcIterator *iter, int n);
+struct EmbedBucket * currentBucket(struct ReebArcIterator *iter);
+struct EmbedBucket * previousBucket(struct ReebArcIterator *iter);
+int iteratorStopped(struct ReebArcIterator *iter);
 
 /* Filtering */
 void filterNullReebGraph(ReebGraph *rg);
+int filterSmartReebGraph(ReebGraph *rg, float threshold);
 int filterExternalReebGraph(ReebGraph *rg, float threshold);
 int filterInternalReebGraph(ReebGraph *rg, float threshold);
 
@@ -110,18 +153,27 @@ void repositionNodes(ReebGraph *rg);
 void postprocessGraph(ReebGraph *rg, char mode);
 void removeNormalNodes(ReebGraph *rg);
 
-/* Graph processing */
-void buildAdjacencyList(ReebGraph *rg);
-
 void sortNodes(ReebGraph *rg);
 void sortArcs(ReebGraph *rg);
 
-int subtreeDepth(ReebNode *node, ReebArc *rootArc);
-int countConnectedArcs(ReebGraph *rg, ReebNode *node);
-int hasAdjacencyList(ReebGraph *rg); 
-int    isGraphCyclic(ReebGraph *rg);
-
-/* Sanity check */
+/*------------ Sanity check ------------*/
 void verifyBuckets(ReebGraph *rg);
+void verifyFaces(ReebGraph *rg);
+
+/*********************** PUBLIC *********************************/
+ReebGraph *BIF_ReebGraphFromEditMesh(void);
+ReebGraph *BIF_ReebGraphMultiFromEditMesh(void);
+void BIF_flagMultiArcs(ReebGraph *rg, int flag);
+
+void BIF_GlobalReebGraphFromEditMesh(void);
+void BIF_GlobalReebFree(void);
+
+ReebNode *BIF_otherNodeFromIndex(ReebArc *arc, ReebNode *node);
+ReebNode *BIF_lowestLevelNode(ReebNode *node);
+
+void REEB_freeGraph(ReebGraph *rg);
+void REEB_exportGraph(ReebGraph *rg, int count);
+void REEB_draw();
+
 
 #endif /*REEB_H_*/
index 75affbfa7f5133e7854fe540c031cc5e68838ce2..c235cfe2b316347aa85f16141f12fc374d90c914 100644 (file)
@@ -433,14 +433,19 @@ typedef struct ToolSettings {
        float skgen_angle_limit;
        float skgen_correlation_limit;
        float skgen_symmetry_limit;
+       float skgen_retarget_angle_weight;
+       float skgen_retarget_length_weight;
+       float skgen_retarget_distance_weight;
        short skgen_options;
        char  skgen_postpro;
        char  skgen_postpro_passes;
        char  skgen_subdivisions[3];
+       char  skgen_multi_level;
+       
+       char tpad[7];
        
        /* Alt+RMB option */
        char edge_mode;
-       char pad3[4];
 } ToolSettings;
 
 /* Used by all brushes to store their properties, which can be directly set
@@ -830,12 +835,19 @@ typedef struct Scene {
 #define RETOPO_ELLIPSE 4
 
 /* toolsettings->skgen_options */
-#define SKGEN_FILTER_INTERNAL  1
-#define SKGEN_FILTER_EXTERNAL  2
-#define        SKGEN_SYMMETRY                  4
-#define        SKGEN_CUT_LENGTH                8
-#define        SKGEN_CUT_ANGLE                 16
-#define        SKGEN_CUT_CORRELATION   32
+#define SKGEN_FILTER_INTERNAL  (1 << 0)
+#define SKGEN_FILTER_EXTERNAL  (1 << 1)
+#define        SKGEN_SYMMETRY                  (1 << 2)
+#define        SKGEN_CUT_LENGTH                (1 << 3)
+#define        SKGEN_CUT_ANGLE                 (1 << 4)
+#define        SKGEN_CUT_CORRELATION   (1 << 5)
+#define        SKGEN_HARMONIC                  (1 << 6)
+#define        SKGEN_STICK_TO_EMBEDDING        (1 << 7)
+#define        SKGEN_ADAPTIVE_DISTANCE         (1 << 8)
+#define SKGEN_FILTER_SMART             (1 << 9)
+#define SKGEN_DISP_LENGTH              (1 << 10)
+#define SKGEN_DISP_WEIGHT              (1 << 11)
+#define SKGEN_DISP_ORIG                        (1 << 12)
 
 #define        SKGEN_SUB_LENGTH                0
 #define        SKGEN_SUB_ANGLE                 1
diff --git a/source/blender/src/autoarmature.c b/source/blender/src/autoarmature.c
new file mode 100644 (file)
index 0000000..1aab8ce
--- /dev/null
@@ -0,0 +1,1647 @@
+/**
+ * $Id:
+ *
+ * ***** BEGIN GPL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software Foundation,
+ * Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+ *
+ * Contributor(s): Martin Poirier
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ * autoarmature.c: Interface for automagically manipulating armature (retarget, created, ...)
+ */
+
+#include <ctype.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h> 
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "MEM_guardedalloc.h"
+
+#include "DNA_ID.h"
+#include "DNA_armature_types.h"
+#include "DNA_mesh_types.h"
+#include "DNA_meshdata_types.h"
+#include "DNA_object_types.h"
+#include "DNA_scene_types.h"
+#include "DNA_view3d_types.h"
+
+#include "BLI_blenlib.h"
+#include "BLI_arithb.h"
+#include "BLI_editVert.h"
+#include "BLI_ghash.h"
+#include "BLI_graph.h"
+#include "BLI_rand.h"
+
+#include "BDR_editobject.h"
+
+#include "BKE_global.h"
+#include "BKE_utildefines.h"
+
+#include "BIF_editarmature.h"
+#include "BIF_space.h"
+
+#include "PIL_time.h"
+
+#include "mydevice.h"
+#include "reeb.h" // FIX ME
+#include "blendef.h"
+
+/************ RIG RETARGET DATA STRUCTURES ***************/
+
+struct RigJoint;
+struct RigGraph;
+struct RigNode;
+struct RigArc;
+struct RigEdge;
+
+typedef struct RigGraph {
+       ListBase        arcs;
+       ListBase        nodes;
+       
+       float length;
+       
+       FreeArc                 free_arc;
+       FreeNode                free_node;
+       RadialSymmetry  radial_symmetry;
+       AxialSymmetry   axial_symmetry;
+       /*********************************/
+
+       struct RigNode *head;
+       ReebGraph *link_mesh;
+} RigGraph;
+
+typedef struct RigNode {
+       void *next, *prev;
+       float p[3];
+       int flag;
+
+       int degree;
+       struct BArc **arcs;
+
+       int symmetry_level;
+       int symmetry_flag;
+       float symmetry_axis[3];
+       /*********************************/
+
+       ReebNode *link_mesh;
+} RigNode;
+
+typedef struct RigArc {
+       void *next, *prev;
+       RigNode *head, *tail;
+       int flag;
+
+       float length;
+
+       int symmetry_level;
+       int symmetry_group;
+       int symmetry_flag;
+       /*********************************/
+       
+       ListBase edges;
+       int count;
+       ReebArc *link_mesh;
+} RigArc;
+
+typedef struct RigEdge {
+       struct RigEdge *next, *prev;
+       float head[3], tail[3];
+       float length;
+       float angle;
+       EditBone *bone;
+} RigEdge;
+
+/*******************************************************************************************************/
+
+static void RIG_calculateEdgeAngle(RigEdge *edge_first, RigEdge *edge_second);
+
+/*********************************** EDITBONE UTILS ****************************************************/
+
+int countEditBoneChildren(ListBase *list, EditBone *parent)
+{
+       EditBone *ebone;
+       int count = 0;
+       
+       for (ebone = list->first; ebone; ebone = ebone->next)
+       {
+               if (ebone->parent == parent)
+               {
+                       count++;
+               }
+       }
+       
+       return count;
+}
+
+EditBone* nextEditBoneChild(ListBase *list, EditBone *parent, int n)
+{
+       EditBone *ebone;
+       
+       for (ebone = list->first; ebone; ebone = ebone->next)
+       {
+               if (ebone->parent == parent)
+               {
+                       if (n == 0)
+                       {
+                               return ebone;
+                       }
+                       n--;
+               }
+       }
+       
+       return NULL;
+}
+
+
+/************************************ DESTRUCTORS ******************************************************/
+
+void RIG_freeRigArc(BArc *arc)
+{
+       BLI_freelistN(&((RigArc*)arc)->edges);
+}
+
+void RIG_freeRigGraph(BGraph *rg)
+{
+       BNode *node;
+       BArc *arc;
+       
+       for (arc = rg->arcs.first; arc; arc = arc->next)
+       {
+               RIG_freeRigArc(arc);
+       }
+       BLI_freelistN(&rg->arcs);
+       
+       for (node = rg->nodes.first; node; node = node->next)
+       {
+               BLI_freeNode((BGraph*)rg, (BNode*)node);
+       }
+       BLI_freelistN(&rg->nodes);
+       
+       MEM_freeN(rg);
+}
+
+/************************************* ALLOCATORS ******************************************************/
+
+static RigGraph *newRigGraph()
+{
+       RigGraph *rg;
+       rg = MEM_callocN(sizeof(RigGraph), "rig graph");
+       
+       rg->head = NULL;
+       
+       rg->free_arc = RIG_freeRigArc;
+       rg->free_node = NULL;
+       
+       return rg;
+}
+
+static RigArc *newRigArc(RigGraph *rg)
+{
+       RigArc *arc;
+       
+       arc = MEM_callocN(sizeof(RigArc), "rig arc");
+       arc->count = 0;
+       
+       BLI_addtail(&rg->arcs, arc);
+       
+       return arc;
+}
+
+static RigNode *newRigNodeHead(RigGraph *rg, RigArc *arc, float p[3])
+{
+       RigNode *node;
+       node = MEM_callocN(sizeof(RigNode), "rig node");
+       BLI_addtail(&rg->nodes, node);
+
+       VECCOPY(node->p, p);
+       node->degree = 1;
+       node->arcs = NULL;
+       
+       arc->head = node;
+       
+       return node;
+}
+
+static void addRigNodeHead(RigGraph *rg, RigArc *arc, RigNode *node)
+{
+       node->degree++;
+
+       arc->head = node;
+}
+
+static RigNode *newRigNodeTail(RigGraph *rg, RigArc *arc, float p[3])
+{
+       RigNode *node;
+       node = MEM_callocN(sizeof(RigNode), "rig node");
+       BLI_addtail(&rg->nodes, node);
+
+       VECCOPY(node->p, p);
+       node->degree = 1;
+       node->arcs = NULL;
+       
+       arc->tail = node;
+
+       return node;
+}
+
+static void RIG_addEdgeToArc(RigArc *arc, float tail[3], EditBone *bone)
+{
+       RigEdge *edge;
+
+       edge = MEM_callocN(sizeof(RigEdge), "rig edge");
+       BLI_addtail(&arc->edges, edge);
+
+
+       VECCOPY(edge->tail, tail);
+       edge->bone = bone;
+       
+       if (edge->prev == NULL)
+       {
+               VECCOPY(edge->head, arc->head->p);
+       }
+       else
+       {
+               RigEdge *last_edge = edge->prev;
+               VECCOPY(edge->head, last_edge->tail);
+               RIG_calculateEdgeAngle(last_edge, edge);
+       }
+       
+       edge->length = VecLenf(edge->head, edge->tail);
+       
+       arc->length += edge->length;
+       
+       arc->count += 1;
+}
+
+
+/*******************************************************************************************************/
+
+static void RIG_calculateEdgeAngle(RigEdge *edge_first, RigEdge *edge_second)
+{
+       float vec_first[3], vec_second[3];
+       
+       VecSubf(vec_first, edge_first->tail, edge_first->head); 
+       VecSubf(vec_second, edge_second->tail, edge_second->head);
+
+       Normalize(vec_first);
+       Normalize(vec_second);
+       
+       edge_first->angle = saacos(Inpf(vec_first, vec_second));
+}
+
+/*******************************************************************************************************/
+
+static void RIG_arcFromBoneChain(RigGraph *rg, ListBase *list, EditBone *root_bone, RigNode *starting_node)
+{
+       EditBone *bone, *last_bone = NULL;
+       RigArc *arc;
+       int contain_head = 0;
+       
+       arc = newRigArc(rg);
+       
+       if (starting_node == NULL)
+       {
+               starting_node = newRigNodeHead(rg, arc, root_bone->head);
+       }
+       else
+       {
+               addRigNodeHead(rg, arc, starting_node);
+       }
+       
+       for(bone = root_bone; bone; bone = nextEditBoneChild(list, bone, 0))
+       {
+               int nb_children;
+               
+               if (bone->parent && (bone->flag & BONE_CONNECTED) == 0)
+               {
+                       RIG_addEdgeToArc(arc, bone->head, NULL);
+               }
+               
+               RIG_addEdgeToArc(arc, bone->tail, bone);
+               
+               if (strcmp(bone->name, "head") == 0)
+               {
+                       contain_head = 1;
+               }
+               
+               nb_children = countEditBoneChildren(list, bone);
+               if (nb_children > 1)
+               {
+                       RigNode *end_node = newRigNodeTail(rg, arc, bone->tail);
+                       int i;
+                       
+                       for (i = 0; i < nb_children; i++)
+                       {
+                               root_bone = nextEditBoneChild(list, bone, i);
+                               RIG_arcFromBoneChain(rg, list, root_bone, end_node);
+                       }
+                       
+                       /* arc ends here, break */
+                       break;
+               }
+               last_bone = bone;
+       }
+       
+       /* If the loop exited without forking */
+       if (bone == NULL)
+       {
+               newRigNodeTail(rg, arc, last_bone->tail);
+       }
+
+       if (contain_head)
+       {
+               rg->head = arc->tail;
+       }
+}
+
+/*******************************************************************************************************/
+static void RIG_findHead(RigGraph *rg)
+{
+       if (rg->head == NULL)
+       {
+               if (BLI_countlist(&rg->arcs) == 1)
+               {
+                       RigArc *arc = rg->arcs.first;
+                       
+                       rg->head = (RigNode*)arc->head;
+               }
+               else
+               {
+                       RigArc *arc;
+                       
+                       for (arc = rg->arcs.first; arc; arc = arc->next)
+                       {
+                               RigEdge *edge = arc->edges.last;
+                               
+                               if (edge->bone->flag & (BONE_TIPSEL|BONE_SELECTED))
+                               {
+                                       rg->head = arc->tail;
+                                       break;
+                               }
+                       }
+               }
+               
+               if (rg->head == NULL)
+               {
+                       rg->head = rg->nodes.first;
+               }
+       }
+}
+
+/*******************************************************************************************************/
+
+void RIG_printNode(RigNode *node, char name[])
+{
+       printf("%s %p %i <%0.3f, %0.3f, %0.3f>\n", name, node, node->degree, node->p[0], node->p[1], node->p[2]);
+       
+       if (node->symmetry_flag & SYM_TOPOLOGICAL)
+       {
+               if (node->symmetry_flag & SYM_AXIAL)
+                       printf("Symmetry AXIAL\n");
+               else if (node->symmetry_flag & SYM_RADIAL)
+                       printf("Symmetry RADIAL\n");
+                       
+               printvecf("symmetry axis", node->symmetry_axis);
+       }
+}
+
+void RIG_printArcBones(RigArc *arc)
+{
+       RigEdge *edge;
+
+       for (edge = arc->edges.first; edge; edge = edge->next)
+       {
+               if (edge->bone)
+                       printf("%s ", edge->bone->name);
+               else
+                       printf("---- ");
+       }
+       printf("\n");
+}
+
+void RIG_printArc(RigArc *arc)
+{
+       RigEdge *edge;
+
+       printf("\n");
+
+       RIG_printNode((RigNode*)arc->head, "head");
+
+       for (edge = arc->edges.first; edge; edge = edge->next)
+       {
+               printf("\tinner joints %0.3f %0.3f %0.3f\n", edge->tail[0], edge->tail[1], edge->tail[2]);
+               printf("\t\tlength %f\n", edge->length);
+               printf("\t\tangle %f\n", edge->angle * 180 / M_PI);
+               if (edge->bone)
+                       printf("\t\t%s\n", edge->bone->name);
+       }       
+       printf("symmetry level: %i\n", arc->symmetry_level);
+
+       RIG_printNode((RigNode*)arc->tail, "tail");
+}
+
+void RIG_printGraph(RigGraph *rg)
+{
+       RigArc *arc;
+
+       for (arc = rg->arcs.first; arc; arc = arc->next)
+       {
+               RIG_printArc(arc);      
+       }
+       
+       if (rg->head)
+       {
+               RIG_printNode(rg->head, "HEAD NODE:");
+       }
+       else
+       {
+               printf("HEAD NODE: NONE\n");
+       }       
+}
+
+/*******************************************************************************************************/
+
+static RigGraph *armatureToGraph(ListBase *list)
+{
+       EditBone *ebone;
+       RigGraph *rg;
+       
+       rg = newRigGraph();
+
+       /* Do the rotations */
+       for (ebone = list->first; ebone; ebone=ebone->next){
+               if (ebone->parent == NULL)
+               {
+                       RIG_arcFromBoneChain(rg, list, ebone, NULL);
+               }
+       }
+       
+       BLI_removeDoubleNodes((BGraph*)rg, 0.001);
+       
+       BLI_buildAdjacencyList((BGraph*)rg);
+       
+       RIG_findHead(rg);
+       
+       return rg;
+}
+
+/************************************ RETARGETTING *****************************************************/
+
+typedef enum 
+{
+       RETARGET_LENGTH,
+       RETARGET_AGGRESSIVE
+} RetargetMode; 
+
+static RetargetMode detectArcRetargetMode(RigArc *arc);
+static void retargetArctoArcLength(RigArc *iarc);
+
+
+static RetargetMode detectArcRetargetMode(RigArc *iarc)
+{
+       RetargetMode mode = RETARGET_AGGRESSIVE;
+       ReebArc *earc = iarc->link_mesh;
+       RigEdge *edge;
+       int large_angle = 0;
+       float avg_angle = 0;
+       float avg_length = 0;
+       int nb_edges = 0;
+       
+       
+       for (edge = iarc->edges.first; edge; edge = edge->next)
+       {
+               avg_angle += edge->angle;
+               nb_edges++;
+       }
+       
+       avg_angle /= nb_edges - 1; /* -1 because last edge doesn't have an angle */
+
+       avg_length = iarc->length / nb_edges;
+       
+       
+       if (nb_edges > 2)
+       {
+               for (edge = iarc->edges.first; edge; edge = edge->next)
+               {
+                       if (fabs(edge->angle - avg_angle) > M_PI / 6)
+                       {
+                               large_angle = 1;
+                       }
+               }
+       }
+       else if (nb_edges == 2 && avg_angle > 0)
+       {
+               large_angle = 1;
+       }
+               
+       
+       if (large_angle == 0)
+       {
+               mode = RETARGET_LENGTH;
+       }
+       
+       if (earc->bcount <= (iarc->count - 1))
+       {
+               mode = RETARGET_LENGTH;
+       }
+       
+       mode = RETARGET_AGGRESSIVE;
+       
+       return mode;
+}
+
+static void printPositions(int *positions, int nb_positions)
+{
+       int i;
+       
+       for (i = 0; i < nb_positions; i++)
+       {
+               printf("%i ", positions[i]);
+       }
+       printf("\n");
+}
+
+static float costDistance(ReebArcIterator *iter, float *vec0, float *vec1, int i0, int i1)
+{
+       EmbedBucket *bucket = NULL;
+       float max_dist = 0;
+       float v1[3], v2[3], c[3];
+       float v1_inpf;
+
+       if (G.scene->toolsettings->skgen_retarget_distance_weight > 0)
+       {
+               VecSubf(v1, vec0, vec1);
+               
+               v1_inpf = Inpf(v1, v1);
+               
+               if (v1_inpf > 0)
+               {
+                       int j;
+                       for (j = i0 + 1; j < i1 - 1; j++)
+                       {
+                               float dist;
+                               
+                               bucket = peekBucket(iter, j);
+       
+                               VecSubf(v2, bucket->p, vec1);
+               
+                               Crossf(c, v1, v2);
+                               
+                               dist = Inpf(c, c) / v1_inpf;
+                               
+                               max_dist = dist > max_dist ? dist : max_dist;
+                       }
+                       
+                       return max_dist;
+               }
+               else
+               {
+                       return FLT_MAX;
+               }
+               
+               return G.scene->toolsettings->skgen_retarget_distance_weight * max_dist;
+       }
+       else
+       {
+               return 0;
+       }
+}
+
+static float costAngle(float original_angle, float current_angle)
+{
+       return 0;
+}
+
+static float costLength(float original_length, float current_length)
+{
+       float length_ratio = fabs((current_length - original_length) / original_length);
+       return G.scene->toolsettings->skgen_retarget_length_weight * length_ratio * length_ratio;
+}
+
+static float calcCost(ReebArcIterator *iter, RigEdge *e1, RigEdge *e2, float *vec0, float *vec1, float *vec2, int i0, int i1, int i2)
+{
+       float vec_second[3], vec_first[3];
+       float angle = e1->angle;
+       float test_angle, length1, length2;
+       float new_cost = 0;
+
+       VecSubf(vec_second, vec2, vec1);
+       length2 = Normalize(vec_second);
+
+       VecSubf(vec_first, vec1, vec0); 
+       length1 = Normalize(vec_first);
+       
+       if (length1 > 0 && length2 > 0)
+       {
+               test_angle = saacos(Inpf(vec_first, vec_second));
+               /* ANGLE COST HERE */
+               if (angle > 0)
+               {
+                       new_cost += G.scene->toolsettings->skgen_retarget_angle_weight * fabs((test_angle - angle) / angle);
+               }
+               else
+               {
+                       new_cost += G.scene->toolsettings->skgen_retarget_angle_weight * fabs(test_angle);
+               }
+       }
+       else
+       {
+               new_cost += M_PI;
+       }
+
+       /* Length cost */
+       new_cost += costLength(e1->length, length1);
+       new_cost += costLength(e2->length, length2);
+
+       /* Distance cost */
+       new_cost += costDistance(iter, vec0, vec1, i0, i1);
+       new_cost += costDistance(iter, vec1, vec2, i1, i2);
+
+       return new_cost;
+}
+
+#define MAX_COST 100 /* FIX ME */
+
+static void calcGradient(RigEdge *e1, RigEdge *e2, ReebArcIterator *iter, int index, int nb_joints, float *cost_cube, int *positions, float **vec_cache)
+{
+       EmbedBucket *bucket = NULL;
+       float *vec0, *vec1, *vec2;
+       float current_cost;
+       int i0, i1, i2;
+       int next_position;
+
+       vec0 = vec_cache[index];
+       vec1 = vec_cache[index + 1];
+       vec2 = vec_cache[index + 2];
+       
+       if (index == 0)
+       {
+               i0 = iter->start;
+       }
+       else
+       {
+               i0 = positions[index - 1];
+       }
+       
+       i1 = positions[index];
+       
+       if (index +1 == nb_joints)
+       {
+               i2 = iter->end;
+       }
+       else
+       {
+               i2 = positions[index + 1];
+       }
+
+
+       current_cost = calcCost(iter, e1, e2, vec0, vec1, vec2, i0, i1, i2);
+       cost_cube[index * 3 + 1] = current_cost;
+       
+       next_position = positions[index] + 1;
+       
+       if (index + 1 < nb_joints && next_position == positions[index + 1])
+       {
+               cost_cube[index * 3 + 2] = MAX_COST;
+       }
+       else
+       {
+               bucket = peekBucket(iter, next_position);
+               
+               if (bucket == NULL)
+               {
+                       cost_cube[index * 3 + 2] = MAX_COST;
+               }
+               else
+               {
+                       vec1 = bucket->p;
+                       
+                       cost_cube[index * 3 + 2] = calcCost(iter, e1, e2, vec0, vec1, vec2, i0, next_position, i2) - current_cost;
+               }
+       }
+
+       next_position = positions[index] - 1;
+       
+       if (index - 1 > -1 && next_position == positions[index - 1])
+       {
+               cost_cube[index * 3] = MAX_COST;
+       }
+       else
+       {
+               bucket = peekBucket(iter, next_position);
+               
+               if (bucket == NULL)
+               {
+                       cost_cube[index * 3] = MAX_COST;
+               }
+               else
+               {
+                       vec1 = bucket->p;
+                       
+                       cost_cube[index * 3] = calcCost(iter, e1, e2, vec0, vec1, vec2, i0, next_position, i2) - current_cost;
+               }
+       }
+}
+
+static float probability(float delta_cost, float iterations)
+{
+       if (delta_cost < 0)
+       {
+               return 1;
+       }
+       else
+       {
+               float temperature = (1 - iterations);
+               return (float)exp(delta_cost) * temperature;
+       }
+}
+
+static int neighbour(int nb_joints, float *cost_cube, int *moving_joint, int *moving_direction)
+{
+       int total = 0;
+       int chosen = 0;
+       int i;
+       
+       for (i = 0; i < nb_joints; i++)
+       {
+               if (cost_cube[i * 3] < MAX_COST)
+               {
+                       total++;
+               }
+               
+               if (cost_cube[i * 3 + 2] < MAX_COST)
+               {
+                       total++;
+               }
+       }
+       
+       if (total == 0)
+       {
+               return 0;
+       }
+       
+       chosen = (int)(BLI_drand() * total);
+       
+       for (i = 0; i < nb_joints; i++)
+       {
+               if (cost_cube[i * 3] < MAX_COST)
+               {
+                       if (chosen == 0)
+                       {
+                               *moving_joint = i;
+                               *moving_direction = -1;
+                               break;
+                       }
+                       chosen--;
+               }
+               
+               if (cost_cube[i * 3 + 2] < MAX_COST)
+               {
+                       if (chosen == 0)
+                       {
+                               *moving_joint = i;
+                               *moving_direction = 1;
+                               break;
+                       }
+                       chosen--;
+               }
+       }
+       
+       return 1;
+}
+
+static void retargetArctoArcAggresive(RigArc *iarc)
+{
+       ReebArcIterator iter;
+       RigEdge *edge;
+       EmbedBucket *bucket = NULL;
+       ReebNode *node_start, *node_end;
+       ReebArc *earc = iarc->link_mesh;
+       float min_cost = FLT_MAX;
+       float *vec0, *vec1, *vec2;
+       float **vec_cache;
+       float *cost_cache;
+       int *best_positions;
+       int *positions;
+       int nb_edges = BLI_countlist(&iarc->edges);
+       int nb_joints = nb_edges - 1;
+       int symmetry_axis = 0;
+       int last_index = 0;
+       int first_pass = 1;
+       int must_move = nb_joints - 1;
+       int i;
+       
+       printf("aggressive\n");
+
+       positions = MEM_callocN(sizeof(int) * nb_joints, "Aggresive positions");
+       best_positions = MEM_callocN(sizeof(int) * nb_joints, "Best Aggresive positions");
+       cost_cache = MEM_callocN(sizeof(float) * nb_edges, "Cost cache");
+       vec_cache = MEM_callocN(sizeof(float*) * (nb_edges + 1), "Vec cache");
+       
+       /* symmetry axis */
+       if (earc->symmetry_level == 1 && iarc->symmetry_level == 1)
+       {
+               symmetry_axis = 1;
+               node_start = earc->tail;
+               node_end = earc->head;
+       }
+       else
+       {
+               node_start = earc->head;
+               node_end = earc->tail;
+       }
+       
+       /* init with first values */
+       for (i = 0; i < nb_joints; i++)
+       {
+               positions[i] = i + 1;
+               //positions[i] = (earc->bcount / nb_edges) * (i + 1);
+       }
+       
+       /* init cost cache */
+       for (i = 0; i < nb_edges; i++)
+       {
+               cost_cache[i] = 0;
+       }
+       
+       vec_cache[0] = node_start->p;
+       vec_cache[nb_edges] = node_end->p;
+
+#if 0 /* BRUTE FORCE */
+       while(1)
+       {
+               float cost = 0;
+               int need_calc = 0;
+               
+               /* increment to next possible solution */
+               
+               i = nb_joints - 1;
+
+               /* increment positions, starting from the last one
+                * until a valid increment is found
+                * */
+               for (i = must_move; i >= 0; i--)
+               {
+                       int remaining_joints = nb_joints - (i + 1); 
+                       
+                       positions[i] += 1;
+                       need_calc = i;
+                       
+                       if (positions[i] + remaining_joints < earc->bcount)
+                       {
+                               break;
+                       }
+               }
+               
+               if (first_pass)
+               {
+                       need_calc = 0;
+                       first_pass = 0;
+               }
+
+               if (i == -1)
+               {
+                       break;
+               }
+               
+               /* reset joints following the last increment*/
+               for (i = i + 1; i < nb_joints; i++)
+               {
+                       positions[i] = positions[i - 1] + 1;
+               }
+       
+               /* calculating cost */
+               initArcIterator(&iter, earc, node_start);
+               
+               vec0 = NULL;
+               vec1 = node_start->p;
+               vec2 = NULL;
+               
+               for (edge = iarc->edges.first, i = 0, last_index = 0;
+                        edge;
+                        edge = edge->next, i += 1)
+               {
+
+                       if (i >= need_calc)
+                       { 
+                               float vec_first[3], vec_second[3];
+                               float length1, length2;
+                               float new_cost = 0;
+                               int i1, i2;
+                               
+                               if (i < nb_joints)
+                               {
+                                       i2 = positions[i];
+                                       bucket = peekBucket(&iter, positions[i]);
+                                       vec2 = bucket->p;
+                                       vec_cache[i + 1] = vec2; /* update cache for updated position */
+                               }
+                               else
+                               {
+                                       i2 = iter.length;
+                                       vec2 = node_end->p;
+                               }
+                               
+                               if (i > 0)
+                               {
+                                       i1 = positions[i - 1];
+                               }
+                               else
+                               {
+                                       i1 = 1;
+                               }
+                               
+                               vec1 = vec_cache[i];
+                               
+
+                               VecSubf(vec_second, vec2, vec1);
+                               length2 = Normalize(vec_second);
+       
+                               /* check angle */
+                               if (i != 0 && G.scene->toolsettings->skgen_retarget_angle_weight > 0)
+                               {
+                                       RigEdge *previous = edge->prev;
+                                       float angle = previous->angle;
+                                       float test_angle;
+                                       
+                                       vec0 = vec_cache[i - 1];
+                                       VecSubf(vec_first, vec1, vec0); 
+                                       length1 = Normalize(vec_first);
+                                       
+                                       if (length1 > 0 && length2 > 0)
+                                       {
+                                               test_angle = saacos(Inpf(vec_first, vec_second));
+                                               /* ANGLE COST HERE */
+                                               if (angle > 0)
+                                               {
+                                                       new_cost += G.scene->toolsettings->skgen_retarget_angle_weight * fabs((test_angle - angle) / angle);
+                                               }
+                                               else
+                                               {
+                                                       new_cost += G.scene->toolsettings->skgen_retarget_angle_weight * fabs(test_angle);
+                                               }
+                                       }
+                                       else
+                                       {
+                                               new_cost += G.scene->toolsettings->skgen_retarget_angle_weight;
+                                       }
+                               }
+       
+                               /* Length Cost */
+                               new_cost += costLength(edge->length, length2);
+                               
+                               /* Distance Cost */
+                               new_cost += calcMaximumDistance(&iter, vec1, vec2, i1, i2);
+                               
+                               cost_cache[i] = new_cost;
+                       }
+                       
+                       cost += cost_cache[i];
+                       
+                       if (cost > min_cost)
+                       {
+                               must_move = i;
+                               break;
+                       }
+               }
+               
+               if (must_move != i || must_move > nb_joints - 1)
+               {
+                       must_move = nb_joints - 1;
+               }
+
+               /* cost optimizing */
+               if (cost < min_cost)
+               {
+                       min_cost = cost;
+                       memcpy(best_positions, positions, sizeof(int) * nb_joints);
+               }
+       }
+#elif 1 /* SIMULATED ANNEALING */
+       {
+               RigEdge *previous;
+               float *cost_cube;
+               int k, kmax;
+               
+               kmax = 10000;
+               
+               BLI_srand(nb_joints);
+               
+               /* [joint: index][position: -1, 0, +1] */
+               cost_cube = MEM_callocN(sizeof(float) * 3 * nb_joints, "Cost Cube");
+               
+               initArcIterator(&iter, earc, node_start);
+
+               /* init vec_cache */
+               for (i = 0; i < nb_joints; i++)
+               {
+                       bucket = peekBucket(&iter, positions[i]);
+                       vec_cache[i + 1] = bucket->p;
+               }
+               
+               min_cost = 0;
+
+               /* init cost cube */
+               for (previous = iarc->edges.first, edge = previous->next, i = 0;
+                        edge;
+                        previous = edge, edge = edge->next, i += 1)
+               {
+                       calcGradient(previous, edge, &iter, i, nb_joints, cost_cube, positions, vec_cache);
+                       
+                       min_cost += cost_cube[3 * i + 1];
+               }
+               
+               for (k = 0; k < kmax; k++)
+               {
+                       int status;
+                       int moving_joint = -1;
+                       int move_direction = -1;
+                       float delta_cost;
+                       
+                       status = neighbour(nb_joints, cost_cube, &moving_joint, &move_direction);
+                       
+                       if (status == 0)
+                       {
+                               break;
+                       }
+                       
+                       delta_cost = cost_cube[moving_joint * 3 + (1 + move_direction)];
+
+                       if (probability(delta_cost, (float)k / (float)kmax) > BLI_frand())
+                       {
+                               /* update position */                   
+                               positions[moving_joint] += move_direction;
+                               
+                               /* update vector cache */
+                               bucket = peekBucket(&iter, positions[moving_joint]);
+                               vec_cache[moving_joint + 1] = bucket->p;
+                               
+                               min_cost += delta_cost;
+
+                               printf("%i: %0.3f\n", k, delta_cost);
+       
+                               /* update cost cube */                  
+                               for (previous = iarc->edges.first, edge = previous->next, i = 0;
+                                        edge;
+                                        previous = edge, edge = edge->next, i += 1)
+                               {
+                                       if (i == moving_joint - 1 ||
+                                               i == moving_joint ||
+                                               i == moving_joint + 1)
+                                       {
+                                               calcGradient(previous, edge, &iter, i, nb_joints, cost_cube, positions, vec_cache);
+                                       }
+                               }
+                       }
+               }
+               
+               memcpy(best_positions, positions, sizeof(int) * nb_joints);
+               
+               MEM_freeN(cost_cube);
+       }       
+#else    /* GRADIENT DESCENT*/
+       {
+               RigEdge *previous;
+               float *cost_cube;
+               
+               /* [joint: index][position: -1, 0, +1] */
+               cost_cube = MEM_callocN(sizeof(float) * 3 * nb_joints, "Cost Cube");
+               
+               initArcIterator(&iter, earc, node_start);
+
+               /* init vec_cache */
+               for (i = 0; i < nb_joints; i++)
+               {
+                       bucket = peekBucket(&iter, positions[i]);
+                       vec_cache[i + 1] = bucket->p;
+               }
+
+               /* init cost cube */
+               for (previous = iarc->edges.first, edge = previous->next, i = 0;
+                        edge;
+                        previous = edge, edge = edge->next, i += 1)
+               {
+                       calcGradient(previous, edge, &iter, i, nb_joints, cost_cube, positions, vec_cache);
+               }
+
+               while(1)
+               {
+                       float min_gradient = 0;
+                       int moving_joint = -1;
+                       int move_direction = -1;
+               
+                       printf("-----------------\n");
+               
+                       for (i = 0; i < nb_joints; i++)
+                       {
+                               printf("%i[%i]: %f\t\t(%f)\t\t%f\n", i, positions[i], cost_cube[i * 3], cost_cube[i * 3 + 1], cost_cube[i * 3 + 2]); 
+                               if (cost_cube[i * 3] < min_gradient)
+                               {
+                                       min_gradient = cost_cube[i * 3]; 
+                                       moving_joint = i;
+                                       move_direction = -1;
+                               }
+                               
+                               if  (cost_cube[i * 3 + 2] < min_gradient)
+                               {
+                                       min_gradient = cost_cube[i * 3 + 2]; 
+                                       moving_joint = i;
+                                       move_direction = 1;
+                               }
+                       }
+                       
+                       if (moving_joint == -1)
+                       {
+                               break;
+                       }
+                       
+                       positions[moving_joint] += move_direction;
+                       
+                       /* update vector cache */
+                       bucket = peekBucket(&iter, positions[moving_joint]);
+                       vec_cache[moving_joint + 1] = bucket->p;
+
+                       /* update cost cube */                  
+                       for (previous = iarc->edges.first, edge = previous->next, i = 0;
+                                edge;
+                                previous = edge, edge = edge->next, i += 1)
+                       {
+                               if (i == moving_joint - 1 ||
+                                       i == moving_joint ||
+                                       i == moving_joint + 1)
+                               {
+                                       calcGradient(previous, edge, &iter, i, nb_joints, cost_cube, positions, vec_cache);
+                               }
+                       }
+                       
+                       
+               }
+               
+               memcpy(best_positions, positions, sizeof(int) * nb_joints);
+               
+               MEM_freeN(cost_cube);
+       }
+#endif
+
+       vec0 = node_start->p;
+       initArcIterator(&iter, earc, node_start);
+       
+       printPositions(best_positions, nb_joints);
+       printf("min_cost %f\n", min_cost);
+       printf("buckets: %i\n", earc->bcount);
+
+       /* set joints to best position */
+       for (edge = iarc->edges.first, i = 0, last_index = 0;
+                edge;
+                edge = edge->next, i++)
+       {
+               EditBone *bone = edge->bone;
+               
+               if (i < nb_joints)
+               {
+                       bucket = peekBucket(&iter, best_positions[i]);
+                       vec1 = bucket->p;
+               }
+               else
+               {
+                       vec1 = node_end->p;
+               }
+               
+               if (bone)
+               {
+                       VECCOPY(bone->head, vec0);
+                       VECCOPY(bone->tail, vec1);
+                       printf("===\n");
+                       printvecf("vec0", vec0);
+                       printvecf("vec1", vec1);
+                       if (i < nb_joints)
+                               printf("position: %i\n", best_positions[i]);
+                       printf("last_index: %i\n", last_index);
+               }
+               
+               vec0 = vec1;
+       }
+       
+       MEM_freeN(positions);
+       MEM_freeN(best_positions);
+       MEM_freeN(cost_cache);
+       MEM_freeN(vec_cache);
+}
+
+static void retargetArctoArcLength(RigArc *iarc)
+{
+       ReebArcIterator iter;
+       ReebArc *earc = iarc->link_mesh;
+       ReebNode *node_start, *node_end;
+       RigEdge *edge;
+       EmbedBucket *bucket = NULL;
+       float embedding_length = 0;
+       float *vec0 = NULL;
+       float *vec1 = NULL;
+       float *previous_vec = NULL;
+       int symmetry_axis = 0;
+
+       
+       /* symmetry axis */
+       if (earc->symmetry_level == 1 && iarc->symmetry_level == 1)
+       {
+               symmetry_axis = 1;
+               node_start = (ReebNode*)earc->tail;
+               node_end = (ReebNode*)earc->head;
+       }
+       else
+       {
+               node_start = (ReebNode*)earc->head;
+               node_end = (ReebNode*)earc->tail;
+       }
+       
+       initArcIterator(&iter, earc, node_start);
+
+       bucket = nextBucket(&iter);
+       
+       vec0 = node_start->p;
+       
+       while (bucket != NULL)
+       {
+               vec1 = bucket->p;
+               
+               embedding_length += VecLenf(vec0, vec1);
+               
+               vec0 = vec1;
+               bucket = nextBucket(&iter);
+       }
+       
+       embedding_length += VecLenf(node_end->p, vec1);
+       
+       /* fit bones */
+       initArcIterator(&iter, earc, node_start);
+
+       bucket = nextBucket(&iter);
+
+       vec0 = node_start->p;
+       previous_vec = vec0;
+       vec1 = bucket->p;
+       
+       printf("arc: %f embedding %f\n",  iarc->length, embedding_length);
+       
+       for (edge = iarc->edges.first; edge; edge = edge->next)
+       {
+               EditBone *bone = edge->bone;
+               float new_bone_length = edge->length / iarc->length * embedding_length;
+
+#if 0          
+               while (bucket && new_bone_length > VecLenf(vec0, vec1))
+               {
+                       bucket = nextBucket(&iter);
+                       previous_vec = vec1;
+                       vec1 = bucket->p;
+               }
+               
+               if (bucket == NULL)
+               {
+                       vec1 = node_end->p;
+               }
+               
+               if (embedding_length < VecLenf(vec0, vec1))
+               {
+                       float dv[3], off[3];
+                       float a, b, c, f;
+                       
+                       /* Solve quadratic distance equation */
+                       VecSubf(dv, vec1, previous_vec);
+                       a = Inpf(dv, dv);
+                       
+                       VecSubf(off, previous_vec, vec0);
+                       b = 2 * Inpf(dv, off);
+                       
+                       c = Inpf(off, off) - (new_bone_length * new_bone_length);
+                       
+                       f = (-b + (float)sqrt(b * b - 4 * a * c)) / (2 * a);
+                       
+                       if (isnan(f) == 0 && f < 1.0f)
+                       {
+                               VECCOPY(vec1, dv);
+                               VecMulf(vec1, f);
+                               VecAddf(vec1,vec1, vec0);
+                       }
+               }
+#else
+               float length = 0;
+
+               while (bucket && new_bone_length > length)
+               {
+                       length += VecLenf(previous_vec, vec1);
+                       bucket = nextBucket(&iter);
+                       previous_vec = vec1;
+                       vec1 = bucket->p;
+               }
+               
+               if (bucket == NULL)
+               {
+                       vec1 = node_end->p;
+               }
+#endif
+
+               /* no need to move virtual edges (space between unconnected bones) */           
+               if (bone)
+               {
+                       printf("BONE: %s\n", bone->name);
+                       VECCOPY(bone->head, vec0);
+                       VECCOPY(bone->tail, vec1);
+               }
+               printvecf("vec0", vec0);
+               printvecf("vec1", vec1);
+               printf("old: %f target: %f new: %f\n", edge->length, new_bone_length, VecLenf(vec0, vec1));
+               
+               vec0 = vec1;
+               previous_vec = vec1;
+       }
+}
+
+static void retargetArctoArc(RigArc *iarc)
+{
+       ReebArc *earc = iarc->link_mesh;
+       
+       if (BLI_countlist(&iarc->edges) == 1)
+       {
+               RigEdge *edge = iarc->edges.first;
+               EditBone *bone = edge->bone;
+               
+               /* symmetry axis */
+               if (earc->symmetry_level == 1 && iarc->symmetry_level == 1)
+               {
+                       VECCOPY(bone->head, earc->tail->p);
+                       VECCOPY(bone->tail, earc->head->p);
+               }
+               /* or not */
+               else
+               {
+                       VECCOPY(bone->head, earc->head->p);
+                       VECCOPY(bone->tail, earc->tail->p);
+               }
+       }
+       else
+       {
+               RetargetMode mode = detectArcRetargetMode(iarc);
+               
+               if (mode == RETARGET_AGGRESSIVE)
+               {
+                       retargetArctoArcAggresive(iarc);
+               }
+               else
+               {               
+                       retargetArctoArcLength(iarc);
+               }
+       }
+}
+
+static void matchMultiResolutionArc(RigNode *start_node, RigArc *next_iarc, ReebArc *next_earc)
+{
+       ReebNode *enode = next_earc->head;
+       int ishape, eshape;
+       int MAGIC_NUMBER = 100; /* FIXME */
+
+       ishape = BLI_subtreeShape((BNode*)start_node, (BArc*)next_iarc, 1) % MAGIC_NUMBER;
+       eshape = BLI_subtreeShape((BNode*)enode, (BArc*)next_earc, 1) % MAGIC_NUMBER;
+       
+       while (ishape != eshape && next_earc->link_up)
+       {
+               next_earc->flag = 1; // mark previous as taken, to prevent backtrack on lower levels
+               
+               next_earc = next_earc->link_up;
+               enode = next_earc->head;
+               eshape = BLI_subtreeShape((BNode*)enode, (BArc*)next_earc, 1) % MAGIC_NUMBER;
+       } 
+
+       next_earc->flag = 1; // mark as taken
+       next_iarc->link_mesh = next_earc;
+}
+
+static void matchMultiResolutionStartingArc(ReebGraph *reebg, RigArc *iarc, RigNode *inode)
+{
+       ReebArc *earc;
+       ReebNode *enode;
+       int ishape, eshape;
+       int MAGIC_NUMBER = 100; /* FIXME */
+       
+       earc = reebg->arcs.first;
+       enode = earc->head;
+       
+       ishape = BLI_subtreeShape((BNode*)inode, (BArc*)iarc, 1) % MAGIC_NUMBER;
+       eshape = BLI_subtreeShape((BNode*)enode, (BArc*)earc, 1) % MAGIC_NUMBER;
+       
+       while (ishape != eshape && reebg->link_up)
+       {
+               earc->flag = 1; // mark previous as taken, to prevent backtrack on lower levels
+               
+               reebg = reebg->link_up;
+               
+               earc = reebg->arcs.first;
+               enode = earc->head;
+               
+               eshape = BLI_subtreeShape((BNode*)enode, (BArc*)earc, 1) % MAGIC_NUMBER;
+       } 
+
+       earc->flag = 1; // mark as taken
+       iarc->link_mesh = earc;
+       inode->link_mesh = enode;
+}
+
+static void findCorrespondingArc(RigArc *start_arc, RigNode *start_node, RigArc *next_iarc)
+{
+       ReebNode *enode = start_node->link_mesh;
+       ReebArc *next_earc;
+       int symmetry_level = next_iarc->symmetry_level;
+       int symmetry_group = next_iarc->symmetry_group;
+       int symmetry_flag = next_iarc->symmetry_flag;
+       int i;
+       
+       next_iarc->link_mesh = NULL;
+               
+       for(i = 0; i < enode->degree; i++)
+       {
+               next_earc = (ReebArc*)enode->arcs[i];
+               if (next_earc->flag == 0 && /* not already taken */
+                       next_earc->symmetry_flag == symmetry_flag &&
+                       next_earc->symmetry_group == symmetry_group &&
+                       next_earc->symmetry_level == symmetry_level)
+               {
+                       printf("-----------------------\n");
+                       printf("CORRESPONDING ARC FOUND\n");
+                       RIG_printArcBones(next_iarc);
+                       printf("flag %i -- symmetry level %i -- symmetry flag %i\n", next_earc->flag, next_earc->symmetry_level, next_earc->symmetry_flag);
+                       
+                       matchMultiResolutionArc(start_node, next_iarc, next_earc);
+                       break;
+               }
+       }
+       
+       /* not found, try at higher nodes (lower node might have filtered internal arcs, messing shape of tree */
+       if (next_iarc->link_mesh == NULL)
+       {
+               if (enode->link_up)
+               {
+                       start_node->link_mesh = enode->link_up;
+                       findCorrespondingArc(start_arc, start_node, next_iarc);
+               }
+       }
+
+       /* still not found, print debug info */
+       if (next_iarc->link_mesh == NULL)
+       {
+               printf("--------------------------\n");
+               printf("NO CORRESPONDING ARC FOUND\n");
+               RIG_printArcBones(next_iarc);
+               
+               printf("LOOKING FOR\n");
+               printf("flag %i -- symmetry level %i -- symmetry flag %i\n", 0, symmetry_level, symmetry_flag);
+               
+               printf("CANDIDATES\n");
+               for(i = 0; i < enode->degree; i++)
+               {
+                       next_earc = (ReebArc*)enode->arcs[i];
+                       printf("flag %i -- symmetry level %i -- symmetry flag %i\n", next_earc->flag, next_earc->symmetry_level, next_earc->symmetry_flag);
+               }
+       }
+
+}
+
+static void retargetSubgraph(RigGraph *rigg, RigArc *start_arc, RigNode *start_node)
+{
+       RigArc *iarc = start_arc;
+       ReebArc *earc = start_arc->link_mesh;
+       RigNode *inode = start_node;
+       ReebNode *enode = start_node->link_mesh;
+       int i;
+               
+       retargetArctoArc(iarc);
+       
+       enode = BIF_otherNodeFromIndex(earc, enode);
+       inode = (RigNode*)BLI_otherNode((BArc*)iarc, (BNode*)inode);
+       
+       /* Link with lowest possible node
+        * Enabling going back to lower levels for each arc
+        * */
+       inode->link_mesh = BIF_lowestLevelNode(enode);
+       
+       for(i = 0; i < inode->degree; i++)
+       {
+               RigArc *next_iarc = (RigArc*)inode->arcs[i];
+               
+               /* no back tracking */
+               if (next_iarc != iarc)
+               {
+                       findCorrespondingArc(iarc, inode, next_iarc);
+                       if (next_iarc->link_mesh)
+                       {
+                               retargetSubgraph(rigg, next_iarc, inode);
+                       }
+               }
+       }
+}
+
+static void retargetGraphs(RigGraph *rigg)
+{
+       ReebGraph *reebg = rigg->link_mesh;
+       ReebArc *earc;
+       RigArc *iarc;
+       ReebNode *enode;
+       RigNode *inode;
+       
+       /* flag all ReebArcs as not taken */
+       BIF_flagMultiArcs(reebg, 0);
+       
+       /* return to first level */
+       reebg = rigg->link_mesh;
+       
+       iarc = (RigArc*)rigg->head->arcs[0];
+       inode = iarc->tail;
+       
+       matchMultiResolutionStartingArc(reebg, iarc, inode);
+
+       earc = iarc->link_mesh; /* has been set earlier */
+       enode = earc->head;
+
+       inode->link_mesh = enode;
+
+       retargetSubgraph(rigg, iarc, inode);
+}
+
+void BIF_retargetArmature()
+{
+       Object *ob;
+       Base *base;
+       ReebGraph *reebg;
+       
+       reebg = BIF_ReebGraphMultiFromEditMesh();
+       
+       
+       printf("Reeb Graph created\n");
+
+       base= FIRSTBASE;
+       for (base = FIRSTBASE; base; base = base->next)
+       {
+               if TESTBASELIB(base) {
+                       ob = base->object;
+
+                       if (ob->type==OB_ARMATURE)
+                       {
+                               RigGraph *rigg;
+                               ListBase  list;
+                               bArmature *arm;
+                               
+                               arm = ob->data;
+                       
+                               /* Put the armature into editmode */
+                               list.first= list.last = NULL;
+                               make_boneList(&list, &arm->bonebase, NULL);
+                       
+                               rigg = armatureToGraph(&list);
+                               
+                               BLI_markdownSymmetry((BGraph*)rigg, (BNode*)rigg->head, G.scene->toolsettings->skgen_symmetry_limit);
+                               
+                               printf("Armature graph created\n");
+               
+                               RIG_printGraph(rigg);
+                               
+                               rigg->link_mesh = reebg;
+                               
+                               printf("retargetting %s\n", ob->id.name);
+                               
+                               retargetGraphs(rigg);
+                               
+                               /* Turn the list into an armature */
+                               editbones_to_armature(&list, ob);
+                               
+                               BLI_freelistN(&list);
+
+                               RIG_freeRigGraph((BGraph*)rigg);
+                       }
+               }
+       }
+
+       REEB_freeGraph(reebg);
+       
+       BIF_undo_push("Retarget Skeleton");
+       
+       exit_editmode(EM_FREEDATA|EM_FREEUNDO|EM_WAITCURSOR); // freedata, and undo
+
+       allqueue(REDRAWVIEW3D, 0);
+}
index 529a795a101e90c7167b07838b04ae01e5d1e48f..4df4826de28f4a3d3f58945f0d90154becabc2d8 100644 (file)
 #include "butspace.h" // own module
 #include "multires.h"
 
+#include "reeb.h"
+
 static float editbutweight= 1.0;
 float editbutvweight= 1;
 static int actmcol= 0, acttface= 0, acttface_rnd = 0, actmcol_rnd = 0;
@@ -4894,6 +4896,9 @@ void do_meshbuts(unsigned short event)
        case B_GEN_SKELETON:
                generateSkeleton();
                break;
+       case B_RETARGET_SKELETON:
+               BIF_retargetArmature();
+               break;
        }
 
        /* WATCH IT: previous events only in editmode! */
@@ -4992,6 +4997,80 @@ static void skgen_reorder(void *option, void *arg2)
        }
 }
 
+static void skgen_graphgen(void *arg1, void *arg2)
+{
+       BIF_GlobalReebGraphFromEditMesh();
+       allqueue(REDRAWVIEW3D, 0);
+}
+
+static void skgen_graphfree(void *arg1, void *arg2)
+{
+       BIF_GlobalReebFree();
+       allqueue(REDRAWVIEW3D, 0);
+}
+
+static void skgen_graph_block(uiBlock *block)
+{
+       uiBlockBeginAlign(block);
+       uiDefButS(block, NUM, B_DIFF, "Resolution:",                                                    1025,150,225,19, &G.scene->toolsettings->skgen_resolution,10.0,1000.0, 0, 0,            "Specifies the resolution of the graph's embedding");
+       uiDefButBitS(block, TOG, SKGEN_HARMONIC, B_DIFF,                "H",                    1250,150, 25,19, &G.scene->toolsettings->skgen_options, 0, 0, 0, 0,                                     "Apply harmonic smoothing to the weighting");
+       uiDefButBitS(block, TOG, SKGEN_FILTER_INTERNAL, B_DIFF, "Filter In",    1025,130, 83,19, &G.scene->toolsettings->skgen_options, 0, 0, 0, 0,                                     "Filter internal small arcs from graph");
+       uiDefButF(block, NUM, B_DIFF,                                                   "",                             1111,130,164,19, &G.scene->toolsettings->skgen_threshold_internal,0.0, 10.0, 10, 0,     "Specify the threshold ratio for filtering internal arcs");
+       uiDefButBitS(block, TOG, SKGEN_FILTER_EXTERNAL, B_DIFF, "Filter Ex",    1025,110, 53,19, &G.scene->toolsettings->skgen_options, 0, 0, 0, 0,                                     "Filter external small arcs from graph");
+       uiDefButBitS(block, TOG, SKGEN_FILTER_SMART,    B_DIFF, "Sm",                   1078,110, 30,19, &G.scene->toolsettings->skgen_options, 0, 0, 0, 0,                                     "Smart Filtering");
+       uiDefButF(block, NUM, B_DIFF,                                                   "",                             1111,110,164,19, &G.scene->toolsettings->skgen_threshold_external,0.0, 10.0, 10, 0,     "Specify the threshold ratio for filtering external arcs");
+       uiBlockEndAlign(block);
+}
+
+static void editing_panel_mesh_skgen_display(Object *ob, Mesh *me)
+{
+       uiBlock *block;
+       uiBut *but;
+
+       block= uiNewBlock(&curarea->uiblocks, "editing_panel_mesh_skgen_display", UI_EMBOSS, UI_HELV, curarea->win);
+       uiNewPanelTabbed("Mesh Tools More", "Skgen");
+       if(uiNewPanel(curarea, block, "Graph", "Editing", 960, 0, 318, 204)==0) return;
+       
+       but = uiDefBut(block, BUT, B_DIFF, "Generate",                          1025,170,125,19, 0, 0, 0, 0, 0, "Generate Graph from Mesh");
+       uiButSetFunc(but, skgen_graphgen, NULL, NULL);
+       but = uiDefBut(block, BUT, B_DIFF, "Free",                                      1150,170,125,19, 0, 0, 0, 0, 0, "Free Graph from Mesh");
+       uiButSetFunc(but, skgen_graphfree, NULL, NULL);
+       
+       skgen_graph_block(block);
+
+       uiDefButBitS(block, TOG, SKGEN_DISP_LENGTH, REDRAWVIEW3D,       "Length",                       1025, 60, 83,19, &G.scene->toolsettings->skgen_options, 0, 0, 0, 0,             "Show Length");
+       uiDefButBitS(block, TOG, SKGEN_DISP_WEIGHT, REDRAWVIEW3D,       "Weight",                       1108, 60, 83,19, &G.scene->toolsettings->skgen_options, 0, 0, 0, 0,             "Show Weight");
+       uiDefButBitS(block, TOG, SKGEN_DISP_ORIG, REDRAWVIEW3D,         "Original",                     1191, 60, 84,19, &G.scene->toolsettings->skgen_options, 0, 0, 0, 0,             "Show Original Graph");
+
+       uiDefButC(block, NUM, REDRAWVIEW3D,                                             "Level:",                       1025, 40, 125,19, &G.scene->toolsettings->skgen_multi_level, 0, 5, 1, 0,"Specify the level to draw");
+}
+
+static void editing_panel_mesh_skgen_retarget(Object *ob, Mesh *me)
+{
+       uiBlock *block;
+
+       block= uiNewBlock(&curarea->uiblocks, "editing_panel_mesh_skgen_retarget", UI_EMBOSS, UI_HELV, curarea->win);
+       uiNewPanelTabbed("Mesh Tools More", "Skgen");
+       if(uiNewPanel(curarea, block, "Retarget", "Editing", 960, 0, 318, 204)==0) return;
+       
+       uiDefBut(block, BUT, B_RETARGET_SKELETON, "Retarget Skeleton",                  1025,170,250,19, 0, 0, 0, 0, 0, "Retarget Selected Armature to this Mesh");
+
+       skgen_graph_block(block);
+
+       uiDefButF(block, NUM, B_DIFF,                                                   "Ang:",                 1025, 60, 83,19, &G.scene->toolsettings->skgen_retarget_angle_weight, 0, 10, 1, 0,              "Angle Weight");
+       uiDefButF(block, NUM, B_DIFF,                                                   "Len:",                 1108, 60, 83,19, &G.scene->toolsettings->skgen_retarget_length_weight, 0, 10, 1, 0,             "Length Weight");
+       uiDefButF(block, NUM, B_DIFF,                                                   "Dist:",                1191, 60, 84,19, &G.scene->toolsettings->skgen_retarget_distance_weight, 0, 10, 1, 0,           "Distance Weight");
+
+       uiBlockBeginAlign(block);
+       uiDefButBitS(block, TOG, SKGEN_SYMMETRY, B_DIFF,                "Symmetry",             1025, 30,125,19, &G.scene->toolsettings->skgen_options, 0, 0, 0, 0,                                     "Restore symmetries based on topology");
+       uiDefButF(block, NUM, B_DIFF,                                                   "T:",                   1150, 30,125,19, &G.scene->toolsettings->skgen_symmetry_limit,0.0, 1.0, 10, 0,  "Specify the threshold distance for considering potential symmetric arcs");
+       uiDefButC(block, NUM, B_DIFF,                                                   "P:",                   1025, 10, 62,19, &G.scene->toolsettings->skgen_postpro_passes, 0, 10, 10, 0,            "Specify the number of processing passes on the embeddings");
+       uiDefButC(block, ROW, B_DIFF,                                                   "Smooth",               1087, 10, 63,19, &G.scene->toolsettings->skgen_postpro, 5.0, (float)SKGEN_SMOOTH, 0, 0, "Smooth embeddings");
+       uiDefButC(block, ROW, B_DIFF,                                                   "Average",              1150, 10, 62,19, &G.scene->toolsettings->skgen_postpro, 5.0, (float)SKGEN_AVERAGE, 0, 0, "Average embeddings");
+       uiDefButC(block, ROW, B_DIFF,                                                   "Sharpen",              1212, 10, 63,19, &G.scene->toolsettings->skgen_postpro, 5.0, (float)SKGEN_SHARPEN, 0, 0, "Sharpen embeddings");
+       uiBlockEndAlign(block);
+}
+
 static void editing_panel_mesh_skgen(Object *ob, Mesh *me)
 {
        uiBlock *block;
@@ -4999,17 +5078,14 @@ static void editing_panel_mesh_skgen(Object *ob, Mesh *me)
        int i;
 
        block= uiNewBlock(&curarea->uiblocks, "editing_panel_mesh_skgen", UI_EMBOSS, UI_HELV, curarea->win);
-       if(uiNewPanel(curarea, block, "Skeleton Generator", "Editing", 960, 0, 318, 204)==0) return;
+       uiNewPanelTabbed("Mesh Tools More", "Skgen");
+       if(uiNewPanel(curarea, block, "Generator", "Editing", 960, 0, 318, 204)==0) return;
        
-       uiDefBut(block, BUT, B_GEN_SKELETON, "Generate Skeleton",                       1025,170,250,19, 0, 0, 0, 0, 0, "Generate Skeleton from Mesh");
+       uiDefBut(block, BUT, B_GEN_SKELETON, "Generate",                        1025,170,250,19, 0, 0, 0, 0, 0, "Generate Skeleton from Mesh");
 
-       uiBlockBeginAlign(block);
-       uiDefButS(block, NUM, B_DIFF, "Resolution:",                                                    1025,150,250,19, &G.scene->toolsettings->skgen_resolution,10.0,1000.0, 0, 0,            "Specifies the resolution of the graph's embedding");
-       uiDefButBitS(block, TOG, SKGEN_FILTER_INTERNAL, B_DIFF, "Filter In",    1025,130, 83,19, &G.scene->toolsettings->skgen_options, 0, 0, 0, 0,                                     "Filter internal small arcs from graph");
-       uiDefButF(block, NUM, B_DIFF,                                                   "T:",                   1111,130,164,19, &G.scene->toolsettings->skgen_threshold_internal,0.0, 1.0, 10, 0,      "Specify the threshold ratio for filtering internal arcs");
-       uiDefButBitS(block, TOG, SKGEN_FILTER_EXTERNAL, B_DIFF, "Filter Ex",    1025,110, 83,19, &G.scene->toolsettings->skgen_options, 0, 0, 0, 0,                                     "Filter external small arcs from graph");
-       uiDefButF(block, NUM, B_DIFF,                                                   "T:",                   1111,110,164,19, &G.scene->toolsettings->skgen_threshold_external,0.0, 1.0, 10, 0,      "Specify the threshold ratio for filtering external arcs");
+       skgen_graph_block(block);
 
+       uiBlockBeginAlign(block);
        for(i = 0; i < SKGEN_SUB_TOTAL; i++)
        {
                int y = 90 - 20 * i;
@@ -5029,8 +5105,10 @@ static void editing_panel_mesh_skgen(Object *ob, Mesh *me)
                                uiDefButF(block, NUM, B_DIFF,                                                   "T:",                   1111, y,164,19, &G.scene->toolsettings->skgen_angle_limit,0.0, 90.0, 10, 0,                     "Specify the threshold angle in degrees for subdivision");
                                break;
                        case SKGEN_SUB_CORRELATION:
-                               uiDefButBitS(block, TOG, SKGEN_CUT_CORRELATION, B_DIFF, "Correlation",  1041, y, 67,19, &G.scene->toolsettings->skgen_options, 0, 0, 0, 0,                                      "Subdivide arcs based on correlation");
-                               uiDefButF(block, NUM, B_DIFF,                                                   "T:",                   1111, y,164,19, &G.scene->toolsettings->skgen_correlation_limit,0.0, 1.0, 0.01, 0,      "Specify the threshold correlation for subdivision");
+                               uiDefButBitS(block, TOG, SKGEN_CUT_CORRELATION, B_DIFF, "Adaptative",   1041, y, 67,19, &G.scene->toolsettings->skgen_options, 0, 0, 0, 0,                                      "Subdivide arcs adaptatively");
+                               uiDefButF(block, NUM, B_DIFF,                                                   "T:",                   1111, y,114,19, &G.scene->toolsettings->skgen_correlation_limit,0.0, 1.0, 0.01, 0,      "Specify the adaptive threshold for subdivision");
+                               uiDefButBitS(block, TOG, SKGEN_STICK_TO_EMBEDDING, B_DIFF,              "E",                    1225, y, 25,19, &G.scene->toolsettings->skgen_options, 0, 0, 0, 0,                                      "Stick endpoint to embedding");
+                               uiDefButBitS(block, TOG, SKGEN_ADAPTIVE_DISTANCE, B_DIFF,               "D",                    1250, y, 25,19, &G.scene->toolsettings->skgen_options, 0, 0, 0, 0,                                      "Adaptive distance (on) or variance(off)");
                                break;
                }
        }
@@ -6471,8 +6549,9 @@ void editing_panels()
                        editing_panel_mesh_tools1(ob, ob->data);
                        uiNewPanelTabbed("Mesh Tools 1", "Editing");
                        
-                       if (G.rt == 42) /* hidden for now, no time for docs */
-                               editing_panel_mesh_skgen(ob, ob->data);
+                       editing_panel_mesh_skgen(ob, ob->data);
+                       editing_panel_mesh_skgen_retarget(ob, ob->data);
+                       editing_panel_mesh_skgen_display(ob, ob->data);
                        
                        editing_panel_mesh_uvautocalculation();
                        if (EM_texFaceCheck())
index f595a101f63e3484a069430faa44f09bf2de62f8..780e7847bfbdbd0b68ed8d0a2f2fc19094982b59 100644 (file)
 
 #include "RE_pipeline.h"       // make_stars
 
+#include "reeb.h"
+
 #include "multires.h"
 
 /* For MULTISAMPLE_ARB #define.
@@ -3147,6 +3149,8 @@ void drawview3dspace(ScrArea *sa, void *spacedata)
                        BIF_drawPropCircle(); // only editmode and particles have proportional edit
                BIF_drawSnap();
        }
+       
+       REEB_draw();
 
        if(G.scene->radio) RAD_drawall(v3d->drawtype>=OB_SOLID);
        
index 6310dd0a262c8e6a85bd0139af08c004d19b36f3..d2557d9050f7068b1c3a89c8e08ba3733e9c2ff7 100644 (file)
@@ -4189,542 +4189,7 @@ void transform_armature_mirror_update(void)
 /*************************************** SKELETON GENERATOR ******************************************/
 /*****************************************************************************************************/
 
-/**************************************** SYMMETRY HANDLING ******************************************/
 
-void markdownSymmetryArc(ReebArc *arc, ReebNode *node, int level);
-
-void mirrorAlongAxis(float v[3], float center[3], float axis[3])
-{
-       float dv[3], pv[3];
-       
-       VecSubf(dv, v, center);
-       Projf(pv, dv, axis);
-       VecMulf(pv, -2);
-       VecAddf(v, v, pv);
-}
-
-/* Helper structure for radial symmetry */
-typedef struct RadialArc
-{
-       ReebArc *arc; 
-       float n[3]; /* normalized vector joining the nodes of the arc */
-} RadialArc;
-
-void reestablishRadialSymmetry(ReebNode *node, int depth, float axis[3])
-{
-       RadialArc *ring = NULL;
-       RadialArc *unit;
-       float limit = G.scene->toolsettings->skgen_symmetry_limit;
-       int symmetric = 1;
-       int count = 0;
-       int i;
-
-       /* count the number of arcs in the symmetry ring */
-       for (i = 0; node->arcs[i] != NULL; i++)
-       {
-               ReebArc *connectedArc = node->arcs[i];
-               
-               /* depth is store as a negative in flag. symmetry level is positive */
-               if (connectedArc->flags == -depth)
-               {
-                       count++;
-               }
-       }
-
-       ring = MEM_callocN(sizeof(RadialArc) * count, "radial symmetry ring");
-       unit = ring;
-
-       /* fill in the ring */
-       for (unit = ring, i = 0; node->arcs[i] != NULL; i++)
-       {
-               ReebArc *connectedArc = node->arcs[i];
-               
-               /* depth is store as a negative in flag. symmetry level is positive */
-               if (connectedArc->flags == -depth)
-               {
-                       ReebNode *otherNode = OTHER_NODE(connectedArc, node);
-                       float vec[3];
-
-                       unit->arc = connectedArc;
-
-                       /* project the node to node vector on the symmetry plane */
-                       VecSubf(unit->n, otherNode->p, node->p);
-                       Projf(vec, unit->n, axis);
-                       VecSubf(unit->n, unit->n, vec);
-
-                       Normalize(unit->n);
-
-                       unit++;
-               }
-       }
-
-       /* sort ring */
-       for (i = 0; i < count - 1; i++)
-       {
-               float minAngle = 3; /* arbitrary high value, higher than 2, at least */
-               int minIndex = -1;
-               int j;
-
-               for (j = i + 1; j < count; j++)
-               {
-                       float angle = Inpf(ring[i].n, ring[j].n);
-
-                       /* map negative values to 1..2 */
-                       if (angle < 0)
-                       {
-                               angle = 1 - angle;
-                       }
-
-                       if (angle < minAngle)
-                       {
-                               minIndex = j;
-                               minAngle = angle;
-                       }
-               }
-
-               /* swap if needed */
-               if (minIndex != i + 1)
-               {
-                       RadialArc tmp;
-                       tmp = ring[i + 1];
-                       ring[i + 1] = ring[minIndex];
-                       ring[minIndex] = tmp;
-               }
-       }
-
-       for (i = 0; i < count && symmetric; i++)
-       {
-               ReebNode *node1, *node2;
-               float tangent[3];
-               float normal[3];
-               float p[3];
-               int j = (i + 1) % count; /* next arc in the circular list */
-
-               VecAddf(tangent, ring[i].n, ring[j].n);
-               Crossf(normal, tangent, axis);
-               
-               node1 = OTHER_NODE(ring[i].arc, node);
-               node2 = OTHER_NODE(ring[j].arc, node);
-
-               VECCOPY(p, node2->p);
-               mirrorAlongAxis(p, node->p, normal);
-               
-               /* check if it's within limit before continuing */
-               if (VecLenf(node1->p, p) > limit)
-               {
-                       symmetric = 0;
-               }
-
-       }
-
-       if (symmetric)
-       {
-               /* first pass, merge incrementally */
-               for (i = 0; i < count - 1; i++)
-               {
-                       ReebNode *node1, *node2;
-                       float tangent[3];
-                       float normal[3];
-                       int j = i + 1;
-       
-                       VecAddf(tangent, ring[i].n, ring[j].n);
-                       Crossf(normal, tangent, axis);
-                       
-                       node1 = OTHER_NODE(ring[i].arc, node);
-                       node2 = OTHER_NODE(ring[j].arc, node);
-       
-                       /* mirror first node and mix with the second */
-                       mirrorAlongAxis(node1->p, node->p, normal);
-                       VecLerpf(node2->p, node2->p, node1->p, 1.0f / (j + 1));
-                       
-                       /* Merge buckets
-                        * there shouldn't be any null arcs here, but just to be safe 
-                        * */
-                       if (ring[i].arc->bcount > 0 && ring[j].arc->bcount > 0)
-                       {
-                               ReebArcIterator iter1, iter2;
-                               EmbedBucket *bucket1 = NULL, *bucket2 = NULL;
-                               
-                               initArcIterator(&iter1, ring[i].arc, node);
-                               initArcIterator(&iter2, ring[j].arc, node);
-                               
-                               bucket1 = nextBucket(&iter1);
-                               bucket2 = nextBucket(&iter2);
-                       
-                               /* Make sure they both start at the same value */       
-                               while(bucket1 && bucket1->val < bucket2->val)
-                               {
-                                       bucket1 = nextBucket(&iter1);
-                               }
-                               
-                               while(bucket2 && bucket2->val < bucket1->val)
-                               {
-                                       bucket2 = nextBucket(&iter2);
-                               }
-               
-               
-                               for ( ;bucket1 && bucket2; bucket1 = nextBucket(&iter1), bucket2 = nextBucket(&iter2))
-                               {
-                                       bucket2->nv += bucket1->nv; /* add counts */
-                                       
-                                       /* mirror on axis */
-                                       mirrorAlongAxis(bucket1->p, node->p, normal);
-                                       /* add bucket2 in bucket1 */
-                                       VecLerpf(bucket2->p, bucket2->p, bucket1->p, (float)bucket1->nv / (float)(bucket2->nv));
-                               }
-                       }
-               }
-               
-               /* second pass, mirror back on previous arcs */
-               for (i = count - 1; i > 0; i--)
-               {
-                       ReebNode *node1, *node2;
-                       float tangent[3];
-                       float normal[3];
-                       int j = i - 1;
-       
-                       VecAddf(tangent, ring[i].n, ring[j].n);
-                       Crossf(normal, tangent, axis);
-                       
-                       node1 = OTHER_NODE(ring[i].arc, node);
-                       node2 = OTHER_NODE(ring[j].arc, node);
-       
-                       /* copy first node than mirror */
-                       VECCOPY(node2->p, node1->p);
-                       mirrorAlongAxis(node2->p, node->p, normal);
-                       
-                       /* Copy buckets
-                        * there shouldn't be any null arcs here, but just to be safe 
-                        * */
-                       if (ring[i].arc->bcount > 0 && ring[j].arc->bcount > 0)
-                       {
-                               ReebArcIterator iter1, iter2;
-                               EmbedBucket *bucket1 = NULL, *bucket2 = NULL;
-                               
-                               initArcIterator(&iter1, ring[i].arc, node);
-                               initArcIterator(&iter2, ring[j].arc, node);
-                               
-                               bucket1 = nextBucket(&iter1);
-                               bucket2 = nextBucket(&iter2);
-                       
-                               /* Make sure they both start at the same value */       
-                               while(bucket1 && bucket1->val < bucket2->val)
-                               {
-                                       bucket1 = nextBucket(&iter1);
-                               }
-                               
-                               while(bucket2 && bucket2->val < bucket1->val)
-                               {
-                                       bucket2 = nextBucket(&iter2);
-                               }
-               
-               
-                               for ( ;bucket1 && bucket2; bucket1 = nextBucket(&iter1), bucket2 = nextBucket(&iter2))
-                               {
-                                       /* copy and mirror back to bucket2 */                   
-                                       bucket2->nv = bucket1->nv;
-                                       VECCOPY(bucket2->p, bucket1->p);
-                                       mirrorAlongAxis(bucket2->p, node->p, normal);
-                               }
-                       }
-               }
-       }
-
-       MEM_freeN(ring);
-}
-
-void reestablishAxialSymmetry(ReebNode *node, int depth, float axis[3])
-{
-       ReebArc *arc1 = NULL;
-       ReebArc *arc2 = NULL;
-       ReebNode *node1 = NULL, *node2 = NULL;
-       float limit = G.scene->toolsettings->skgen_symmetry_limit;
-       float nor[3], vec[3], p[3];
-       int i;
-       
-       for (i = 0; node->arcs[i] != NULL; i++)
-       {
-               ReebArc *connectedArc = node->arcs[i];
-               
-               /* depth is store as a negative in flag. symmetry level is positive */
-               if (connectedArc->flags == -depth)
-               {
-                       if (arc1 == NULL)
-                       {
-                               arc1 = connectedArc;
-                               node1 = OTHER_NODE(arc1, node);
-                       }
-                       else
-                       {
-                               arc2 = connectedArc;
-                               node2 = OTHER_NODE(arc2, node);
-                               break; /* Can stop now, the two arcs have been found */
-                       }
-               }
-       }
-       
-       /* shouldn't happen, but just to be sure */
-       if (node1 == NULL || node2 == NULL)
-       {
-               return;
-       }
-       
-       VecSubf(p, node1->p, node->p);
-       Crossf(vec, p, axis);
-       Crossf(nor, vec, axis);
-       
-       /* mirror node2 along axis */
-       VECCOPY(p, node2->p);
-       mirrorAlongAxis(p, node->p, nor);
-       
-       /* check if it's within limit before continuing */
-       if (VecLenf(node1->p, p) <= limit)
-       {
-       
-               /* average with node1 */
-               VecAddf(node1->p, node1->p, p);
-               VecMulf(node1->p, 0.5f);
-               
-               /* mirror back on node2 */
-               VECCOPY(node2->p, node1->p);
-               mirrorAlongAxis(node2->p, node->p, nor);
-               
-               /* Merge buckets
-                * there shouldn't be any null arcs here, but just to be safe 
-                * */
-               if (arc1->bcount > 0 && arc2->bcount > 0)
-               {
-                       ReebArcIterator iter1, iter2;
-                       EmbedBucket *bucket1 = NULL, *bucket2 = NULL;
-                       
-                       initArcIterator(&iter1, arc1, node);
-                       initArcIterator(&iter2, arc2, node);
-                       
-                       bucket1 = nextBucket(&iter1);
-                       bucket2 = nextBucket(&iter2);
-               
-                       /* Make sure they both start at the same value */       
-                       while(bucket1 && bucket1->val < bucket2->val)
-                       {
-                               bucket1 = nextBucket(&iter1);
-                       }
-                       
-                       while(bucket2 && bucket2->val < bucket1->val)
-                       {
-                               bucket2 = nextBucket(&iter2);
-                       }
-       
-       
-                       for ( ;bucket1 && bucket2; bucket1 = nextBucket(&iter1), bucket2 = nextBucket(&iter2))
-                       {
-                               bucket1->nv += bucket2->nv; /* add counts */
-                               
-                               /* mirror on axis */
-                               mirrorAlongAxis(bucket2->p, node->p, nor);
-                               /* add bucket2 in bucket1 */
-                               VecLerpf(bucket1->p, bucket1->p, bucket2->p, (float)bucket2->nv / (float)(bucket1->nv));
-       
-                               /* copy and mirror back to bucket2 */                   
-                               bucket2->nv = bucket1->nv;
-                               VECCOPY(bucket2->p, bucket1->p);
-                               mirrorAlongAxis(bucket2->p, node->p, nor);
-                       }
-               }
-       }
-}
-
-void markdownSecondarySymmetry(ReebNode *node, int depth, int level)
-{
-       float axis[3] = {0, 0, 0};
-       int count = 0;
-       int i;
-
-       /* Only reestablish spatial symmetry if needed */
-       if (G.scene->toolsettings->skgen_options & SKGEN_SYMMETRY)
-       {
-               /* count the number of branches in this symmetry group
-                * and determinte the axis of symmetry
-                *  */  
-               for (i = 0; node->arcs[i] != NULL; i++)
-               {
-                       ReebArc *connectedArc = node->arcs[i];
-                       
-                       /* depth is store as a negative in flag. symmetry level is positive */
-                       if (connectedArc->flags == -depth)
-                       {
-                               count++;
-                       }
-                       /* If arc is on the axis */
-                       else if (connectedArc->flags == level)
-                       {
-                               VecAddf(axis, axis, connectedArc->v1->p);
-                               VecSubf(axis, axis, connectedArc->v2->p);
-                       }
-               }
-       
-               Normalize(axis);
-       
-               /* Split between axial and radial symmetry */
-               if (count == 2)
-               {
-                       reestablishAxialSymmetry(node, depth, axis);
-               }
-               else
-               {
-                       reestablishRadialSymmetry(node, depth, axis);
-               }
-       }
-
-       /* markdown secondary symetries */      
-       for (i = 0; node->arcs[i] != NULL; i++)
-       {
-               ReebArc *connectedArc = node->arcs[i];
-               
-               if (connectedArc->flags == -depth)
-               {
-                       /* markdown symmetry for branches corresponding to the depth */
-                       markdownSymmetryArc(connectedArc, node, level + 1);
-               }
-       }
-}
-
-void markdownSymmetryArc(ReebArc *arc, ReebNode *node, int level)
-{
-       int i;
-       arc->flags = level;
-       
-       node = OTHER_NODE(arc, node);
-       
-       for (i = 0; node->arcs[i] != NULL; i++)
-       {
-               ReebArc *connectedArc = node->arcs[i];
-               
-               if (connectedArc != arc)
-               {
-                       ReebNode *connectedNode = OTHER_NODE(connectedArc, node);
-                       
-                       /* symmetry level is positive value, negative values is subtree depth */
-                       connectedArc->flags = -subtreeDepth(connectedNode, connectedArc);
-               }
-       }
-
-       arc = NULL;
-
-       for (i = 0; node->arcs[i] != NULL; i++)
-       {
-               int issymmetryAxis = 0;
-               ReebArc *connectedArc = node->arcs[i];
-               
-               /* only arcs not already marked as symetric */
-               if (connectedArc->flags < 0)
-               {
-                       int j;
-                       
-                       /* true by default */
-                       issymmetryAxis = 1;
-                       
-                       for (j = 0; node->arcs[j] != NULL && issymmetryAxis == 1; j++)
-                       {
-                               ReebArc *otherArc = node->arcs[j];
-                               
-                               /* different arc, same depth */
-                               if (otherArc != connectedArc && otherArc->flags == connectedArc->flags)
-                               {
-                                       /* not on the symmetry axis */
-                                       issymmetryAxis = 0;
-                               } 
-                       }
-               }
-               
-               /* arc could be on the symmetry axis */
-               if (issymmetryAxis == 1)
-               {
-                       /* no arc as been marked previously, keep this one */
-                       if (arc == NULL)
-                       {
-                               arc = connectedArc;
-                       }
-                       else
-                       {
-                               /* there can't be more than one symmetry arc */
-                               arc = NULL;
-                               break;
-                       }
-               }
-       }
-       
-       /* go down the arc continuing the symmetry axis */
-       if (arc)
-       {
-               markdownSymmetryArc(arc, node, level);
-       }
-
-       
-       /* secondary symmetry */
-       for (i = 0; node->arcs[i] != NULL; i++)
-       {
-               ReebArc *connectedArc = node->arcs[i];
-               
-               /* only arcs not already marked as symetric and is not the next arc on the symmetry axis */
-               if (connectedArc->flags < 0)
-               {
-                       /* subtree depth is store as a negative value in the flag */
-                       markdownSecondarySymmetry(node, -connectedArc->flags, level);
-               }
-       }
-}
-
-void markdownSymmetry(ReebGraph *rg)
-{
-       ReebNode *node;
-       ReebArc *arc;
-       /* only for Acyclic graphs */
-       int cyclic = isGraphCyclic(rg);
-       
-       /* mark down all arcs as non-symetric */
-       for (arc = rg->arcs.first; arc; arc = arc->next)
-       {
-               arc->flags = 0;
-       }
-       
-       /* mark down all nodes as not on the symmetry axis */
-       for (node = rg->nodes.first; node; node = node->next)
-       {
-               node->flags = 0;
-       }
-
-       /* node list is sorted, so lowest node is always the head (by design) */
-       node = rg->nodes.first;
-       
-       /* only work on acyclic graphs and if only one arc is incident on the first node */
-       if (cyclic == 0 && countConnectedArcs(rg, node) == 1)
-       {
-               arc = node->arcs[0];
-               
-               markdownSymmetryArc(arc, node, 1);
-
-               /* mark down non-symetric arcs */
-               for (arc = rg->arcs.first; arc; arc = arc->next)
-               {
-                       if (arc->flags < 0)
-                       {
-                               arc->flags = 0;
-                       }
-                       else
-                       {
-                               /* mark down nodes with the lowest level symmetry axis */
-                               if (arc->v1->flags == 0 || arc->v1->flags > arc->flags)
-                               {
-                                       arc->v1->flags = arc->flags;
-                               }
-                               if (arc->v2->flags == 0 || arc->v2->flags > arc->flags)
-                               {
-                                       arc->v2->flags = arc->flags;
-                               }
-                       }
-               }
-       }
-}
 
 /**************************************** SUBDIVISION ALGOS ******************************************/
 
@@ -4789,7 +4254,7 @@ EditBone * subdivideByAngle(ReebArc *arc, ReebNode *head, ReebNode *tail)
        return lastBone;
 }
 
-float calcCorrelation(ReebArc *arc, int start, int end, float v0[3], float n[3])
+float calcVariance(ReebArc *arc, int start, int end, float v0[3], float n[3])
 {
        int len = 2 + abs(end - start);
        
@@ -4837,19 +4302,47 @@ float calcCorrelation(ReebArc *arc, int start, int end, float v0[3], float n[3])
                /* adding start(0) and end(1) values to s_t */
                s_t += (avg_t * avg_t) + (1 - avg_t) * (1 - avg_t);
                
-               return 1.0f - s_xyz / s_t; 
+               return s_xyz / s_t; 
        }
        else
        {
-               return 1.0f;
+               return 0;
+       }
+}
+
+float calcDistance(ReebArc *arc, int start, int end, float head[3], float tail[3])
+{
+       ReebArcIterator iter;
+       EmbedBucket *bucket = NULL;
+       float max_dist = 0;
+       
+       /* calculate maximum distance */
+       for (initArcIterator2(&iter, arc, start, end), bucket = nextBucket(&iter);
+               bucket;
+               bucket = nextBucket(&iter))
+       {
+               float v1[3], v2[3], c[3];
+               float dist;
+               
+               VecSubf(v1, head, tail);
+               VecSubf(v2, bucket->p, tail);
+
+               Crossf(c, v1, v2);
+               
+               dist = Inpf(c, c) / Inpf(v1, v1);
+               
+               max_dist = dist > max_dist ? dist : max_dist;
        }
+       
+       
+       return max_dist; 
 }
 
 EditBone * subdivideByCorrelation(ReebArc *arc, ReebNode *head, ReebNode *tail)
 {
        ReebArcIterator iter;
        float n[3];
-       float CORRELATION_THRESHOLD = G.scene->toolsettings->skgen_correlation_limit;
+       float ADAPTIVE_THRESHOLD = G.scene->toolsettings->skgen_correlation_limit;
        EditBone *lastBone = NULL;
        
        /* init iterator to get start and end from head */
@@ -4858,15 +4351,17 @@ EditBone * subdivideByCorrelation(ReebArc *arc, ReebNode *head, ReebNode *tail)
        /* Calculate overall */
        VecSubf(n, arc->buckets[iter.end].p, head->p);
        
-       if (G.scene->toolsettings->skgen_options & SKGEN_CUT_CORRELATION && 
-               calcCorrelation(arc, iter.start, iter.end, head->p, n) < CORRELATION_THRESHOLD)
+       if (G.scene->toolsettings->skgen_options & SKGEN_CUT_CORRELATION)
        {
                EmbedBucket *bucket = NULL;
                EmbedBucket *previous = NULL;
                EditBone *child = NULL;
                EditBone *parent = NULL;
+               float normal[3] = {0, 0, 0};
+               float avg_normal[3];
+               int total = 0;
                int boneStart = iter.start;
-
+               
                parent = add_editbone("Bone");
                parent->flag = BONE_SELECTED|BONE_TIPSEL|BONE_ROOTSEL;
                VECCOPY(parent->head, head->p);
@@ -4875,12 +4370,46 @@ EditBone * subdivideByCorrelation(ReebArc *arc, ReebNode *head, ReebNode *tail)
                        bucket;
                        previous = bucket, bucket = nextBucket(&iter))
                {
-                       /* Calculate normal */
-                       VecSubf(n, bucket->p, parent->head);
+                       float btail[3];
+                       float value = 0;
 
-                       if (calcCorrelation(arc, boneStart, iter.index, parent->head, n) < CORRELATION_THRESHOLD)
+                       if (G.scene->toolsettings->skgen_options & SKGEN_STICK_TO_EMBEDDING)
                        {
-                               VECCOPY(parent->tail, previous->p);
+                               VECCOPY(btail, bucket->p);
+                       }
+                       else
+                       {
+                               float length;
+                               
+                               /* Calculate normal */
+                               VecSubf(n, bucket->p, parent->head);
+                               length = Normalize(n);
+                               
+                               total += 1;
+                               VecAddf(normal, normal, n);
+                               VECCOPY(avg_normal, normal);
+                               VecMulf(avg_normal, 1.0f / total);
+                                
+                               VECCOPY(btail, avg_normal);
+                               VecMulf(btail, length);
+                               VecAddf(btail, btail, parent->head);
+                       }
+
+                       if (G.scene->toolsettings->skgen_options & SKGEN_ADAPTIVE_DISTANCE)
+                       {
+                               value = calcDistance(arc, boneStart, iter.index, parent->head, btail);
+                       }
+                       else
+                       {
+                               float n[3];
+                               
+                               VecSubf(n, btail, parent->head);
+                               value = calcVariance(arc, boneStart, iter.index, parent->head, n);
+                       }
+
+                       if (value > ADAPTIVE_THRESHOLD)
+                       {
+                               VECCOPY(parent->tail, btail);
 
                                child = add_editbone("Bone");
                                VECCOPY(child->head, parent->tail);
@@ -4889,6 +4418,9 @@ EditBone * subdivideByCorrelation(ReebArc *arc, ReebNode *head, ReebNode *tail)
                                
                                parent = child; // new child is next parent
                                boneStart = iter.index; // start from end
+                               
+                               normal[0] = normal[1] = normal[2] = 0;
+                               total = 0;
                        }
                }
 
@@ -4906,7 +4438,7 @@ float arcLengthRatio(ReebArc *arc)
        float embedLength = 0.0f;
        int i;
        
-       arcLength = VecLenf(arc->v1->p, arc->v2->p);
+       arcLength = VecLenf(arc->head->p, arc->tail->p);
        
        if (arc->bcount > 0)
        {
@@ -4916,8 +4448,8 @@ float arcLengthRatio(ReebArc *arc)
                        embedLength += VecLenf(arc->buckets[i - 1].p, arc->buckets[i].p);
                }
                /* Add head and tail -> embedding vectors */
-               embedLength += VecLenf(arc->v1->p, arc->buckets[0].p);
-               embedLength += VecLenf(arc->v2->p, arc->buckets[arc->bcount - 1].p);
+               embedLength += VecLenf(arc->head->p, arc->buckets[0].p);
+               embedLength += VecLenf(arc->tail->p, arc->buckets[arc->bcount - 1].p);
        }
        else
        {
@@ -5049,8 +4581,6 @@ void generateSkeletonFromReebGraph(ReebGraph *rg)
        {
                exit_editmode(EM_FREEDATA|EM_FREEUNDO|EM_WAITCURSOR); // freedata, and undo
        }
-
-       setcursor_space(SPACE_VIEW3D, CURSOR_WAIT);
        
        dst = add_object(OB_ARMATURE);
        base_init_from_view3d(BASACT, G.vd);
@@ -5067,7 +4597,7 @@ void generateSkeletonFromReebGraph(ReebGraph *rg)
 
        arcBoneMap = BLI_ghash_new(BLI_ghashutil_ptrhash, BLI_ghashutil_ptrcmp);
        
-       markdownSymmetry(rg);
+       BLI_markdownSymmetry((BGraph*)rg, rg->nodes.first, G.scene->toolsettings->skgen_symmetry_limit);
        
        for (arc = rg->arcs.first; arc; arc = arc->next) 
        {
@@ -5077,43 +4607,43 @@ void generateSkeletonFromReebGraph(ReebGraph *rg)
 
                /* Find out the direction of the arc through simple heuristics (in order of priority) :
                 * 
-                * 1- Arcs on primary symmetry axis (flags == 1) point up (head: high weight -> tail: low weight)
+                * 1- Arcs on primary symmetry axis (symmetry == 1) point up (head: high weight -> tail: low weight)
                 * 2- Arcs starting on a primary axis point away from it (head: node on primary axis)
                 * 3- Arcs point down (head: low weight -> tail: high weight)
                 *
-                * Finally, the arc direction is stored in its flags: 1 (low -> high), -1 (high -> low)
+                * Finally, the arc direction is stored in its flag: 1 (low -> high), -1 (high -> low)
                 */
 
                /* if arc is a symmetry axis, internal bones go up the tree */          
-               if (arc->flags == 1 && arc->v2->degree != 1)
+               if (arc->symmetry_level == 1 && arc->tail->degree != 1)
                {
-                       head = arc->v2;
-                       tail = arc->v1;
+                       head = arc->tail;
+                       tail = arc->head;
                        
-                       arc->flags = -1; /* mark arc direction */
+                       arc->flag = -1; /* mark arc direction */
                }
                /* Bones point AWAY from the symmetry axis */
-               else if (arc->v1->flags == 1)
+               else if (arc->head->symmetry_level == 1)
                {
-                       head = arc->v1;
-                       tail = arc->v2;
+                       head = arc->head;
+                       tail = arc->tail;
                        
-                       arc->flags = 1; /* mark arc direction */
+                       arc->flag = 1; /* mark arc direction */
                }
-               else if (arc->v2->flags == 1)
+               else if (arc->tail->symmetry_level == 1)
                {
-                       head = arc->v2;
-                       tail = arc->v1;
+                       head = arc->tail;
+                       tail = arc->head;
                        
-                       arc->flags = -1; /* mark arc direction */
+                       arc->flag = -1; /* mark arc direction */
                }
                /* otherwise, always go from low weight to high weight */
                else
                {
-                       head = arc->v1;
-                       tail = arc->v2;
+                       head = arc->head;
+                       tail = arc->tail;
                        
-                       arc->flags = 1; /* mark arc direction */
+                       arc->flag = 1; /* mark arc direction */
                }
                
                /* Loop over subdivision methods */     
@@ -5155,12 +4685,12 @@ void generateSkeletonFromReebGraph(ReebGraph *rg)
                ReebArc *incomingArc = NULL;
                int i;
 
-               for (i = 0; node->arcs[i] != NULL; i++)
+               for (i = 0; i < node->degree; i++)
                {
-                       arc = node->arcs[i];
+                       arc = (ReebArc*)node->arcs[i];
 
                        /* if arc is incoming into the node */
-                       if ((arc->v1 == node && arc->flags == -1) || (arc->v2 == node && arc->flags == 1))
+                       if ((arc->head == node && arc->flag == -1) || (arc->tail == node && arc->flag == 1))
                        {
                                if (incomingArc == NULL)
                                {
@@ -5181,12 +4711,12 @@ void generateSkeletonFromReebGraph(ReebGraph *rg)
                        EditBone *parentBone = BLI_ghash_lookup(arcBoneMap, incomingArc);
 
                        /* Look for outgoing arcs and parent their bones */
-                       for (i = 0; node->arcs[i] != NULL; i++)
+                       for (i = 0; i < node->degree; i++)
                        {
                                arc = node->arcs[i];
 
                                /* if arc is outgoing from the node */
-                               if ((arc->v1 == node && arc->flags == 1) || (arc->v2 == node && arc->flags == -1))
+                               if ((arc->head == node && arc->flag == 1) || (arc->tail == node && arc->flag == -1))
                                {
                                        EditBone *childBone = BLI_ghash_lookup(arcBoneMap, arc);
 
@@ -5204,89 +4734,21 @@ void generateSkeletonFromReebGraph(ReebGraph *rg)
        }
        
        BLI_ghash_free(arcBoneMap, NULL, NULL);
-
-       setcursor_space(SPACE_VIEW3D, CURSOR_EDIT);
        
        BIF_undo_push("Generate Skeleton");
 }
 
 void generateSkeleton(void)
 {
-       EditMesh *em = G.editMesh;
-       ReebGraph *rg = NULL;
-       int i;
+       ReebGraph *reebg;
        
-       if (em == NULL)
-               return;
-
        setcursor_space(SPACE_VIEW3D, CURSOR_WAIT);
-
-       if (weightFromDistance(em) == 0)
-       {
-               error("No selected vertex\n");
-               return;
-       }
-               
-       renormalizeWeight(em, 1.0f);
-       
-       weightToHarmonic(em);
-       
-#ifdef DEBUG_REEB
-       weightToVCol(em);
-#endif
-       
-       rg = generateReebGraph(em, G.scene->toolsettings->skgen_resolution);
-
-       verifyBuckets(rg);
-       
-       /* Remove arcs without embedding */
-       filterNullReebGraph(rg);
-
-       verifyBuckets(rg);
-
-
-       i = 1;
-       /* filter until there's nothing more to do */
-       while (i == 1)
-       {
-               i = 0; /* no work done yet */
-               
-               if (G.scene->toolsettings->skgen_options & SKGEN_FILTER_EXTERNAL)
-               {
-                       i |= filterExternalReebGraph(rg, G.scene->toolsettings->skgen_threshold_external * G.scene->toolsettings->skgen_resolution);
-               }
-       
-               verifyBuckets(rg);
        
-               if (G.scene->toolsettings->skgen_options & SKGEN_FILTER_INTERNAL)
-               {
-                       i |= filterInternalReebGraph(rg, G.scene->toolsettings->skgen_threshold_internal * G.scene->toolsettings->skgen_resolution);
-               }
-       }
+       reebg = BIF_ReebGraphFromEditMesh();
 
-       verifyBuckets(rg);
+       generateSkeletonFromReebGraph(reebg);
 
-       repositionNodes(rg);
-       
-       verifyBuckets(rg);
+       REEB_freeGraph(reebg);
 
-       /* Filtering might have created degree 2 nodes, so remove them */
-       removeNormalNodes(rg);
-       
-       verifyBuckets(rg);
-
-       for(i = 0; i <  G.scene->toolsettings->skgen_postpro_passes; i++)
-       {
-               postprocessGraph(rg, G.scene->toolsettings->skgen_postpro);
-       }
-
-       buildAdjacencyList(rg);
-       
-       sortNodes(rg);
-       
-       sortArcs(rg);
-       
-       generateSkeletonFromReebGraph(rg);
-
-       freeGraph(rg);
+       setcursor_space(SPACE_VIEW3D, CURSOR_EDIT);
 }
index 85fb5815c3e0eeaf2005ab779fc6371657d4b3dd..a174c5c199f99a527c859efc5f723dd8e87d6dec 100644 (file)
@@ -34,6 +34,7 @@
 #include "DNA_scene_types.h"
 #include "DNA_space_types.h"
 #include "DNA_meshdata_types.h"
+#include "DNA_armature_types.h"
 
 #include "MEM_guardedalloc.h"
 
 #include "BLI_arithb.h"
 #include "BLI_editVert.h"
 #include "BLI_edgehash.h"
+#include "BLI_ghash.h"
 
 #include "BDR_editobject.h"
 
+#include "BMF_Api.h"
+
 #include "BIF_editmesh.h"
 #include "BIF_editarmature.h"
 #include "BIF_interface.h"
 #include "BIF_toolbox.h"
 #include "BIF_graphics.h"
+#include "BIF_gl.h"
+#include "BIF_resources.h"
 
 #include "BKE_global.h"
 #include "BKE_utildefines.h"
 
 #include "reeb.h"
 
+/* REPLACE WITH NEW ONE IN UTILDEFINES ONCE PATCH IS APPLIED */
+#define FTOCHAR(val) (val<=0.0f)? 0 : ((val>(1.0f-0.5f/255.0f))? 255 : (char)((255.0f*val)+0.5f))
+
+ReebGraph *GLOBAL_RG = NULL;
+ReebGraph *FILTERED_RG = NULL;
+
 /*
  * Skeleton generation algorithm based on: 
  * "Harmonic Skeleton for Realistic Character Animation"
  * SIGGRAPH 2007
  * 
  * */
+#define DEBUG_REEB
+
+typedef enum {
+       MERGE_LOWER,
+       MERGE_HIGHER,
+       MERGE_APPEND
+} MergeDirection;
 
 int mergeArcs(ReebGraph *rg, ReebArc *a0, ReebArc *a1);
+void mergeArcEdges(ReebGraph *rg, ReebArc *aDst, ReebArc *aSrc, MergeDirection direction);
 int mergeConnectedArcs(ReebGraph *rg, ReebArc *a0, ReebArc *a1);
 EditEdge * NextEdgeForVert(EditMesh *em, EditVert *v);
+void mergeArcFaces(ReebGraph *rg, ReebArc *aDst, ReebArc *aSrc);
+void addFacetoArc(ReebArc *arc, EditFace *efa);
 
-/***************************************** BUCKET UTILS **********************************************/
+void REEB_RadialSymmetry(BNode* root_node, RadialArc* ring, int count);
+void REEB_AxialSymmetry(BNode* root_node, BNode* node1, BNode* node2, struct BArc* barc1, BArc* barc2);
 
-void addVertToBucket(EmbedBucket *b, float co[3])
-{
-       b->nv++;
-       VecLerpf(b->p, b->p, co, 1.0f / b->nv);
-}
 
-void removeVertFromBucket(EmbedBucket *b, float co[3])
+/***************************************** UTILS **********************************************/
+
+void REEB_freeArc(BArc *barc)
 {
-       VecMulf(b->p, (float)b->nv);
-       VecSubf(b->p, b->p, co);
-       b->nv--;
-       VecMulf(b->p, 1.0f / (float)b->nv);
+       ReebArc *arc = (ReebArc*)barc;
+       BLI_freelistN(&arc->edges);
+       
+       if (arc->buckets)
+               MEM_freeN(arc->buckets);
+               
+       if (arc->faces)
+               BLI_ghash_free(arc->faces, NULL, NULL);
+       
+       MEM_freeN(arc);
 }
 
-void mergeBuckets(EmbedBucket *bDst, EmbedBucket *bSrc)
+void REEB_freeGraph(ReebGraph *rg)
 {
-       if (bDst->nv > 0 && bSrc->nv > 0)
+       ReebArc *arc;
+       ReebNode *node;
+       
+       // free nodes
+       for( node = rg->nodes.first; node; node = node->next )
        {
-               bDst->nv += bSrc->nv;
-               VecLerpf(bDst->p, bDst->p, bSrc->p, (float)bSrc->nv / (float)(bDst->nv));
+               BLI_freeNode((BGraph*)rg, (BNode*)node);
        }
-       else if (bSrc->nv > 0)
+       BLI_freelistN(&rg->nodes);
+       
+       // free arcs
+       arc = rg->arcs.first;
+       while( arc )
        {
-               bDst->nv = bSrc->nv;
-               VECCOPY(bDst->p, bSrc->p);
+               ReebArc *next = arc->next;
+               REEB_freeArc((BArc*)arc);
+               arc = next;
+       }
+       
+       // free edge map
+       BLI_edgehash_free(rg->emap, NULL);
+       
+       /* free linked graph */
+       if (rg->link_up)
+       {
+               REEB_freeGraph(rg->link_up);
        }
+       
+       MEM_freeN(rg);
 }
 
-void mergeArcBuckets(ReebArc *aDst, ReebArc *aSrc, float start, float end)
+ReebGraph * newReebGraph()
 {
-       if (aDst->bcount > 0 && aSrc->bcount > 0)
-       {
-               int indexDst = 0, indexSrc = 0;
-               
-               start = MAX3(start, aDst->buckets[0].val, aSrc->buckets[0].val);
-               
-               while(indexDst < aDst->bcount && aDst->buckets[indexDst].val < start)
-               {
-                       indexDst++;
-               }
+       ReebGraph *rg;
+       rg = MEM_callocN(sizeof(ReebGraph), "reeb graph");
+       
+       rg->totnodes = 0;
+       rg->emap = BLI_edgehash_new();
+       
+       
+       rg->free_arc = REEB_freeArc;
+       rg->free_node = NULL;
+       rg->radial_symmetry = REEB_RadialSymmetry;
+       rg->axial_symmetry = REEB_AxialSymmetry;
+       
+       return rg;
+}
 
-               while(indexSrc < aSrc->bcount && aSrc->buckets[indexSrc].val < start)
-               {
-                       indexSrc++;
-               }
-               
-               for( ;  indexDst < aDst->bcount &&
-                               indexSrc < aSrc->bcount &&
-                               aDst->buckets[indexDst].val <= end &&
-                               aSrc->buckets[indexSrc].val <= end
-                               
-                        ;      indexDst++, indexSrc++)
-               {
-                       mergeBuckets(aDst->buckets + indexDst, aSrc->buckets + indexSrc);
-               }
+void BIF_flagMultiArcs(ReebGraph *rg, int flag)
+{
+       for ( ; rg; rg = rg->link_up)
+       {
+               BLI_flagArcs((BGraph*)rg, flag);
        }
 }
 
-void allocArcBuckets(ReebArc *arc)
+ReebNode * addNode(ReebGraph *rg, EditVert *eve, float weight)
 {
-       int i;
-       float start = ceil(arc->v1->weight);
-       arc->bcount = (int)(floor(arc->v2->weight) - start) + 1;
+       ReebNode *node = NULL;
        
-       if (arc->bcount > 0)
+       node = MEM_callocN(sizeof(ReebNode), "reeb node");
+       
+       node->flag = 0; // clear flag on init
+       node->symmetry_level = 0;
+       node->arcs = NULL;
+       node->degree = 0;
+       node->weight = weight;
+       node->index = rg->totnodes;
+       VECCOPY(node->p, eve->co);      
+       
+       BLI_addtail(&rg->nodes, node);
+       rg->totnodes++;
+       
+       return node;
+}
+
+ReebNode * copyNode(ReebGraph *rg, ReebNode *node)
+{
+       ReebNode *cp_node = NULL;
+       
+       cp_node = MEM_callocN(sizeof(ReebNode), "reeb node copy");
+       
+       memcpy(cp_node, node, sizeof(ReebNode));
+       
+       cp_node->prev = NULL;
+       cp_node->next = NULL;
+       cp_node->arcs = NULL;
+       
+       cp_node->link_up = NULL;
+       cp_node->link_down = NULL;
+       
+       BLI_addtail(&rg->nodes, cp_node);
+       rg->totnodes++;
+       
+       return cp_node; 
+}
+
+void relinkNodes(ReebGraph *low_rg, ReebGraph *high_rg)
+{
+       ReebNode *low_node, *high_node;
+       
+       if (low_rg == NULL || high_rg == NULL)
        {
-               arc->buckets = MEM_callocN(sizeof(EmbedBucket) * arc->bcount, "embed bucket");
-               
-               for(i = 0; i < arc->bcount; i++)
+               return;
+       }
+       
+       for (low_node = low_rg->nodes.first; low_node; low_node = low_node->next)
+       {
+               for (high_node = high_rg->nodes.first; high_node; high_node = high_node->next)
                {
-                       arc->buckets[i].val = start + i;
+                       if (low_node->index == high_node->index)
+                       {
+                               high_node->link_down = low_node;
+                               low_node->link_up = high_node;
+                               break;
+                       }
                }
        }
-       else
+}
+
+ReebNode *BIF_otherNodeFromIndex(ReebArc *arc, ReebNode *node)
+{
+       return (arc->head->index == node->index) ? arc->tail : arc->head;
+}
+
+ReebNode *BIF_lowestLevelNode(ReebNode *node)
+{
+       while (node->link_down)
        {
-               arc->buckets = NULL;
+               node = node->link_down;
        }
        
+       return node;
 }
 
-void resizeArcBuckets(ReebArc *arc)
+ReebArc * copyArc(ReebGraph *rg, ReebArc *arc)
 {
-       EmbedBucket *oldBuckets = arc->buckets;
-       int oldBCount = arc->bcount;
+       ReebArc *cp_arc;
+       ReebNode *node;
        
-       allocArcBuckets(arc);
+       cp_arc = MEM_callocN(sizeof(ReebArc), "reeb arc copy");
+
+       memcpy(cp_arc, arc, sizeof(ReebArc));
        
-       if (oldBCount != 0 && arc->bcount != 0)
+       cp_arc->link_up = arc;
+       
+       cp_arc->head = NULL;
+       cp_arc->tail = NULL;
+
+       cp_arc->prev = NULL;
+       cp_arc->next = NULL;
+
+       cp_arc->edges.first = NULL;
+       cp_arc->edges.last = NULL;
+
+       /* copy buckets */      
+       cp_arc->buckets = MEM_callocN(sizeof(EmbedBucket) * cp_arc->bcount, "embed bucket");
+       memcpy(cp_arc->buckets, arc->buckets, sizeof(EmbedBucket) * cp_arc->bcount);
+       
+       /* copy faces map */
+       cp_arc->faces = BLI_ghash_new(BLI_ghashutil_ptrhash, BLI_ghashutil_ptrcmp);
+       mergeArcFaces(rg, cp_arc, arc);
+       
+       /* find corresponding head and tail */
+       for (node = rg->nodes.first; node && (cp_arc->head == NULL || cp_arc->tail == NULL); node = node->next)
        {
-               int oldStart = (int)oldBuckets[0].val;
-               int oldEnd = (int)oldBuckets[oldBCount - 1].val;
-               int newStart = (int)arc->buckets[0].val;
-               int newEnd = (int)arc->buckets[arc->bcount - 1].val;
-               int oldOffset = 0;
-               int newOffset = 0;
-               int len;
-               
-               if (oldStart < newStart)
+               if (node->index == arc->head->index)
                {
-                       oldOffset = newStart - oldStart;
+                       cp_arc->head = node;
                }
-               else
+               else if (node->index == arc->tail->index)
                {
-                       newOffset = oldStart - newStart;
+                       cp_arc->tail = node;
                }
-               
-               len = MIN2(oldEnd - (oldStart + oldOffset) + 1, newEnd - (newStart - newOffset) + 1);
-               
-               memcpy(arc->buckets + newOffset, oldBuckets + oldOffset, len * sizeof(EmbedBucket)); 
        }
+       
+       BLI_addtail(&rg->arcs, cp_arc);
+       
+       return cp_arc;
+}
 
-       if (oldBuckets != NULL)
+ReebGraph * copyReebGraph(ReebGraph *rg)
+{
+       ReebNode *node;
+       ReebArc *arc;
+       ReebGraph *cp_rg = newReebGraph();
+       
+       cp_rg->resolution = rg->resolution;
+       cp_rg->length = rg->length;
+       cp_rg->link_up = rg;
+
+       /* Copy nodes */        
+       for (node = rg->nodes.first; node; node = node->next)
        {
-               MEM_freeN(oldBuckets);
+               copyNode(cp_rg, node);
+       }
+       
+       /* Copy arcs */
+       for (arc = rg->arcs.first; arc; arc = arc->next)
+       {
+               copyArc(cp_rg, arc);
        }
+       
+       BLI_rebuildAdjacencyList((BGraph*)cp_rg);
+       
+       return cp_rg;
 }
-/***************************************** UTILS **********************************************/
 
 ReebEdge * copyEdge(ReebEdge *edge)
 {
@@ -213,7 +346,8 @@ ReebEdge * copyEdge(ReebEdge *edge)
 void printArc(ReebArc *arc)
 {
        ReebEdge *edge;
-       printf("arc: (%i)%f -> (%i)%f\n", arc->v1->index, arc->v1->weight, arc->v2->index, arc->v2->weight);
+       ReebNode *head = (ReebNode*)arc->head;
+       printf("arc: (%i)%f -> (%i)%f\n", head->index, head->weight, head->index, head->weight);
        
        for(edge = arc->edges.first; edge ; edge = edge->next)
        {
@@ -221,51 +355,36 @@ void printArc(ReebArc *arc)
        }
 }
 
-void freeArc(ReebArc *arc)
+#ifdef DEBUG_REEB_NODE
+void NodeDegreeDecrement(ReebGraph *rg, ReebNode *node)
 {
-       BLI_freelistN(&arc->edges);
-       
-       if (arc->buckets)
-               MEM_freeN(arc->buckets);
-       
-       MEM_freeN(arc);
-}
+       node->degree--;
 
-void freeGraph(ReebGraph *rg)
-{
-       ReebArc *arc;
-       ReebNode *node;
-       
-       // free nodes
-       for( node = rg->nodes.first; node; node = node->next )
+       if (node->degree == 0)
        {
-               // Free adjacency lists
-               if (node->arcs != NULL)
-               {
-                       MEM_freeN(node->arcs);
-               }
+               printf("would remove node %i\n", node->index);
        }
-       BLI_freelistN(&rg->nodes);
-       
-       // free arcs
-       arc = rg->arcs.first;
-       while( arc )
+}
+
+void NodeDegreeIncrement(ReebGraph *rg, ReebNode *node)
+{
+       if (node->degree == 0)
        {
-               ReebArc *next = arc->next;
-               freeArc(arc);
-               arc = next;
+               printf("first connect node %i\n", node->index);
        }
-       
-       // free edge map
-       BLI_edgehash_free(rg->emap, NULL);
-       
-       MEM_freeN(rg);
+
+       node->degree++;
 }
 
+#else
+#define NodeDegreeDecrement(rg, node) {node->degree--;}
+#define NodeDegreeIncrement(rg, node) {node->degree++;}
+#endif
+
 void repositionNodes(ReebGraph *rg)
 {
-       ReebArc *arc = NULL;
-       ReebNode *node = NULL;
+       BArc *arc = NULL;
+       BNode *node = NULL;
        
        // Reset node positions
        for(node = rg->nodes.first; node; node = node->next)
@@ -275,23 +394,24 @@ void repositionNodes(ReebGraph *rg)
        
        for(arc = rg->arcs.first; arc; arc = arc->next)
        {
-               if (arc->bcount > 0)
+               if (((ReebArc*)arc)->bcount > 0)
                {
                        float p[3];
                        
-                       VECCOPY(p, arc->buckets[0].p);
-                       VecMulf(p, 1.0f / arc->v1->degree);
-                       VecAddf(arc->v1->p, arc->v1->p, p);
+                       VECCOPY(p, ((ReebArc*)arc)->buckets[0].p);
+                       VecMulf(p, 1.0f / arc->head->degree);
+                       VecAddf(arc->head->p, arc->head->p, p);
                        
-                       VECCOPY(p, arc->buckets[arc->bcount - 1].p);
-                       VecMulf(p, 1.0f / arc->v2->degree);
-                       VecAddf(arc->v2->p, arc->v2->p, p);
+                       VECCOPY(p, ((ReebArc*)arc)->buckets[((ReebArc*)arc)->bcount - 1].p);
+                       VecMulf(p, 1.0f / arc->tail->degree);
+                       VecAddf(arc->tail->p, arc->tail->p, p);
                }
        }
 }
 
 void verifyNodeDegree(ReebGraph *rg)
 {
+#ifdef DEBUG_REEB
        ReebNode *node = NULL;
        ReebArc *arc = NULL;
 
@@ -300,7 +420,7 @@ void verifyNodeDegree(ReebGraph *rg)
                int count = 0;
                for(arc = rg->arcs.first; arc; arc = arc->next)
                {
-                       if (arc->v1 == node || arc->v2 == node)
+                       if (arc->head == node || arc->tail == node)
                        {
                                count++;
                        }
@@ -309,7 +429,12 @@ void verifyNodeDegree(ReebGraph *rg)
                {
                        printf("degree error in node %i: expected %i got %i\n", node->index, count, node->degree);
                }
+               if (node->degree == 0)
+               {
+                       printf("zero degree node %i with weight %f\n", node->index, node->weight);
+               }
        }
+#endif
 }
 
 void verifyBuckets(ReebGraph *rg)
@@ -318,6 +443,9 @@ void verifyBuckets(ReebGraph *rg)
        ReebArc *arc = NULL;
        for(arc = rg->arcs.first; arc; arc = arc->next)
        {
+               ReebNode *head = (ReebNode*)arc->head;
+               ReebNode *tail = (ReebNode*)arc->tail;
+
                if (arc->bcount > 0)
                {
                        int i;
@@ -330,99 +458,424 @@ void verifyBuckets(ReebGraph *rg)
                                }
                        }
                        
-                       if (ceil(arc->v1->weight) < arc->buckets[0].val)
+                       if (ceil(head->weight) < arc->buckets[0].val)
                        {
                                printArc(arc);
-                               printf("alloc error in first bucket: %f should be %f \n", arc->buckets[0].val, ceil(arc->v1->weight));
+                               printf("alloc error in first bucket: %f should be %f \n", arc->buckets[0].val, ceil(head->weight));
                        }
-                       if (floor(arc->v2->weight) < arc->buckets[arc->bcount - 1].val)
+                       if (floor(tail->weight) < arc->buckets[arc->bcount - 1].val)
                        {
                                printArc(arc);
-                               printf("alloc error in last bucket: %f should be %f \n", arc->buckets[arc->bcount - 1].val, floor(arc->v2->weight));
+                               printf("alloc error in last bucket: %f should be %f \n", arc->buckets[arc->bcount - 1].val, floor(tail->weight));
                        }
                }
        }
 #endif
 }
 
-/************************************** ADJACENCY LIST *************************************************/
-
-void addArcToNodeAdjacencyList(ReebNode *node, ReebArc *arc)
+void verifyFaces(ReebGraph *rg)
 {
-       ReebArc **arclist;
-
-       for(arclist = node->arcs; *arclist; arclist++)
-       {       }
+#ifdef DEBUG_REEB
+       int total = 0;
+       ReebArc *arc = NULL;
+       for(arc = rg->arcs.first; arc; arc = arc->next)
+       {
+               total += BLI_ghash_size(arc->faces);
+       }
        
-       *arclist = arc;
+#endif
 }
 
-void buildAdjacencyList(ReebGraph *rg)
+void verifyMultiResolutionLinks(ReebGraph *rg)
 {
-       ReebNode *node = NULL;
-       ReebArc *arc = NULL;
-
-       for(node = rg->nodes.first; node; node = node->next)
+#ifdef DEBUG_REEB
+       ReebGraph *lower_rg = rg->link_up;
+       
+       if (lower_rg)
        {
-               if (node->arcs != NULL)
+               ReebArc *arc;
+               
+               for (arc = rg->arcs.first; arc; arc = arc->next)
                {
-                       MEM_freeN(node->arcs);
+                       if (BLI_findindex(&lower_rg->arcs, arc->link_up) == -1)
+                       {
+                               printf("missing arc %p\n", arc->link_up);
+                       }
                }
                
-               node->arcs = MEM_callocN((node->degree + 1) * sizeof(ReebArc*), "adjacency list");
+               
+               verifyMultiResolutionLinks(lower_rg);
        }
+#endif
+}
+/***************************************** BUCKET UTILS **********************************************/
 
-       for(arc = rg->arcs.first; arc; arc= arc->next)
+void addVertToBucket(EmbedBucket *b, float co[3])
+{
+       b->nv++;
+       VecLerpf(b->p, b->p, co, 1.0f / b->nv);
+}
+
+void removeVertFromBucket(EmbedBucket *b, float co[3])
+{
+       VecMulf(b->p, (float)b->nv);
+       VecSubf(b->p, b->p, co);
+       b->nv--;
+       VecMulf(b->p, 1.0f / (float)b->nv);
+}
+
+void mergeBuckets(EmbedBucket *bDst, EmbedBucket *bSrc)
+{
+       if (bDst->nv > 0 && bSrc->nv > 0)
+       {
+               bDst->nv += bSrc->nv;
+               VecLerpf(bDst->p, bDst->p, bSrc->p, (float)bSrc->nv / (float)(bDst->nv));
+       }
+       else if (bSrc->nv > 0)
        {
-               addArcToNodeAdjacencyList(arc->v1, arc);
-               addArcToNodeAdjacencyList(arc->v2, arc);
+               bDst->nv = bSrc->nv;
+               VECCOPY(bDst->p, bSrc->p);
        }
 }
 
-int hasAdjacencyList(ReebGraph *rg)
+void mergeArcBuckets(ReebArc *aDst, ReebArc *aSrc, float start, float end)
 {
-       ReebNode *node;
+       if (aDst->bcount > 0 && aSrc->bcount > 0)
+       {
+               int indexDst = 0, indexSrc = 0;
+               
+               start = MAX3(start, aDst->buckets[0].val, aSrc->buckets[0].val);
+               
+               while(indexDst < aDst->bcount && aDst->buckets[indexDst].val < start)
+               {
+                       indexDst++;
+               }
+
+               while(indexSrc < aSrc->bcount && aSrc->buckets[indexSrc].val < start)
+               {
+                       indexSrc++;
+               }
+               
+               for( ;  indexDst < aDst->bcount &&
+                               indexSrc < aSrc->bcount &&
+                               aDst->buckets[indexDst].val <= end &&
+                               aSrc->buckets[indexSrc].val <= end
+                               
+                        ;      indexDst++, indexSrc++)
+               {
+                       mergeBuckets(aDst->buckets + indexDst, aSrc->buckets + indexSrc);
+               }
+       }
+}
+
+void flipArcBuckets(ReebArc *arc)
+{
+       int i, j;
        
-       for(node = rg->nodes.first; node; node = node->next)
+       for (i = 0, j = arc->bcount - 1; i < j; i++, j--)
+       {
+               EmbedBucket tmp;
+               
+               tmp = arc->buckets[i];
+               arc->buckets[i] = arc->buckets[j];
+               arc->buckets[j] = tmp;
+       }
+}
+
+void allocArcBuckets(ReebArc *arc)
+{
+       int i;
+       float start = ceil(((ReebNode*)arc->head)->weight);
+       arc->bcount = (int)(floor(((ReebNode*)arc->tail)->weight) - start) + 1;
+       
+       if (arc->bcount > 0)
        {
-               if (node->arcs == NULL)
+               arc->buckets = MEM_callocN(sizeof(EmbedBucket) * arc->bcount, "embed bucket");
+               
+               for(i = 0; i < arc->bcount; i++)
                {
-                       return 0;
+                       arc->buckets[i].val = start + i;
                }
        }
+       else
+       {
+               arc->buckets = NULL;
+       }
        
-       return 1;
 }
 
-int countConnectedArcs(ReebGraph *rg, ReebNode *node)
+void resizeArcBuckets(ReebArc *arc)
+{
+       EmbedBucket *oldBuckets = arc->buckets;
+       int oldBCount = arc->bcount;
+       
+       allocArcBuckets(arc);
+       
+       if (oldBCount != 0 && arc->bcount != 0)
+       {
+               int oldStart = (int)oldBuckets[0].val;
+               int oldEnd = (int)oldBuckets[oldBCount - 1].val;
+               int newStart = (int)arc->buckets[0].val;
+               int newEnd = (int)arc->buckets[arc->bcount - 1].val;
+               int oldOffset = 0;
+               int newOffset = 0;
+               int len;
+               
+               if (oldStart < newStart)
+               {
+                       oldOffset = newStart - oldStart;
+               }
+               else
+               {
+                       newOffset = oldStart - newStart;
+               }
+               
+               len = MIN2(oldEnd - (oldStart + oldOffset) + 1, newEnd - (newStart - newOffset) + 1);
+               
+               memcpy(arc->buckets + newOffset, oldBuckets + oldOffset, len * sizeof(EmbedBucket)); 
+       }
+
+       if (oldBuckets != NULL)
+       {
+               MEM_freeN(oldBuckets);
+       }
+}
+
+void calculateArcLength(ReebArc *arc)
 {
-       int count = 0;
+       ReebArcIterator iter;
+       EmbedBucket *bucket = NULL;
+       float *vec0, *vec1;
+
+       arc->length = 0;
+       
+       initArcIterator(&iter, arc, arc->head);
+
+       bucket = nextBucket(&iter);
+       
+       vec0 = arc->head->p;
+       vec1 = arc->head->p; /* in case there's no embedding */
+       
+       while (bucket != NULL)
+       {
+               vec1 = bucket->p;
+               
+               arc->length += VecLenf(vec0, vec1);
+               
+               vec0 = vec1;
+               bucket = nextBucket(&iter);
+       }
+       
+       arc->length += VecLenf(arc->tail->p, vec1);     
+}
+
+void calculateGraphLength(ReebGraph *rg)
+{
+       ReebArc *arc;
        
-       /* use adjacency list if present */
-       if (node->arcs)
+       for (arc = rg->arcs.first; arc; arc = arc->next)
        {
-               ReebArc **arcs;
+               calculateArcLength(arc);
+       }
+}
+
+/**************************************** SYMMETRY HANDLING ******************************************/
+
+void REEB_RadialSymmetry(BNode* root_node, RadialArc* ring, int count)
+{
+       ReebNode *node = (ReebNode*)root_node;
+       float axis[3];
+       int i;
        
-               for(arcs = node->arcs; *arcs; arcs++)
+       VECCOPY(axis, root_node->symmetry_axis);
+       
+       /* first pass, merge incrementally */
+       for (i = 0; i < count - 1; i++)
+       {
+               ReebNode *node1, *node2;
+               ReebArc *arc1, *arc2;
+               float tangent[3];
+               float normal[3];
+               int j = i + 1;
+
+               VecAddf(tangent, ring[i].n, ring[j].n);
+               Crossf(normal, tangent, axis);
+               
+               node1 = (ReebNode*)BLI_otherNode(ring[i].arc, root_node);
+               node2 = (ReebNode*)BLI_otherNode(ring[j].arc, root_node);
+               
+               arc1 = (ReebArc*)ring[i].arc;
+               arc2 = (ReebArc*)ring[j].arc;
+
+               /* mirror first node and mix with the second */
+               BLI_mirrorAlongAxis(node1->p, root_node->p, normal);
+               VecLerpf(node2->p, node2->p, node1->p, 1.0f / (j + 1));
+               
+               /* Merge buckets
+                * there shouldn't be any null arcs here, but just to be safe 
+                * */
+               if (arc1->bcount > 0 && arc2->bcount > 0)
                {
-                       count++;
+                       ReebArcIterator iter1, iter2;
+                       EmbedBucket *bucket1 = NULL, *bucket2 = NULL;
+                       
+                       initArcIterator(&iter1, arc1, (ReebNode*)root_node);
+                       initArcIterator(&iter2, arc2, (ReebNode*)root_node);
+                       
+                       bucket1 = nextBucket(&iter1);
+                       bucket2 = nextBucket(&iter2);
+               
+                       /* Make sure they both start at the same value */       
+                       while(bucket1 && bucket1->val < bucket2->val)
+                       {
+                               bucket1 = nextBucket(&iter1);
+                       }
+                       
+                       while(bucket2 && bucket2->val < bucket1->val)
+                       {
+                               bucket2 = nextBucket(&iter2);
+                       }
+       
+       
+                       for ( ;bucket1 && bucket2; bucket1 = nextBucket(&iter1), bucket2 = nextBucket(&iter2))
+                       {
+                               bucket2->nv += bucket1->nv; /* add counts */
+                               
+                               /* mirror on axis */
+                               BLI_mirrorAlongAxis(bucket1->p, root_node->p, normal);
+                               /* add bucket2 in bucket1 */
+                               VecLerpf(bucket2->p, bucket2->p, bucket1->p, (float)bucket1->nv / (float)(bucket2->nv));
+                       }
                }
        }
-       else
+       
+       /* second pass, mirror back on previous arcs */
+       for (i = count - 1; i > 0; i--)
        {
-               ReebArc *arc;
-               for(arc = rg->arcs.first; arc; arc = arc->next)
+               ReebNode *node1, *node2;
+               ReebArc *arc1, *arc2;
+               float tangent[3];
+               float normal[3];
+               int j = i - 1;
+
+               VecAddf(tangent, ring[i].n, ring[j].n);
+               Crossf(normal, tangent, axis);
+               
+               node1 = (ReebNode*)BLI_otherNode(ring[i].arc, root_node);
+               node2 = (ReebNode*)BLI_otherNode(ring[j].arc, root_node);
+               
+               arc1 = (ReebArc*)ring[i].arc;
+               arc2 = (ReebArc*)ring[j].arc;
+
+               /* copy first node than mirror */
+               VECCOPY(node2->p, node1->p);
+               BLI_mirrorAlongAxis(node2->p, root_node->p, normal);
+               
+               /* Copy buckets
+                * there shouldn't be any null arcs here, but just to be safe 
+                * */
+               if (arc1->bcount > 0 && arc2->bcount > 0)
                {
-                       if (arc->v1 == node || arc->v2 == node)
+                       ReebArcIterator iter1, iter2;
+                       EmbedBucket *bucket1 = NULL, *bucket2 = NULL;
+                       
+                       initArcIterator(&iter1, arc1, node);
+                       initArcIterator(&iter2, arc2, node);
+                       
+                       bucket1 = nextBucket(&iter1);
+                       bucket2 = nextBucket(&iter2);
+               
+                       /* Make sure they both start at the same value */       
+                       while(bucket1 && bucket1->val < bucket2->val)
                        {
-                               count++;
+                               bucket1 = nextBucket(&iter1);
+                       }
+                       
+                       while(bucket2 && bucket2->val < bucket1->val)
+                       {
+                               bucket2 = nextBucket(&iter2);
+                       }
+       
+       
+                       for ( ;bucket1 && bucket2; bucket1 = nextBucket(&iter1), bucket2 = nextBucket(&iter2))
+                       {
+                               /* copy and mirror back to bucket2 */                   
+                               bucket2->nv = bucket1->nv;
+                               VECCOPY(bucket2->p, bucket1->p);
+                               BLI_mirrorAlongAxis(bucket2->p, node->p, normal);
                        }
                }
        }
+}
+
+void REEB_AxialSymmetry(BNode* root_node, BNode* node1, BNode* node2, struct BArc* barc1, BArc* barc2)
+{
+       ReebArc *arc1, *arc2;
+       float nor[3], p[3];
+
+       arc1 = (ReebArc*)barc1;
+       arc2 = (ReebArc*)barc2;
+
+       VECCOPY(nor, root_node->symmetry_axis);
+       
+       /* mirror node2 along axis */
+       VECCOPY(p, node2->p);
+       BLI_mirrorAlongAxis(p, root_node->p, nor);
+
+       /* average with node1 */
+       VecAddf(node1->p, node1->p, p);
+       VecMulf(node1->p, 0.5f);
+       
+       /* mirror back on node2 */
+       VECCOPY(node2->p, node1->p);
+       BLI_mirrorAlongAxis(node2->p, root_node->p, nor);
+       
+       /* Merge buckets
+        * there shouldn't be any null arcs here, but just to be safe 
+        * */
+       if (arc1->bcount > 0 && arc2->bcount > 0)
+       {
+               ReebArcIterator iter1, iter2;
+               EmbedBucket *bucket1 = NULL, *bucket2 = NULL;
+               
+               initArcIterator(&iter1, arc1, (ReebNode*)root_node);
+               initArcIterator(&iter2, arc2, (ReebNode*)root_node);
+               
+               bucket1 = nextBucket(&iter1);
+               bucket2 = nextBucket(&iter2);
        
-       return count;
+               /* Make sure they both start at the same value */       
+               while(bucket1 && bucket1->val < bucket2->val)
+               {
+                       bucket1 = nextBucket(&iter1);
+               }
+               
+               while(bucket2 && bucket2->val < bucket1->val)
+               {
+                       bucket2 = nextBucket(&iter2);
+               }
+
+
+               for ( ;bucket1 && bucket2; bucket1 = nextBucket(&iter1), bucket2 = nextBucket(&iter2))
+               {
+                       bucket1->nv += bucket2->nv; /* add counts */
+                       
+                       /* mirror on axis */
+                       BLI_mirrorAlongAxis(bucket2->p, root_node->p, nor);
+                       /* add bucket2 in bucket1 */
+                       VecLerpf(bucket1->p, bucket1->p, bucket2->p, (float)bucket2->nv / (float)(bucket1->nv));
+
+                       /* copy and mirror back to bucket2 */                   
+                       bucket2->nv = bucket1->nv;
+                       VECCOPY(bucket2->p, bucket1->p);
+                       BLI_mirrorAlongAxis(bucket2->p, root_node->p, nor);
+               }
+       }
 }
 
+/************************************** ADJACENCY LIST *************************************************/
+
+
 /****************************************** SMOOTHING **************************************************/
 
 void postprocessGraph(ReebGraph *rg, char mode)
@@ -492,12 +945,14 @@ int compareArcsWeight(void *varc1, void *varc2)
 {
        ReebArc *arc1 = (ReebArc*)varc1;
        ReebArc *arc2 = (ReebArc*)varc2;
+       ReebNode *node1 = (ReebNode*)arc1->head; 
+       ReebNode *node2 = (ReebNode*)arc2->head; 
        
-       if (arc1->v1->weight < arc2->v1->weight)
+       if (node1->weight < node2->weight)
        {
                return -1;
        }
-       if (arc1->v1->weight > arc2->v1->weight)
+       if (node1->weight > node2->weight)
        {
                return 1;
        }
@@ -514,12 +969,24 @@ void sortArcs(ReebGraph *rg)
 
 /****************************************** FILTERING **************************************************/
 
+float lengthArc(ReebArc *arc)
+{
+#if 0
+       ReebNode *head = (ReebNode*)arc->head;
+       ReebNode *tail = (ReebNode*)arc->tail;
+       
+       return tail->weight - head->weight;
+#else
+       return arc->length;
+#endif
+}
+
 int compareArcs(void *varc1, void *varc2)
 {
        ReebArc *arc1 = (ReebArc*)varc1;
        ReebArc *arc2 = (ReebArc*)varc2;
-       float len1 = arc1->v2->weight - arc1->v1->weight;
-       float len2 = arc2->v2->weight - arc2->v1->weight;
+       float len1 = lengthArc(arc1);
+       float len2 = lengthArc(arc2);
        
        if (len1 < len2)
        {
@@ -539,12 +1006,17 @@ void filterArc(ReebGraph *rg, ReebNode *newNode, ReebNode *removedNode, ReebArc
 {
        ReebArc *arc = NULL, *nextArc = NULL;
 
-       /* first pass, merge buckets for arcs that spawned the two nodes into the source arc*/
-       for(arc = rg->arcs.first; arc; arc = arc->next)
+       if (merging)
        {
-               if (arc->v1 == srcArc->v1 && arc->v2 == srcArc->v2 && arc != srcArc)
+               /* first pass, merge buckets for arcs that spawned the two nodes into the source arc*/
+               for(arc = rg->arcs.first; arc; arc = arc->next)
                {
-                       mergeArcBuckets(srcArc, arc, srcArc->v1->weight, srcArc->v2->weight);
+                       if (arc->head == srcArc->head && arc->tail == srcArc->tail && arc != srcArc)
+                       {
+                               ReebNode *head = srcArc->head;
+                               ReebNode *tail = srcArc->tail;
+                               mergeArcBuckets(srcArc, arc, head->weight, tail->weight);
+                       }
                }
        }
 
@@ -554,48 +1026,57 @@ void filterArc(ReebGraph *rg, ReebNode *newNode, ReebNode *removedNode, ReebArc
        {
                nextArc = arc->next;
                
-               if (arc->v1 == removedNode || arc->v2 == removedNode)
+               if (arc->head == removedNode || arc->tail == removedNode)
                {
-                       if (arc->v1 == removedNode)
+                       if (arc->head == removedNode)
                        {
-                               arc->v1 = newNode;
+                               arc->head = newNode;
                        }
                        else
                        {
-                               arc->v2 = newNode;
+                               arc->tail = newNode;
                        }
 
                        // Remove looped arcs                   
-                       if (arc->v1 == arc->v2)
+                       if (arc->head == arc->tail)
                        {
                                // v1 or v2 was already newNode, since we're removing an arc, decrement degree
-                               newNode->degree--;
+                               NodeDegreeDecrement(rg, newNode);
                                
-                               // If it's safeArc, it'll be removed later, so keep it for now
+                               // If it's srcArc, it'll be removed later, so keep it for now
                                if (arc != srcArc)
                                {
                                        BLI_remlink(&rg->arcs, arc);
-                                       freeArc(arc);
+                                       REEB_freeArc((BArc*)arc);
                                }
                        }
-                       // Remove flipped arcs
-                       else if (arc->v1->weight > arc->v2->weight)
-                       {
-                               // Decrement degree from the other node
-                               OTHER_NODE(arc, newNode)->degree--;
-                               
-                               BLI_remlink(&rg->arcs, arc);
-                               freeArc(arc);
-                       }
                        else
                        {
-                               newNode->degree++; // incrementing degree since we're adding an arc
+                               /* flip arcs that flipped, can happen on diamond shapes, mostly on null arcs */
+                               if (arc->head->weight > arc->tail->weight)
+                               {
+                                       ReebNode *tmp;
+                                       tmp = arc->head;
+                                       arc->head = arc->tail;
+                                       arc->tail = tmp;
+                                       
+                                       flipArcBuckets(arc);
+                               }
+                               //newNode->degree++; // incrementing degree since we're adding an arc
+                               NodeDegreeIncrement(rg, newNode);
+                               mergeArcFaces(rg, arc, srcArc);
 
                                if (merging)
                                {
+                                       ReebNode *head = arc->head;
+                                       ReebNode *tail = arc->tail;
+
                                        // resize bucket list
                                        resizeArcBuckets(arc);
-                                       mergeArcBuckets(arc, srcArc, arc->v1->weight, arc->v2->weight);
+                                       mergeArcBuckets(arc, srcArc, head->weight, tail->weight);
+                                       
+                                       /* update length */
+                                       arc->length += srcArc->length;
                                }
                        }
                }
@@ -615,13 +1096,12 @@ void filterNullReebGraph(ReebGraph *rg)
                // Only collapse arcs too short to have any embed bucket
                if (arc->bcount == 0)
                {
-                       ReebNode *newNode = arc->v1;
-                       ReebNode *removedNode = arc->v2;
+                       ReebNode *newNode = (ReebNode*)arc->head;
+                       ReebNode *removedNode = (ReebNode*)arc->tail;
                        float blend;
                        
                        blend = (float)newNode->degree / (float)(newNode->degree + removedNode->degree); // blending factors
                        
-                       //newNode->weight = FloatLerpf(newNode->weight, removedNode->weight, blend);
                        VecLerpf(newNode->p, newNode->p, removedNode->p, blend);
                        
                        filterArc(rg, newNode, removedNode, arc, 0);
@@ -630,9 +1110,9 @@ void filterNullReebGraph(ReebGraph *rg)
                        nextArc = arc->next;
                        
                        BLI_remlink(&rg->arcs, arc);
-                       freeArc(arc);
+                       REEB_freeArc((BArc*)arc);
                        
-                       BLI_freelinkN(&rg->nodes, removedNode);
+                       BLI_removeNode((BGraph*)rg, (BNode*)removedNode);
                }
                
                arc = nextArc;
@@ -652,22 +1132,28 @@ int filterInternalReebGraph(ReebGraph *rg, float threshold)
                nextArc = arc->next;
 
                // Only collapse non-terminal arcs that are shorter than threshold
-               if ((arc->v1->degree > 1 && arc->v2->degree > 1 && arc->v2->weight - arc->v1->weight < threshold))
+               if (arc->head->degree > 1 && arc->tail->degree > 1 && (lengthArc(arc) < threshold))
                {
                        ReebNode *newNode = NULL;
                        ReebNode *removedNode = NULL;
                        
+#if 0 // Old method
                        /* Keep the node with the highestn number of connected arcs */
-                       if (arc->v1->degree >= arc->v2->degree)
+                       if (arc->head->degree >= arc->tail->degree)
                        {
-                               newNode = arc->v1;
-                               removedNode = arc->v2;
+                               newNode = arc->head;
+                               removedNode = arc->tail;
                        }
                        else
                        {
-                               newNode = arc->v2;
-                               removedNode = arc->v1;
+                               newNode = arc->tail;
+                               removedNode = arc->head;
                        }
+#else
+                       /* Always remove lower node, so arcs don't flip */
+                       newNode = arc->head;
+                       removedNode = arc->tail;
+#endif
                        
                        filterArc(rg, newNode, removedNode, arc, 1);
 
@@ -675,9 +1161,9 @@ int filterInternalReebGraph(ReebGraph *rg, float threshold)
                        nextArc = arc->next;
                        
                        BLI_remlink(&rg->arcs, arc);
-                       freeArc(arc);
+                       REEB_freeArc((BArc*)arc);
                        
-                       BLI_freelinkN(&rg->nodes, removedNode);
+                       BLI_removeNode((BGraph*)rg, (BNode*)removedNode);
                        value = 1;
                }
                
@@ -694,220 +1180,388 @@ int filterExternalReebGraph(ReebGraph *rg, float threshold)
        
        BLI_sortlist(&rg->arcs, compareArcs);
 
-       arc = rg->arcs.first;
-       while(arc)
+       
+       for (arc = rg->arcs.first; arc; arc = nextArc)
        {
                nextArc = arc->next;
 
                // Only collapse terminal arcs that are shorter than threshold
-               if ((arc->v1->degree == 1 || arc->v2->degree == 1) && arc->v2->weight - arc->v1->weight < threshold)
+               if ((arc->head->degree == 1 || arc->tail->degree == 1) && (lengthArc(arc) < threshold))
                {
                        ReebNode *terminalNode = NULL;
                        ReebNode *middleNode = NULL;
-                       ReebNode *newNode = NULL;
                        ReebNode *removedNode = NULL;
-                       int merging = 0;
                        
                        // Assign terminal and middle nodes
-                       if (arc->v1->degree == 1)
+                       if (arc->head->degree == 1)
                        {
-                               terminalNode = arc->v1;
-                               middleNode = arc->v2;
+                               terminalNode = arc->head;
+                               middleNode = arc->tail;
                        }
                        else
                        {
-                               terminalNode = arc->v2;
-                               middleNode = arc->v1;
+                               terminalNode = arc->tail;
+                               middleNode = arc->head;
                        }
                        
-                       // If middle node is a normal node, merge to terminal node
+                       // If middle node is a normal node, it will be removed later
                        if (middleNode->degree == 2)
                        {
-                               merging = 1;
-                               newNode = terminalNode;
-                               removedNode = middleNode;
+                               continue;
+//                             removedNode = middleNode;
+//
+//                             filterArc(rg, terminalNode, removedNode, arc, 1);
                        }
                        // Otherwise, just plain remove of the arc
                        else
                        {
-                               merging = 0;
-                               newNode = middleNode;
                                removedNode = terminalNode;
-                       }
-                       
-                       // Merging arc
-                       if (merging)
-                       {
-                               filterArc(rg, newNode, removedNode, arc, 1);
-                       }
-                       else
-                       {
+
                                // removing arc, so we need to decrease the degree of the remaining node
-                               newNode->degree--;
+                               NodeDegreeDecrement(rg, middleNode);
                        }
 
                        // Reset nextArc, it might have changed
                        nextArc = arc->next;
                        
                        BLI_remlink(&rg->arcs, arc);
-                       freeArc(arc);
+                       REEB_freeArc((BArc*)arc);
                        
-                       BLI_freelinkN(&rg->nodes, removedNode);
+                       BLI_removeNode((BGraph*)rg, (BNode*)removedNode);
                        value = 1;
                }
-               
-               arc = nextArc;
        }
+
+       /* join on normal nodes */      
+       removeNormalNodes(rg);
        
        return value;
 }
 
-/************************************** WEIGHT SPREADING ***********************************************/
-
-int compareVerts( const void* a, const void* b )
+int filterCyclesReebGraph(ReebGraph *rg, float distance_threshold)
 {
-       EditVert *va = *(EditVert**)a;
-       EditVert *vb = *(EditVert**)b;
-       int value = 0;
+       int filtered = 0;
        
-       if (va->tmp.fp < vb->tmp.fp)
-       {
-               value = -1;
-       }
-       else if (va->tmp.fp > vb->tmp.fp)
+       if (BLI_isGraphCyclic((BGraph*)rg))
        {
-               value = 1;
-       }
-
-       return value;           
-}
+               ReebArc *arc1, *arc2;
+               ReebArc *next2;
+               
+               for (arc1 = rg->arcs.first; arc1; arc1 = arc1->next)
+               {
+                       for (arc2 = rg->arcs.first; arc2; arc2 = next2)
+                       {
+                               next2 = arc2->next;
+                               if (arc1 != arc2 && arc1->head == arc2->head && arc1->tail == arc2->tail)
+                               {
+                                       mergeArcEdges(rg, arc1, arc2, MERGE_APPEND);
+                                       mergeArcFaces(rg, arc1, arc2);
+                                       mergeArcBuckets(arc1, arc2, arc1->head->weight, arc1->tail->weight);
 
-void spreadWeight(EditMesh *em)
-{
-       EditVert **verts, *eve;
-       float lastWeight = 0.0f;
-       int totvert = BLI_countlist(&em->verts);
-       int i;
-       int work_needed = 1;
+                                       NodeDegreeDecrement(rg, arc1->head);
+                                       NodeDegreeDecrement(rg, arc1->tail);
+
+                                       BLI_remlink(&rg->arcs, arc2);
+                                       REEB_freeArc((BArc*)arc2);
+                                       
+                                       filtered = 1;
+                               }
+                       }
+               }
+       }
        
-       verts = MEM_callocN(sizeof(EditVert*) * totvert, "verts array");
+       return filtered;
+}
+
+int filterSmartReebGraph(ReebGraph *rg, float threshold)
+{
+       ReebArc *arc = NULL, *nextArc = NULL;
+       int value = 0;
        
-       for(eve = em->verts.first, i = 0; eve; eve = eve->next, i++)
-       {
-               verts[i] = eve;
+       BLI_sortlist(&rg->arcs, compareArcs);
+
+#ifdef DEBUG_REEB
+       {       
+               EditFace *efa;
+               for(efa=G.editMesh->faces.first; efa; efa=efa->next) {
+                       efa->tmp.fp = -1;
+               }
        }
-       
-       while(work_needed == 1)
+#endif
+
+       arc = rg->arcs.first;
+       while(arc)
        {
-               work_needed = 0;
-               qsort(verts, totvert, sizeof(EditVert*), compareVerts);
+               nextArc = arc->next;
                
-               for(i = 0; i < totvert; i++)
+               /* need correct normals and center */
+               recalc_editnormals();
+
+               // Only test terminal arcs
+               if (arc->head->degree == 1 || arc->tail->degree == 1)
                {
-                       eve = verts[i];
+                       GHashIterator ghi;
+                       int merging = 0;
+                       int total = BLI_ghash_size(arc->faces);
+                       float avg_angle = 0;
+                       float avg_vec[3] = {0,0,0};
                        
-                       if (i == 0 || (eve->tmp.fp - lastWeight) > FLT_EPSILON)
+                       for(BLI_ghashIterator_init(&ghi, arc->faces);
+                               !BLI_ghashIterator_isDone(&ghi);
+                               BLI_ghashIterator_step(&ghi))
                        {
-                               lastWeight = eve->tmp.fp;
+                               EditFace *efa = BLI_ghashIterator_getValue(&ghi);
+
+#if 0
+                               ReebArcIterator iter;
+                               EmbedBucket *bucket = NULL;
+                               EmbedBucket *previous = NULL;
+                               float min_distance = -1;
+                               float angle = 0;
+               
+                               initArcIterator(&iter, arc, arc->head);
+               
+                               bucket = nextBucket(&iter);
+                               
+                               while (bucket != NULL)
+                               {
+                                       float *vec0 = NULL;
+                                       float *vec1 = bucket->p;
+                                       float midpoint[3], tangent[3];
+                                       float distance;
+               
+                                       /* first bucket. Previous is head */
+                                       if (previous == NULL)
+                                       {
+                                               vec0 = arc->head->p;
+                                       }
+                                       /* Previous is a valid bucket */
+                                       else
+                                       {
+                                               vec0 = previous->p;
+                                       }
+                                       
+                                       VECCOPY(midpoint, vec1);
+                                       
+                                       distance = VecLenf(midpoint, efa->cent);
+                                       
+                                       if (min_distance == -1 || distance < min_distance)
+                                       {
+                                               min_distance = distance;
+                                       
+                                               VecSubf(tangent, vec1, vec0);
+                                               Normalize(tangent);
+                                               
+                                               angle = Inpf(tangent, efa->n);
+                                       }
+                                       
+                                       previous = bucket;
+                                       bucket = nextBucket(&iter);
+                               }
+                               
+                               avg_angle += saacos(fabs(angle));
+#ifdef DEBUG_REEB
+                               efa->tmp.fp = saacos(fabs(angle));
+#endif
+#else
+                               VecAddf(avg_vec, avg_vec, efa->n);              
+#endif
                        }
-                       else
+
+
+#if 0                  
+                       avg_angle /= total;
+#else
+                       VecMulf(avg_vec, 1.0 / total);
+                       avg_angle = Inpf(avg_vec, avg_vec);
+#endif
+                       
+                       arc->angle = avg_angle;
+                       
+                       if (avg_angle > threshold)
+                               merging = 1;
+                       
+                       if (merging)
                        {
-                               work_needed = 1;
-                               eve->tmp.fp = lastWeight + FLT_EPSILON * 2;
-                               lastWeight = eve->tmp.fp;
+                               ReebNode *terminalNode = NULL;
+                               ReebNode *middleNode = NULL;
+                               ReebNode *newNode = NULL;
+                               ReebNode *removedNode = NULL;
+                               int merging = 0;
+                               
+                               // Assign terminal and middle nodes
+                               if (arc->head->degree == 1)
+                               {
+                                       terminalNode = arc->head;
+                                       middleNode = arc->tail;
+                               }
+                               else
+                               {
+                                       terminalNode = arc->tail;
+                                       middleNode = arc->head;
+                               }
+                               
+                               // If middle node is a normal node, merge to terminal node
+                               if (middleNode->degree == 2)
+                               {
+                                       merging = 1;
+                                       newNode = terminalNode;
+                                       removedNode = middleNode;
+                               }
+                               // Otherwise, just plain remove of the arc
+                               else
+                               {
+                                       merging = 0;
+                                       newNode = middleNode;
+                                       removedNode = terminalNode;
+                               }
+                               
+                               // Merging arc
+                               if (merging)
+                               {
+                                       filterArc(rg, newNode, removedNode, arc, 1);
+                               }
+                               else
+                               {
+                                       // removing arc, so we need to decrease the degree of the remaining node
+                                       //newNode->degree--;
+                                       NodeDegreeDecrement(rg, newNode);
+                               }
+       
+                               // Reset nextArc, it might have changed
+                               nextArc = arc->next;
+                               
+                               BLI_remlink(&rg->arcs, arc);
+                               REEB_freeArc((BArc*)arc);
+                               
+                               BLI_freelinkN(&rg->nodes, removedNode);
+                               value = 1;
                        }
                }
+               
+               arc = nextArc;
        }
        
-       MEM_freeN(verts);
+       return value;
 }
-/*********************************** GRAPH AS TREE FUNCTIONS *******************************************/
 
-int subtreeDepth(ReebNode *node, ReebArc *rootArc)
+void filterGraph(ReebGraph *rg, short options, float threshold_internal, float threshold_external)
 {
-       int depth = 0;
+       int done = 1;
        
-       /* Base case, no arcs leading away */
-       if (node->arcs == NULL || *(node->arcs) == NULL)
-       {
-               return 0;
-       }
-       else
-       {
-               ReebArc ** pArc;
+       calculateGraphLength(rg);
+
+       verifyNodeDegree(rg);
 
-               for(pArc = node->arcs; *pArc; pArc++)
+       /* filter until there's nothing more to do */
+       while (done == 1)
+       {
+               done = 0; /* no work done yet */
+               
+               if (options & SKGEN_FILTER_EXTERNAL)
                {
-                       ReebArc *arc = *pArc;
-                       
-                       /* only arcs that go down the tree */
-                       if (arc != rootArc)
-                       {
-                               ReebNode *newNode = OTHER_NODE(arc, node);
-                               depth = MAX2(depth, subtreeDepth(newNode, arc));
-                       }
+//                     done |= filterExternalReebGraph(rg, threshold_external * rg->resolution);
+                       done |= filterExternalReebGraph(rg, threshold_external);
+                       verifyNodeDegree(rg);
+               }
+       
+               if (options & SKGEN_FILTER_INTERNAL)
+               {
+//                     done |= filterInternalReebGraph(rg, threshold_internal * rg->resolution);
+                       done |= filterInternalReebGraph(rg, threshold_internal);
+                       verifyNodeDegree(rg);
                }
        }
+
+       if (options & SKGEN_FILTER_SMART)
+       {
+               filterSmartReebGraph(rg, 0.5);
+               BLI_rebuildAdjacencyList((BGraph*)rg);
+               filterCyclesReebGraph(rg, 0.5);
+       }
        
-       return depth + 1;
-}
+       verifyNodeDegree(rg);
 
-/*************************************** CYCLE DETECTION ***********************************************/
+       repositionNodes(rg);
 
-int detectCycle(ReebNode *node, ReebArc *srcArc)
+       /* Filtering might have created degree 2 nodes, so remove them */
+       removeNormalNodes(rg);
+}
+
+void finalizeGraph(ReebGraph *rg, char passes, char method)
 {
-       int value = 0;
+       int i;
+       
+       BLI_rebuildAdjacencyList((BGraph*)rg);
+
+       sortNodes(rg);
+       
+       sortArcs(rg);
        
-       if (node->flags == 0)
+       for(i = 0; i <  passes; i++)
        {
-               ReebArc ** pArc;
+               postprocessGraph(rg, method);
+       }
+}
 
-               /* mark node as visited */
-               node->flags = 1;
+/************************************** WEIGHT SPREADING ***********************************************/
 
-               for(pArc = node->arcs; *pArc && value == 0; pArc++)
-               {
-                       ReebArc *arc = *pArc;
-                       
-                       /* don't go back on the source arc */
-                       if (arc != srcArc)
-                       {
-                               value = detectCycle(OTHER_NODE(arc, node), arc);
-                       }
-               }
+int compareVerts( const void* a, const void* b )
+{
+       EditVert *va = *(EditVert**)a;
+       EditVert *vb = *(EditVert**)b;
+       int value = 0;
+       
+       if (va->tmp.fp < vb->tmp.fp)
+       {
+               value = -1;
        }
-       else
+       else if (va->tmp.fp > vb->tmp.fp)
        {
                value = 1;
        }
-       
-       return value;
+
+       return value;           
 }
 
-int    isGraphCyclic(ReebGraph *rg)
+void spreadWeight(EditMesh *em)
 {
-       ReebNode *node;
-       int value = 0;
+       EditVert **verts, *eve;
+       float lastWeight = 0.0f;
+       int totvert = BLI_countlist(&em->verts);
+       int i;
+       int work_needed = 1;
        
-       /* NEED TO CHECK IF ADJACENCY LIST EXIST */
+       verts = MEM_callocN(sizeof(EditVert*) * totvert, "verts array");
        
-       /* Mark all nodes as not visited */
-       for(node = rg->nodes.first; node; node = node->next)
+       for(eve = em->verts.first, i = 0; eve; eve = eve->next, i++)
        {
-               node->flags = 0;
+               verts[i] = eve;
        }
-
-       /* detectCycles in subgraphs */ 
-       for(node = rg->nodes.first; node && value == 0; node = node->next)
+       
+       while(work_needed == 1)
        {
-               /* only for nodes in subgraphs that haven't been visited yet */
-               if (node->flags == 0)
+               work_needed = 0;
+               qsort(verts, totvert, sizeof(EditVert*), compareVerts);
+               
+               for(i = 0; i < totvert; i++)
                {
-                       value = value || detectCycle(node, NULL);
-               }               
+                       eve = verts[i];
+                       
+                       if (i == 0 || (eve->tmp.fp - lastWeight) > FLT_EPSILON)
+                       {
+                               lastWeight = eve->tmp.fp;
+                       }
+                       else
+                       {
+                               work_needed = 1;
+                               eve->tmp.fp = lastWeight + FLT_EPSILON * 2;
+                               lastWeight = eve->tmp.fp;
+                       }
+               }
        }
        
-       return value;
+       MEM_freeN(verts);
 }
 
 /******************************************** EXPORT ***************************************************/
@@ -917,9 +1571,8 @@ void exportNode(FILE *f, char *text, ReebNode *node)
        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]);
 }
 
-void exportGraph(ReebGraph *rg, int count)
+void REEB_exportGraph(ReebGraph *rg, int count)
 {
-#ifdef DEBUG_REEB
        ReebArc *arc;
        char filename[128];
        FILE *f;
@@ -937,19 +1590,23 @@ void exportGraph(ReebGraph *rg, int count)
        for(arc = rg->arcs.first; arc; arc = arc->next)
        {
                int i;
+               float p[3];
                
-               exportNode(f, "v1", arc->v1);
+               exportNode(f, "v1", arc->head);
                
                for(i = 0; i < arc->bcount; i++)
                {
                        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]);
                }
                
-               exportNode(f, "v2", arc->v2);
+               VecAddf(p, arc->tail->p, arc->head->p);
+               VecMulf(p, 0.5f);
+               
+               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));
+               exportNode(f, "v2", arc->tail);
        }       
        
        fclose(f);
-#endif
 }
 
 /***************************************** MAIN ALGORITHM **********************************************/
@@ -960,7 +1617,7 @@ ReebArc * findConnectedArc(ReebGraph *rg, ReebArc *arc, ReebNode *v)
        
        for(nextArc = rg->arcs.first; nextArc; nextArc = nextArc->next)
        {
-               if (arc != nextArc && (nextArc->v1 == v || nextArc->v2 == v))
+               if (arc != nextArc && (nextArc->head == v || nextArc->tail == v))
                {
                        break;
                }
@@ -969,46 +1626,71 @@ ReebArc * findConnectedArc(ReebGraph *rg, ReebArc *arc, ReebNode *v)
        return nextArc;
 }
 
+
 void removeNormalNodes(ReebGraph *rg)
 {
-       ReebArc *arc;
+       ReebArc *arc, *nextArc;
        
        // Merge degree 2 nodes
-       for(arc = rg->arcs.first; arc; arc = arc->next)
+       for(arc = rg->arcs.first; arc; arc = nextArc)
        {
-               while (arc->v1->degree == 2 || arc->v2->degree == 2)
+               nextArc = arc->next;
+               
+               while (arc->head->degree == 2 || arc->tail->degree == 2)
                {
                        // merge at v1
-                       if (arc->v1->degree == 2)
+                       if (arc->head->degree == 2)
                        {
-                               ReebArc *nextArc = findConnectedArc(rg, arc, arc->v1);
+                               ReebArc *connectedArc = (ReebArc*)BLI_findConnectedArc((BGraph*)rg, (BArc*)arc, (BNode*)arc->head);
 
-                               // Merge arc only if needed
-                               if (arc->v1 == nextArc->v2)
-                               {                               
-                                       mergeConnectedArcs(rg, arc, nextArc);
+                               /* If arcs are one after the other */
+                               if (arc->head == connectedArc->tail)
+                               {               
+                                       /* remove furthest arc */               
+                                       if (arc->tail->weight < connectedArc->head->weight)
+                                       {
+                                               mergeConnectedArcs(rg, arc, connectedArc);
+                                               nextArc = arc->next;
+                                       }
+                                       else
+                                       {
+                                               mergeConnectedArcs(rg, connectedArc, arc);
+                                               break; /* arc was removed, move to next */
+                                       }
                                }
-                               // Otherwise, mark down vert
+                               /* Otherwise, arcs are side by side */
                                else
                                {
-                                       arc->v1->degree = 3;
+                                       /* Don't do anything, we need to keep the lowest node, even if degree 2 */
+                                       break;
                                }
                        }
                        
                        // merge at v2
-                       if (arc->v2->degree == 2)
+                       if (arc->tail->degree == 2)
                        {
-                               ReebArc *nextArc = findConnectedArc(rg, arc, arc->v2);
+                               ReebArc *connectedArc = (ReebArc*)BLI_findConnectedArc((BGraph*)rg, (BArc*)arc, (BNode*)arc->tail);
                                
-                               // Merge arc only if needed
-                               if (arc->v2 == nextArc->v1)
+                               /* If arcs are one after the other */
+                               if (arc->tail == connectedArc->head)
                                {                               
-                                       mergeConnectedArcs(rg, arc, nextArc);
+                                       /* remove furthest arc */               
+                                       if (arc->head->weight < connectedArc->tail->weight)
+                                       {
+                                               mergeConnectedArcs(rg, arc, connectedArc);
+                                               nextArc = arc->next;
+                                       }
+                                       else
+                                       {
+                                               mergeConnectedArcs(rg, connectedArc, arc);
+                                               break; /* arc was removed, move to next */
+                                       }
                                }
-                               // Otherwise, mark down vert
+                               /* Otherwise, arcs are side by side */
                                else
                                {
-                                       arc->v2->degree = 3;
+                                       /* Don't do anything, we need to keep the lowest node, even if degree 2 */
+                                       break;
                                }
                        }
                }
@@ -1041,11 +1723,23 @@ ReebArc *nextArcMappedToEdge(ReebArc *arc, ReebEdge *e)
        return result;
 }
 
-typedef enum {
-       MERGE_LOWER,
-       MERGE_HIGHER,
-       MERGE_APPEND
-} MergeDirection;
+void addFacetoArc(ReebArc *arc, EditFace *efa)
+{
+       BLI_ghash_insert(arc->faces, efa, efa);
+}
+
+void mergeArcFaces(ReebGraph *rg, ReebArc *aDst, ReebArc *aSrc)
+{
+       GHashIterator ghi;
+       
+       for(BLI_ghashIterator_init(&ghi, aSrc->faces);
+               !BLI_ghashIterator_isDone(&ghi);
+               BLI_ghashIterator_step(&ghi))
+       {
+               EditFace *efa = BLI_ghashIterator_getValue(&ghi);
+               BLI_ghash_insert(aDst->faces, efa, efa);
+       }
+} 
 
 void mergeArcEdges(ReebGraph *rg, ReebArc *aDst, ReebArc *aSrc, MergeDirection direction)
 {
@@ -1109,29 +1803,32 @@ int mergeConnectedArcs(ReebGraph *rg, ReebArc *a0, ReebArc *a1)
        int result = 0;
        ReebNode *removedNode = NULL;
        
+       a0->length += a1->length;
+       
        mergeArcEdges(rg, a0, a1, MERGE_APPEND);
+       mergeArcFaces(rg, a0, a1);
        
        // Bring a0 to the combine length of both arcs
-       if (a0->v2 == a1->v1)
+       if (a0->tail == a1->head)
        {
-               removedNode = a0->v2;
-               a0->v2 = a1->v2;
+               removedNode = a0->tail;
+               a0->tail = a1->tail;
        }
-       else if (a0->v1 == a1->v2)
+       else if (a0->head == a1->tail)
        {
-               removedNode = a0->v1;
-               a0->v1 = a1->v1;
+               removedNode = a0->head;
+               a0->head = a1->head;
        }
        
        resizeArcBuckets(a0);
        // Merge a1 in a0
-       mergeArcBuckets(a0, a1, a0->v1->weight, a0->v2->weight);
+       mergeArcBuckets(a0, a1, a0->head->weight, a0->tail->weight);
        
        // remove a1 from graph
        BLI_remlink(&rg->arcs, a1);
-       freeArc(a1);
+       REEB_freeArc((BArc*)a1);
        
-       BLI_freelinkN(&rg->nodes, removedNode);
+       BLI_removeNode((BGraph*)rg, (BNode*)removedNode);
        result = 1;
        
        return result;
@@ -1141,74 +1838,89 @@ int mergeArcs(ReebGraph *rg, ReebArc *a0, ReebArc *a1)
 {
        int result = 0;
        // TRIANGLE POINTS DOWN
-       if (a0->v1->weight == a1->v1->weight) // heads are the same
+       if (a0->head->weight == a1->head->weight) // heads are the same
        {
-               if (a0->v2->weight == a1->v2->weight) // tails also the same, arcs can be totally merge together
+               if (a0->tail->weight == a1->tail->weight) // tails also the same, arcs can be totally merge together
                {
                        mergeArcEdges(rg, a0, a1, MERGE_APPEND);
+                       mergeArcFaces(rg, a0, a1);
                        
-                       mergeArcBuckets(a0, a1, a0->v1->weight, a0->v2->weight);
+                       mergeArcBuckets(a0, a1, a0->head->weight, a0->tail->weight);
                        
                        // Adjust node degree
-                       a1->v1->degree--;
-                       a1->v2->degree--;
+                       //a1->head->degree--;
+                       NodeDegreeDecrement(rg, a1->head);
+                       //a1->tail->degree--;
+                       NodeDegreeDecrement(rg, a1->tail);
                        
                        // remove a1 from graph
                        BLI_remlink(&rg->arcs, a1);
                        
-                       freeArc(a1);
+                       REEB_freeArc((BArc*)a1);
                        result = 1;
                }
-               else if (a0->v2->weight > a1->v2->weight) // a1->v2->weight is in the middle
+               else if (a0->tail->weight > a1->tail->weight) // a1->tail->weight is in the middle
                {
                        mergeArcEdges(rg, a1, a0, MERGE_LOWER);
+                       mergeArcFaces(rg, a1, a0);
 
                        // Adjust node degree
-                       a0->v1->degree--;
-                       a1->v2->degree++;
+                       //a0->head->degree--;
+                       NodeDegreeDecrement(rg, a0->head);
+                       //a1->tail->degree++;
+                       NodeDegreeIncrement(rg, a1->tail);
                        
-                       mergeArcBuckets(a1, a0, a1->v1->weight, a1->v2->weight);
-                       a0->v1 = a1->v2;
+                       mergeArcBuckets(a1, a0, a1->head->weight, a1->tail->weight);
+                       a0->head = a1->tail;
                        resizeArcBuckets(a0);
                }
                else // a0>n2 is in the middle
                {
                        mergeArcEdges(rg, a0, a1, MERGE_LOWER);
+                       mergeArcFaces(rg, a0, a1);
                        
                        // Adjust node degree
-                       a1->v1->degree--;
-                       a0->v2->degree++;
+                       //a1->head->degree--;
+                       NodeDegreeDecrement(rg, a1->head);
+                       //a0->tail->degree++;
+                       NodeDegreeIncrement(rg, a0->tail);
                        
-                       mergeArcBuckets(a0, a1, a0->v1->weight, a0->v2->weight);
-                       a1->v1 = a0->v2;
+                       mergeArcBuckets(a0, a1, a0->head->weight, a0->tail->weight);
+                       a1->head = a0->tail;
                        resizeArcBuckets(a1);
                }
        }
        // TRIANGLE POINTS UP
-       else if (a0->v2->weight == a1->v2->weight) // tails are the same
+       else if (a0->tail->weight == a1->tail->weight) // tails are the same
        {
-               if (a0->v1->weight > a1->v1->weight) // a0->v1->weight is in the middle
+               if (a0->head->weight > a1->head->weight) // a0->head->weight is in the middle
                {
                        mergeArcEdges(rg, a0, a1, MERGE_HIGHER);
+                       mergeArcFaces(rg, a0, a1);
                        
                        // Adjust node degree
-                       a1->v2->degree--;
-                       a0->v1->degree++;
+                       //a1->tail->degree--;
+                       NodeDegreeDecrement(rg, a1->tail);
+                       //a0->head->degree++;
+                       NodeDegreeIncrement(rg, a0->head);
                        
-                       mergeArcBuckets(a0, a1, a0->v1->weight, a0->v2->weight);
-                       a1->v2 = a0->v1;
+                       mergeArcBuckets(a0, a1, a0->head->weight, a0->tail->weight);
+                       a1->tail = a0->head;
                        resizeArcBuckets(a1);
                }
-               else // a1->v1->weight is in the middle
+               else // a1->head->weight is in the middle
                {
                        mergeArcEdges(rg, a1, a0, MERGE_HIGHER);
+                       mergeArcFaces(rg, a1, a0);
 
                        // Adjust node degree
-                       a0->v2->degree--;
-                       a1->v1->degree++;
+                       //a0->tail->degree--;
+                       NodeDegreeDecrement(rg, a0->tail);
+                       //a1->head->degree++;
+                       NodeDegreeIncrement(rg, a1->head);
 
-                       mergeArcBuckets(a1, a0, a1->v1->weight, a1->v2->weight);
-                       a0->v2 = a1->v1;
+                       mergeArcBuckets(a1, a0, a1->head->weight, a1->tail->weight);
+                       a0->tail = a1->head;
                        resizeArcBuckets(a0);
                }
        }
@@ -1229,7 +1941,7 @@ void glueByMergeSort(ReebGraph *rg, ReebArc *a0, ReebArc *a1, ReebEdge *e0, Reeb
                
                if (total == 0) // if it wasn't a total merge, go forward
                {
-                       if (a0->v2->weight < a1->v2->weight)
+                       if (a0->tail->weight < a1->tail->weight)
                        {
                                a0 = nextArcMappedToEdge(a0, e0);
                        }
@@ -1252,25 +1964,6 @@ void mergePaths(ReebGraph *rg, ReebEdge *e0, ReebEdge *e1, ReebEdge *e2)
        glueByMergeSort(rg, a0, a2, e0, e2);
 } 
 
-ReebNode * addNode(ReebGraph *rg, EditVert *eve, float weight)
-{
-       ReebNode *node = NULL;
-       
-       node = MEM_callocN(sizeof(ReebNode), "reeb node");
-       
-       node->flags = 0; // clear flags on init
-       node->arcs = NULL;
-       node->degree = 0;
-       node->weight = weight;
-       node->index = rg->totnodes;
-       VECCOPY(node->p, eve->co);      
-       
-       BLI_addtail(&rg->nodes, node);
-       rg->totnodes++;
-       
-       return node;
-}
-
 ReebEdge * createArc(ReebGraph *rg, ReebNode *node1, ReebNode *node2)
 {
        ReebEdge *edge;
@@ -1288,7 +1981,9 @@ ReebEdge * createArc(ReebGraph *rg, ReebNode *node1, ReebNode *node2)
                arc = MEM_callocN(sizeof(ReebArc), "reeb arc");
                edge = MEM_callocN(sizeof(ReebEdge), "reeb edge");
                
-               arc->flags = 0; // clear flags on init
+               arc->flag = 0; // clear flag on init
+               arc->symmetry_level = 0;
+               arc->faces = BLI_ghash_new(BLI_ghashutil_ptrhash, BLI_ghashutil_ptrcmp);
                
                if (node1->weight <= node2->weight)
                {
@@ -1301,12 +1996,14 @@ ReebEdge * createArc(ReebGraph *rg, ReebNode *node1, ReebNode *node2)
                        v2 = node1;     
                }
                
-               arc->v1 = v1;
-               arc->v2 = v2;
+               arc->head = v1;
+               arc->tail = v2;
                
                // increase node degree
-               v1->degree++;
-               v2->degree++;
+               //v1->degree++;
+               NodeDegreeIncrement(rg, v1);
+               //v2->degree++;
+               NodeDegreeIncrement(rg, v2);
 
                BLI_edgehash_insert(rg->emap, node1->index, node2->index, edge);
                
@@ -1321,8 +2018,8 @@ ReebEdge * createArc(ReebGraph *rg, ReebNode *node1, ReebNode *node2)
                /* adding buckets for embedding */
                allocArcBuckets(arc);
                
-               offset = arc->v1->weight;
-               len = arc->v2->weight - arc->v1->weight;
+               offset = arc->head->weight;
+               len = arc->tail->weight - arc->head->weight;
 
 #if 0
                /* This is the actual embedding filling described in the paper
@@ -1330,8 +2027,8 @@ ReebEdge * createArc(ReebGraph *rg, ReebNode *node1, ReebNode *node2)
                 */
                if (arc->bcount > 0)
                {
-                       addVertToBucket(&(arc->buckets[0]), arc->v1->co);
-                       addVertToBucket(&(arc->buckets[arc->bcount - 1]), arc->v2->co);
+                       addVertToBucket(&(arc->buckets[0]), arc->head->co);
+                       addVertToBucket(&(arc->buckets[arc->bcount - 1]), arc->tail->co);
                }
 #else
                for(i = 0; i < arc->bcount; i++)
@@ -1349,7 +2046,7 @@ ReebEdge * createArc(ReebGraph *rg, ReebNode *node1, ReebNode *node2)
        return edge;
 }
 
-void addTriangleToGraph(ReebGraph *rg, ReebNode * n1, ReebNode * n2, ReebNode * n3)
+void addTriangleToGraph(ReebGraph *rg, ReebNode * n1, ReebNode * n2, ReebNode * n3, EditFace *efa)
 {
        ReebEdge *re1, *re2, *re3;
        ReebEdge *e1, *e2, *e3;
@@ -1359,6 +2056,10 @@ void addTriangleToGraph(ReebGraph *rg, ReebNode * n1, ReebNode * n2, ReebNode *
        re2 = createArc(rg, n2, n3);
        re3 = createArc(rg, n3, n1);
        
+       addFacetoArc(re1->arc, efa);
+       addFacetoArc(re2->arc, efa);
+       addFacetoArc(re3->arc, efa);
+       
        len1 = (float)fabs(n1->weight - n2->weight);
        len2 = (float)fabs(n2->weight - n3->weight);
        len3 = (float)fabs(n3->weight - n1->weight);
@@ -1412,10 +2113,9 @@ ReebGraph * generateReebGraph(EditMesh *em, int subdivisions)
        int countfaces = 0;
 #endif
        
-       rg = MEM_callocN(sizeof(ReebGraph), "reeb graph");
+       rg = newReebGraph();
        
-       rg->totnodes = 0;
-       rg->emap = BLI_edgehash_new();
+       rg->resolution = subdivisions;
        
        totvert = BLI_countlist(&em->verts);
        totfaces = BLI_countlist(&em->faces);
@@ -1425,14 +2125,18 @@ ReebGraph * generateReebGraph(EditMesh *em, int subdivisions)
        /* Spread weight to minimize errors */
        spreadWeight(em);
 
-       renormalizeWeight(em, (float)subdivisions);
+       renormalizeWeight(em, (float)rg->resolution);
 
        /* Adding vertice */
-       for(index = 0, eve = em->verts.first; eve; index++, eve = eve->next)
+       for(index = 0, eve = em->verts.first; eve; eve = eve->next)
        {
-               eve->hash = index;
-               eve->f2 = 0;
-               eve->tmp.p = addNode(rg, eve, eve->tmp.fp);
+               if (eve->h == 0)
+               {
+                       eve->hash = index;
+                       eve->f2 = 0;
+                       eve->tmp.p = addNode(rg, eve, eve->tmp.fp);
+                       index++;
+               }
        }
        
        /* Temporarely convert node list to dynamic list, for indexed access */
@@ -1441,30 +2145,35 @@ ReebGraph * generateReebGraph(EditMesh *em, int subdivisions)
        /* Adding face, edge per edge */
        for(efa = em->faces.first; efa; efa = efa->next)
        {
-               ReebNode *n1, *n2, *n3;
-               
-               n1 = (ReebNode*)BLI_dlist_find_link(dlist, efa->v1->hash);
-               n2 = (ReebNode*)BLI_dlist_find_link(dlist, efa->v2->hash);
-               n3 = (ReebNode*)BLI_dlist_find_link(dlist, efa->v3->hash);
-               
-               addTriangleToGraph(rg, n1, n2, n3);
-               
-               if (efa->v4)
-               {
-                       ReebNode *n4 = (ReebNode*)efa->v4->tmp.p;
-                       addTriangleToGraph(rg, n1, n3, n4);
-               }
-
-#ifdef DEBUG_REEB
-               countfaces++;
-               if (countfaces % 100 == 0)
+               if (efa->h == 0)
                {
-                       printf("face %i of %i\n", countfaces, totfaces);
-               }
+                       ReebNode *n1, *n2, *n3;
+                       
+                       n1 = (ReebNode*)BLI_dlist_find_link(dlist, efa->v1->hash);
+                       n2 = (ReebNode*)BLI_dlist_find_link(dlist, efa->v2->hash);
+                       n3 = (ReebNode*)BLI_dlist_find_link(dlist, efa->v3->hash);
+                       
+                       addTriangleToGraph(rg, n1, n2, n3, efa);
+                       
+                       if (efa->v4)
+                       {
+                               ReebNode *n4 = (ReebNode*)efa->v4->tmp.p;
+                               addTriangleToGraph(rg, n1, n3, n4, efa);
+                       }
+#ifdef DEBUG_REEB
+                       countfaces++;
+                       if (countfaces % 100 == 0)
+                       {
+                               printf("\rface %i of %i", countfaces, totfaces);
+                               verifyFaces(rg);
+                       }
 #endif
-               
-               
+               }
        }
+       
+       printf("\n");
+       
+       
        BLI_listbase_from_dlist(dlist, &rg->nodes);
        
        removeNormalNodes(rg);
@@ -1534,6 +2243,35 @@ static float cotan_weight(float *v1, float *v2, float *v3)
        return Inpf(a, b)/clen;
 }
 
+void addTriangle(EditVert *v1, EditVert *v2, EditVert *v3, long e1, long e2, long e3)
+{
+       /* Angle opposite e1 */
+       float t1= cotan_weight(v1->co, v2->co, v3->co) / e2;
+       
+       /* Angle opposite e2 */
+       float t2 = cotan_weight(v2->co, v3->co, v1->co) / e3;
+
+       /* Angle opposite e3 */
+       float t3 = cotan_weight(v3->co, v1->co, v2->co) / e1;
+       
+       int i1 = v1->hash;
+       int i2 = v2->hash;
+       int i3 = v3->hash;
+       
+       nlMatrixAdd(i1, i1, t2+t3);
+       nlMatrixAdd(i2, i2, t1+t3);
+       nlMatrixAdd(i3, i3, t1+t2);
+
+       nlMatrixAdd(i1, i2, -t3);
+       nlMatrixAdd(i2, i1, -t3);
+
+       nlMatrixAdd(i2, i3, -t1);
+       nlMatrixAdd(i3, i2, -t1);
+
+       nlMatrixAdd(i3, i1, -t2);
+       nlMatrixAdd(i1, i3, -t2);
+}
+
 int weightToHarmonic(EditMesh *em)
 {
        NLboolean success;
@@ -1561,49 +2299,55 @@ int weightToHarmonic(EditMesh *em)
        /* Find local extrema */
        for(index = 0, eve = em->verts.first; eve; index++, eve = eve->next)
        {
-               EditEdge *eed;
-               int maximum = 1;
-               int minimum = 1;
-               
-               eve->hash = index; /* Assign index to vertex */
-               
-               NextEdgeForVert(NULL, NULL); /* Reset next edge */
-               for(eed = NextEdgeForVert(em, eve); eed && (maximum || minimum); eed = NextEdgeForVert(em, eve))
+               if (eve->h == 0)
                {
-                       EditVert *eve2;
+                       EditEdge *eed;
+                       int maximum = 1;
+                       int minimum = 1;
                        
-                       if (eed->v1 == eve)
-                       {
-                               eve2 = eed->v2;
-                       }
-                       else
+                       eve->hash = index; /* Assign index to vertex */
+                       
+                       NextEdgeForVert(NULL, NULL); /* Reset next edge */
+                       for(eed = NextEdgeForVert(em, eve); eed && (maximum || minimum); eed = NextEdgeForVert(em, eve))
                        {
-                               eve2 = eed->v1;
+                               EditVert *eve2;
+                               
+                               if (eed->v1 == eve)
+                               {
+                                       eve2 = eed->v2;
+                               }
+                               else
+                               {
+                                       eve2 = eed->v1;
+                               }
+                               
+                               if (eve2->h == 0)
+                               {
+                                       /* Adjacent vertex is bigger, not a local maximum */
+                                       if (eve2->tmp.fp > eve->tmp.fp)
+                                       {
+                                               maximum = 0;
+                                       }
+                                       /* Adjacent vertex is smaller, not a local minimum */
+                                       else if (eve2->tmp.fp < eve->tmp.fp)
+                                       {
+                                               minimum = 0;
+                                       }
+                               }
                        }
                        
-                       /* Adjacent vertex is bigger, not a local maximum */
-                       if (eve2->tmp.fp > eve->tmp.fp)
+                       if (maximum || minimum)
                        {
-                               maximum = 0;
+                               float w = eve->tmp.fp;
+                               eve->f1 = 0;
+                               nlSetVariable(0, index, w);
+                               nlLockVariable(index);
                        }
-                       /* Adjacent vertex is smaller, not a local minimum */
-                       else if (eve2->tmp.fp < eve->tmp.fp)
+                       else
                        {
-                               minimum = 0;
+                               eve->f1 = 1;
                        }
                }
-               
-               if (maximum || minimum)
-               {
-                       float w = eve->tmp.fp;
-                       eve->f1 = 0;
-                       nlSetVariable(0, index, w);
-                       nlLockVariable(index);
-               }
-               else
-               {
-                       eve->f1 = 1;
-               }
        }
        
        nlBegin(NL_MATRIX);
@@ -1617,39 +2361,34 @@ int weightToHarmonic(EditMesh *em)
        /* Add faces count to the edge weight */
        for(efa = em->faces.first; efa; efa = efa->next)
        {
-               efa->e1->tmp.l++;
-               efa->e2->tmp.l++;
-               efa->e3->tmp.l++;
+               if (efa->h == 0)
+               {
+                       efa->e1->tmp.l++;
+                       efa->e2->tmp.l++;
+                       efa->e3->tmp.l++;
+                       
+                       if (efa->e4)
+                       {
+                               efa->e4->tmp.l++;
+                       }
+               }
        }
 
        /* Add faces angle to the edge weight */
        for(efa = em->faces.first; efa; efa = efa->next)
        {
-               /* Angle opposite e1 */
-               float t1= cotan_weight(efa->v1->co, efa->v2->co, efa->v3->co) / efa->e2->tmp.l;
-               
-               /* Angle opposite e2 */
-               float t2 = cotan_weight(efa->v2->co, efa->v3->co, efa->v1->co) / efa->e3->tmp.l;
-
-               /* Angle opposite e3 */
-               float t3 = cotan_weight(efa->v3->co, efa->v1->co, efa->v2->co) / efa->e1->tmp.l;
-               
-               int i1 = efa->v1->hash;
-               int i2 = efa->v2->hash;
-               int i3 = efa->v3->hash;
-               
-               nlMatrixAdd(i1, i1, t2+t3);
-               nlMatrixAdd(i2, i2, t1+t3);
-               nlMatrixAdd(i3, i3, t1+t2);
-       
-               nlMatrixAdd(i1, i2, -t3);
-               nlMatrixAdd(i2, i1, -t3);
-       
-               nlMatrixAdd(i2, i3, -t1);
-               nlMatrixAdd(i3, i2, -t1);
-       
-               nlMatrixAdd(i3, i1, -t2);
-               nlMatrixAdd(i1, i3, -t2);
+               if (efa->h == 0)
+               {
+                       if (efa->v4 == NULL)
+                       {
+                               addTriangle(efa->v1, efa->v2, efa->v3, efa->e1->tmp.l, efa->e2->tmp.l, efa->e3->tmp.l);
+                       }
+                       else
+                       {
+                               addTriangle(efa->v1, efa->v2, efa->v3, efa->e1->tmp.l, efa->e2->tmp.l, 2);
+                               addTriangle(efa->v3, efa->v4, efa->v1, efa->e3->tmp.l, efa->e4->tmp.l, 2);
+                       }
+               }
        }
        
        nlEnd(NL_MATRIX);
@@ -1701,7 +2440,7 @@ EditEdge * NextEdgeForVert(EditMesh *em, EditVert *v)
 
        for( ; e ; e = e->next)
        {
-               if (e->v1 == v || e->v2 == v)
+               if ((e->v1 == v || e->v2 == v) && (e->h == 0))
                {
                        break;
                }
@@ -1728,9 +2467,10 @@ int weightFromDistance(EditMesh *em)
                return 0;
        }
 
-       /* Initialize vertice flags and find at least one selected vertex */
-       for(eve = em->verts.first; eve && vCount == 0; eve = eve->next)
+       /* Initialize vertice flag and find at least one selected vertex */
+       for(eve = em->verts.first; eve; eve = eve->next)
        {
+               eve->tmp.fp = 0;
                eve->f1 = 0;
                if (eve->f & SELECT)
                {
@@ -1762,7 +2502,7 @@ int weightFromDistance(EditMesh *em)
                                        
                                        edges = MEM_callocN(totedge * sizeof(EditEdge*), "Edges");
                                        
-                                       /* Calculate edge weight and initialize edge flags */
+                                       /* Calculate edge weight and initialize edge flag */
                                        for(eed= em->edges.first; eed; eed= eed->next)
                                        {
                                                eed->tmp.fp = VecLenf(eed->v1->co, eed->v2->co);
@@ -1826,28 +2566,42 @@ int weightFromDistance(EditMesh *em)
                                                        }                               
                                                        current_eve->tmp.fp = currentWeight;
                                                }
+                                               
+                                       printf("\redge %i / %i", eIndex, totedge);
+                                               
                                        } while (select_eed != NULL);
                                        
                                        MEM_freeN(edges);
+                                       
+                                       printf("\n");
                                }
                        }
                }
        }
 
+       for(eve = em->verts.first; eve && vCount == 0; eve = eve->next)
+       {
+               if (eve->f1 == 0)
+               {
+                       printf("vertex not reached\n");
+                       break;
+               }
+       }
+
        return 1;
 }
 
-MCol MColFromWeight(EditVert *eve)
+MCol MColFromVal(float val)
 {
        MCol col;
        col.a = 255;
-       col.b = (char)(eve->tmp.fp * 255);
+       col.b = (char)(val * 255);
        col.g = 0;
-       col.r = (char)((1.0f - eve->tmp.fp) * 255);
+       col.r = (char)((1.0f - val) * 255);
        return col;
 }
 
-void weightToVCol(EditMesh *em)
+void weightToVCol(EditMesh *em, int index)
 {
        EditFace *efa;
        MCol *mcol;
@@ -1855,15 +2609,149 @@ void weightToVCol(EditMesh *em)
                return;
        }
        
+       for(efa=em->faces.first; efa; efa=efa->next) {
+               mcol = CustomData_em_get_n(&em->fdata, efa->data, CD_MCOL, index);
+
+               if (mcol)
+               {                               
+                       mcol[0] = MColFromVal(efa->v1->tmp.fp);
+                       mcol[1] = MColFromVal(efa->v2->tmp.fp);
+                       mcol[2] = MColFromVal(efa->v3->tmp.fp);
+       
+                       if(efa->v4) {
+                               mcol[3] = MColFromVal(efa->v4->tmp.fp);
+                       }
+               }
+       }
+}
+
+void angleToVCol(EditMesh *em, int index)
+{
+       EditFace *efa;
+       MCol *mcol;
+
+       if (!EM_vertColorCheck()) {
+               return;
+       }
+       
+       for(efa=em->faces.first; efa; efa=efa->next) {
+               MCol col;
+               if (efa->tmp.fp > 0)
+               {
+                       col = MColFromVal(efa->tmp.fp / (M_PI / 2 + 0.1));
+               }
+               else
+               {
+                       col.a = 255;
+                       col.r = 0;
+                       col.g = 255;
+                       col.b = 0;
+               }
+
+               mcol = CustomData_em_get_n(&em->fdata, efa->data, CD_MCOL, index);
+                               
+               if (mcol)
+               {
+                       mcol[0] = col;
+                       mcol[1] = col;
+                       mcol[2] = col;
+       
+                       if(efa->v4) {
+                               mcol[3] = col;
+                       }
+               }
+       }
+}
+
+void blendColor(MCol *dst, MCol *src)
+{
+#if 1
+       float blend_src = (float)src->a / (float)(src->a + dst->a);
+       float blend_dst = (float)dst->a / (float)(src->a + dst->a);
+       dst->a += src->a;
+       dst->r = (char)(dst->r * blend_dst + src->r * blend_src);
+       dst->g = (char)(dst->g * blend_dst + src->g * blend_src);
+       dst->b = (char)(dst->b * blend_dst + src->b * blend_src);
+#else
+       dst->r = src->r;
+       dst->g = src->g;
+       dst->b = src->b;
+#endif
+}
+
+void arcToVCol(ReebGraph *rg, EditMesh *em, int index)
+{
+       GHashIterator ghi;
+       EditFace *efa;
+       ReebArc *arc;
+       MCol *mcol;
+       MCol col;
+       int total = BLI_countlist(&rg->arcs);
+       int i = 0;
+
+       if (!EM_vertColorCheck()) {
+               return;
+       }
+       
+       col.a = 0;
+       
+       col.r = 0;
+       col.g = 0;
+       col.b = 0;
+
+       for(efa=em->faces.first; efa; efa=efa->next) {
+               mcol = CustomData_em_get_n(&em->fdata, efa->data, CD_MCOL, index);
+               
+               if (mcol)
+               {
+                       mcol[0] = col;
+                       mcol[1] = col;
+                       mcol[2] = col;
+       
+                       if(efa->v4) {
+                               mcol[3] = col;
+                       }
+               }
+       }
+
+       for (arc = rg->arcs.first; arc; arc = arc->next, i++)
+       {
+               float r,g,b;
+               col.a = 1;
+               
+               hsv_to_rgb((float)i / (float)total, 1, 1, &r, &g, &b);
+               
+               col.r = FTOCHAR(r);
+               col.g = FTOCHAR(g);
+               col.b = FTOCHAR(b);
+               
+               for(BLI_ghashIterator_init(&ghi, arc->faces);
+                       !BLI_ghashIterator_isDone(&ghi);
+                       BLI_ghashIterator_step(&ghi))
+               {
+                       efa = BLI_ghashIterator_getValue(&ghi);
+
+                       mcol = CustomData_em_get(&em->fdata, efa->data, CD_MCOL);
+                                       
+                       blendColor(&mcol[0], &col);
+                       blendColor(&mcol[1], &col);
+                       blendColor(&mcol[2], &col);
+       
+                       if(efa->v4) {
+                               blendColor(&mcol[3], &col);
+                       }
+               }
+       }
+
        for(efa=em->faces.first; efa; efa=efa->next) {
                mcol = CustomData_em_get(&em->fdata, efa->data, CD_MCOL);
               &n