Skeleton retargetting - Preliminary commit. Results are encouraging but nothing *that* useful yet

Smarter heuristic noise arc filtering for Reeb graph 

8 files changed, 2990 insertions, 664 deletions

diff --git a/source/blender/include/BIF_editarmature.h b/source/blender/include/BIF_editarmature.h
index 0e1557ac378..5ed02471978 100644
--- a/source/blender/include/BIF_editarmature.h
+++ b/source/blender/include/BIF_editarmature.h
@@ -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
 
 
diff --git a/source/blender/include/butspace.h b/source/blender/include/butspace.h
index 7571d64be91..bd1cc46c211 100644
--- a/source/blender/include/butspace.h
+++ b/source/blender/include/butspace.h
@@ -444,7 +444,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
diff --git a/source/blender/include/reeb.h b/source/blender/include/reeb.h
index c8352aedec5..e26b4080249 100644
--- a/source/blender/include/reeb.h
+++ b/source/blender/include/reeb.h
@@ -30,6 +30,7 @@
 
 #include "DNA_listBase.h"
 
+struct GHash;
 struct EdgeHash;
 struct ReebArc;
 struct ReebEdge;
@@ -55,7 +56,11 @@ typedef struct ReebNode {
 	int degree;
 	float weight;
 	float p[3];
-	int flags;
+	int flag;
+
+	int symmetry_level;
+	int symmetry_flag;
+	float symmetry_axis[3];
 } ReebNode;
 
 typedef struct ReebEdge {
@@ -63,6 +68,7 @@ typedef struct ReebEdge {
 	struct ReebArc  *arc;
 	struct ReebNode *v1, *v2;
 	struct ReebEdge *nextEdge;
+	int flag;
 } ReebEdge;
 
 typedef struct ReebArc {
@@ -71,7 +77,13 @@ typedef struct ReebArc {
 	struct ReebNode *v1, *v2;
 	struct EmbedBucket *buckets;
 	int	bcount;
-	int flags;
+	int flag;
+
+	int symmetry_level;
+	int symmetry_flag;
+
+	struct GHash *faces;	
+	float angle;
 } ReebArc;
 
 typedef struct ReebArcIterator {
@@ -87,21 +99,28 @@ 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(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);
+ReebGraph * newReebGraph();
 
 #define OTHER_NODE(arc, node) ((arc->v1 == node) ? arc->v2 : arc->v1)
 
 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 * 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);
 
@@ -121,7 +140,32 @@ int countConnectedArcs(ReebGraph *rg, ReebNode *node);
 int hasAdjacencyList(ReebGraph *rg); 
 int	isGraphCyclic(ReebGraph *rg);
 
-/* Sanity check */
+/*------------ Symmetry handling ------------*/
+void markdownSymmetry(ReebGraph *rg);
+
+/* ReebNode symmetry flags */
+#define SYM_TOPOLOGICAL	1
+#define SYM_PHYSICAL	2
+
+/* the following two are exclusive */
+#define SYM_AXIAL		4
+#define SYM_RADIAL		8
+
+/* ReebArc symmetry flags
+ * 
+ * axial symetry sides */
+#define SYM_SIDE_POSITIVE		1
+#define SYM_SIDE_NEGATIVE		2
+
+
+
+/*------------ Sanity check ------------*/
 void verifyBuckets(ReebGraph *rg);
+void verifyFaces(ReebGraph *rg);
+
+/*********************** PUBLIC *********************************/
+ReebGraph *BIF_ReebGraphFromEditMesh(void);
+void REEB_freeGraph(ReebGraph *rg);
+void REEB_exportGraph(ReebGraph *rg, int count);
 
 #endif /*REEB_H_*/
diff --git a/source/blender/makesdna/DNA_scene_types.h b/source/blender/makesdna/DNA_scene_types.h
index 7b1b979b777..1ddebf5d6e4 100644
--- a/source/blender/makesdna/DNA_scene_types.h
+++ b/source/blender/makesdna/DNA_scene_types.h
@@ -433,6 +433,8 @@ 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;
 	short skgen_options;
 	char  skgen_postpro;
 	char  skgen_postpro_passes;
@@ -831,6 +833,7 @@ typedef struct Scene {
 #define	SKGEN_CUT_LENGTH		8
 #define	SKGEN_CUT_ANGLE			16
 #define	SKGEN_CUT_CORRELATION	32
+#define	SKGEN_HARMONIC			64
 
 #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
index 00000000000..d3cea1be786
--- /dev/null
+++ b/source/blender/src/autoarmature.c
@@ -0,0 +1,1690 @@
+/**
+ * $Id: editarmature.c 14848 2008-05-15 08:05:56Z aligorith $
+ *
+ * ***** 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 "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;
+	struct RigNode *head;
+	
+	ReebGraph *link;
+} RigGraph;
+
+typedef struct RigNode {
+	struct RigNode *next, *prev;
+	float p[3];
+	int degree;
+	struct RigArc **arcs;
+	int flag;
+	
+	int symmetry_level;
+	int symmetry_flag;
+	float symmetry_axis[3];
+
+	ReebNode *link;
+} RigNode;
+
+typedef struct RigArc {
+	struct RigArc *next, *prev;
+	RigNode *head, *tail;
+	ListBase edges;
+	float length;
+	int flag;
+
+	int symmetry_level;
+	int symmetry_flag;
+	
+	int count;
+	ReebArc *link;
+} 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);
+void RIG_markdownSymmetry(RigGraph *rg);
+void RIG_markdownSymmetryArc(RigArc *arc, RigNode *node, int level);
+void RIG_markdownSecondarySymmetry(RigNode *node, int depth, int level);
+
+
+/*******************************************************************************************************/
+static RigNode *RIG_otherNode(RigArc *arc, RigNode *node)
+{
+	if (arc->head == node)
+		return arc->tail;
+	else
+		return arc->head;
+}
+
+static void RIG_flagNodes(RigGraph *rg, int flag)
+{
+	RigNode *node;
+	
+	for(node = rg->nodes.first; node; node = node->next)
+	{
+		node->flag = flag;
+	}
+}
+
+static void RIG_flagArcs(RigGraph *rg, int flag)
+{
+	RigArc *arc;
+	
+	for(arc = rg->arcs.first; arc; arc = arc->next)
+	{
+		arc->flag = flag;
+	}
+}
+
+static void RIG_addArcToNodeAdjacencyList(RigNode *node, RigArc *arc)
+{
+	node->arcs[node->degree] = arc;
+	node->degree++;
+}
+/*********************************** 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;
+}
+
+/************************************* ALLOCATORS ******************************************************/
+
+static RigGraph *newRigGraph()
+{
+	RigGraph *rg;
+	rg = MEM_callocN(sizeof(RigGraph), "rig graph");
+	
+	rg->head = NULL;
+	
+	return rg;
+}
+
+static RigArc *newRigArc(RigGraph *rg)
+{
+	RigArc *arc;
+	
+	arc = MEM_callocN(sizeof(RigArc), "rig arc");
+	arc->length = 0;
+	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;
+}
+
+/************************************ DESTRUCTORS ******************************************************/
+
+static void RIG_freeRigNode(RigNode *node)
+{
+	if (node->arcs)
+	{
+		MEM_freeN(node->arcs);
+	}
+}
+
+static void RIG_freeRigArc(RigArc *arc)
+{
+	BLI_freelistN(&arc->edges);
+}
+
+static void RIG_freeRigGraph(RigGraph *rg)
+{
+	RigNode *node;
+	RigArc *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)
+	{
+		RIG_freeRigNode(node);
+	}
+	BLI_freelistN(&rg->nodes);
+	
+	MEM_freeN(rg);
+}
+
+/*******************************************************************************************************/
+
+static void RIG_buildAdjacencyList(RigGraph *rg)
+{
+	RigNode *node;
+	RigArc *arc;
+
+	for(node = rg->nodes.first; node; node = node->next)
+	{
+		if (node->arcs != NULL)
+		{
+			MEM_freeN(node->arcs);
+		}
+		
+		node->arcs = MEM_callocN((node->degree + 1) * sizeof(RigArc*), "adjacency list");
+		
+		/* temporary use to indicate the first index available in the lists */
+		node->degree = 0;
+	}
+
+	for(arc = rg->arcs.first; arc; arc= arc->next)
+	{
+		RIG_addArcToNodeAdjacencyList(arc->head, arc);
+		RIG_addArcToNodeAdjacencyList(arc->tail, arc);
+	}
+}
+
+static void RIG_replaceNode(RigGraph *rg, RigNode *node_src, RigNode *node_replaced)
+{
+	RigArc *arc, *next_arc;
+	
+	for (arc = rg->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;
+			
+			RIG_freeRigArc(arc);
+			BLI_freelinkN(&rg->arcs, arc);
+		}
+	}
+}
+
+static void RIG_removeDoubleNodes(RigGraph *rg, float limit)
+{
+	RigNode *node_src, *node_replaced;
+	
+	for(node_src = rg->nodes.first; node_src; node_src = node_src->next)
+	{
+		for(node_replaced = rg->nodes.first; node_replaced; node_replaced = node_replaced->next)
+		{
+			if (node_replaced != node_src && VecLenf(node_replaced->p, node_src->p) <= limit)
+			{
+				RIG_replaceNode(rg, node_src, node_replaced);
+			}
+		}
+	}
+	
+}
+
+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));
+}
+
+/*********************************** GRAPH AS TREE FUNCTIONS *******************************************/
+
+int RIG_subtreeDepth(RigNode *node, RigArc *rootArc)
+{
+	int depth = 0;
+	
+	/* Base case, no arcs leading away */
+	if (node->arcs == NULL || *(node->arcs) == NULL)
+	{
+		return 0;
+	}
+	else
+	{
+		RigArc ** pArc;
+
+		for(pArc = node->arcs; *pArc; pArc++)
+		{
+			RigArc *arc = *pArc;
+			
+			/* only arcs that go down the tree */
+			if (arc != rootArc)
+			{
+				RigNode *newNode = RIG_otherNode(arc, node);
+				depth = MAX2(depth, RIG_subtreeDepth(newNode, arc));
+			}
+		}
+	}
+	
+	return depth + BLI_countlist(&rootArc->edges);
+}
+
+int RIG_countConnectedArcs(RigGraph *rg, RigNode *node)
+{
+	int count = 0;
+	
+	/* use adjacency list if present */
+	if (node->arcs)
+	{
+		RigArc **arcs;
+	
+		for(arcs = node->arcs; *arcs; arcs++)
+		{
+			count++;
+		}
+	}
+	else
+	{
+		RigArc *arc;
+		for(arc = rg->arcs.first; arc; arc = arc->next)
+		{
+			if (arc->head == node || arc->tail == node)
+			{
+				count++;
+			}
+		}
+	}
+	
+	return count;
+}
+
+/********************************* SYMMETRY DETECTION **************************************************/
+
+static 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
+{
+	RigArc *arc; 
+	float n[3]; /* normalized vector joining the nodes of the arc */
+} RadialArc;
+
+void RIG_markRadialSymmetry(RigNode *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;
+
+	/* mark topological symmetry */
+	node->symmetry_flag |= SYM_TOPOLOGICAL;
+
+	/* count the number of arcs in the symmetry ring */
+	for (i = 0; node->arcs[i] != NULL; i++)
+	{
+		RigArc *connectedArc = node->arcs[i];
+		
+		/* depth is store as a negative in flag. symmetry level is positive */
+		if (connectedArc->symmetry_level == -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++)
+	{
+		RigArc *connectedArc = node->arcs[i];
+		
+		/* depth is store as a negative in flag. symmetry level is positive */
+		if (connectedArc->symmetry_level == -depth)
+		{
+			RigNode *otherNode = RIG_otherNode(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++)
+	{
+		RigNode *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 = RIG_otherNode(ring[i].arc, node);
+		node2 = RIG_otherNode(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)
+	{
+		/* mark node as symmetric physically */
+		VECCOPY(node->symmetry_axis, axis);
+		node->symmetry_flag |= SYM_PHYSICAL;
+		node->symmetry_flag |= SYM_RADIAL;
+	}
+
+	MEM_freeN(ring);
+}
+
+static void setSideAxialSymmetry(RigNode *root_node, RigNode *end_node, RigArc *arc)
+{
+	float vec[3];
+	
+	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;
+	}
+}
+
+void RIG_markAxialSymmetry(RigNode *node, int depth, float axis[3])
+{
+	RigArc *arc1 = NULL;
+	RigArc *arc2 = NULL;
+	RigNode *node1 = NULL, *node2 = NULL;
+	float limit = G.scene->toolsettings->skgen_symmetry_limit;
+	float nor[3], vec[3], p[3];
+	int i;
+	
+	/* mark topological symmetry */
+	node->symmetry_flag |= SYM_TOPOLOGICAL;
+
+	for (i = 0; node->arcs[i] != NULL; i++)
+	{
+		RigArc *connectedArc = 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 = RIG_otherNode(arc1, node);
+			}
+			else
+			{
+				arc2 = connectedArc;
+				node2 = RIG_otherNode(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(vec, node1->p, node->p);
+	Normalize(vec);
+	VecSubf(p, 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);
+	mirrorAlongAxis(p, node->p, nor);
+	
+	/* check if it's within limit before continuing */
+	if (VecLenf(node1->p, p) <= limit)
+	{
+		/* mark node as symmetric physically */
+		VECCOPY(node->symmetry_axis, nor);
+		node->symmetry_flag |= SYM_PHYSICAL;
+		node->symmetry_flag |= SYM_AXIAL;
+
+		/* set side on arcs */
+		setSideAxialSymmetry(node, node1, arc1);
+		setSideAxialSymmetry(node, node2, arc2);
+		printf("flag: %i <-> %i\n", arc1->symmetry_flag, arc2->symmetry_flag);
+	}
+	else
+	{
+		printf("NOT SYMMETRIC!\n");
+		printf("%f <= %f\n", VecLenf(node1->p, p), limit);
+		printvecf("axis", nor);
+	}
+}
+
+void RIG_markdownSecondarySymmetry(RigNode *node, int depth, int level)
+{
+	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; node->arcs[i] != NULL; i++)
+	{
+		RigArc *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)
+	{
+		RIG_markAxialSymmetry(node, depth, axis);
+	}
+	else
+	{
+		RIG_markRadialSymmetry(node, depth, axis);
+	}
+	
+	/* markdown secondary symetries */	
+	for (i = 0; node->arcs[i] != NULL; i++)
+	{
+		RigArc *connectedArc = node->arcs[i];
+		
+		if (connectedArc->symmetry_level == -depth)
+		{
+			/* markdown symmetry for branches corresponding to the depth */
+			RIG_markdownSymmetryArc(connectedArc, node, level + 1);
+		}
+	}
+}
+
+void RIG_markdownSymmetryArc(RigArc *arc, RigNode *node, int level)
+{
+	int i;
+	arc->symmetry_level = level;
+	
+	node = RIG_otherNode(arc, node);
+	
+	for (i = 0; node->arcs[i] != NULL; i++)
+	{
+		RigArc *connectedArc = node->arcs[i];
+		
+		if (connectedArc != arc)
+		{
+			RigNode *connectedNode = RIG_otherNode(connectedArc, node);
+			
+			/* symmetry level is positive value, negative values is subtree depth */
+			connectedArc->symmetry_level = -RIG_subtreeDepth(connectedNode, connectedArc);
+		}
+	}
+
+	arc = NULL;
+
+	for (i = 0; node->arcs[i] != NULL; i++)
+	{
+		int issymmetryAxis = 0;
+		RigArc *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; node->arcs[j] != NULL && issymmetryAxis == 1; j++)
+			{
+				RigArc *otherArc = node->arcs[j];
+				
+				/* different arc, same depth */
+				if (otherArc != connectedArc && otherArc->symmetry_level == connectedArc->symmetry_level)
+				{
+					/* 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)
+	{
+		RIG_markdownSymmetryArc(arc, node, level);
+	}
+
+	
+	/* secondary symmetry */
+	for (i = 0; node->arcs[i] != NULL; i++)
+	{
+		RigArc *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 */
+			RIG_markdownSecondarySymmetry(node, -connectedArc->symmetry_level, level);
+		}
+	}
+}
+
+void RIG_markdownSymmetry(RigGraph *rg)
+{
+	RigNode *node;
+	RigArc *arc;
+	
+	/* mark down all arcs as non-symetric */
+	RIG_flagArcs(rg, 0);
+	
+	/* mark down all nodes as not on the symmetry axis */
+	RIG_flagNodes(rg, 0);
+
+	if (rg->head)
+	{
+		node = rg->head;
+	}
+	else
+	{
+		/* !TODO! DO SOMETHING SMART HERE */
+		return;
+	}
+	
+	/* only work on acyclic graphs and if only one arc is incident on the first node */
+	if (RIG_countConnectedArcs(rg, node) == 1)
+	{
+		arc = node->arcs[0];
+		
+		RIG_markdownSymmetryArc(arc, node, 1);
+
+		/* mark down non-symetric arcs */
+		for (arc = rg->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;
+				}
+			}
+		}
+	}
+}
+
+/*******************************************************************************************************/
+
+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 = arc->head;
+		}
+	}
+}
+
+/*******************************************************************************************************/
+
+static 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);
+	}
+}
+
+static 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");
+}
+
+static void RIG_printArc(RigArc *arc)
+{
+	RigEdge *edge;
+
+	printf("\n");
+
+	RIG_printNode(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(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);
+		}
+	}
+	
+	RIG_removeDoubleNodes(rg, 0);
+	
+	RIG_buildAdjacencyList(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;
+	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;
+	}
+	
+	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 void retargetArctoArcAggresive(RigArc *iarc)
+{
+	ReebArcIterator iter;
+	RigEdge *edge;
+	EmbedBucket *bucket = NULL;
+	ReebNode *node_start, *node_end;
+	ReebArc *earc = iarc->link;
+	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;
+
+	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->v2;
+		node_end = earc->v1;
+	}
+	else
+	{
+		node_start = earc->v1;
+		node_end = earc->v2;
+	}
+	
+	/* init with first values */
+	for (i = 0; i < nb_joints; i++)
+	{
+		positions[i] = 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;
+
+	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;
+				
+				if (i < nb_joints)
+				{
+					bucket = nextNBucket(&iter, positions[i] - last_index);
+					vec2 = bucket->p;
+					vec_cache[i + 1] = vec2; /* update cache for updated position */
+				}
+				else
+				{
+					vec2 = node_end->p;
+				}
+				
+				vec1 = vec_cache[i];
+				
+
+				VecSubf(vec_second, vec2, vec1);
+				length2 = Normalize(vec_second);
+	
+				/* check angle */
+				if (i != 0)
+				{
+					RigEdge *previous = edge->prev;
+					float angle = previous->angle;
+					float test_angle = previous->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 */
+						new_cost += G.scene->toolsettings->skgen_retarget_angle_weight * fabs((test_angle - angle) / test_angle);
+					}
+					else
+					{
+						new_cost += M_PI;
+					}
+				}
+	
+				/* LENGTH COST HERE */
+				new_cost += G.scene->toolsettings->skgen_retarget_length_weight * fabs((length2 - edge->length) / edge->length);
+				cost_cache[i] = new_cost;
+				
+				last_index =  positions[i];
+			}
+			
+			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);
+		}
+	}
+
+	vec0 = node_start->p;
+	initArcIterator(&iter, earc, node_start);
+	
+	printPositions(best_positions, nb_joints);
+	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 = nextNBucket(&iter, best_positions[i] - last_index);
+			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);
+			printf("position: %i\n", best_positions[i]);
+			printf("last_index: %i\n", last_index);
+		}
+		
+		vec0 = vec1;
+		last_index =  best_positions[i];
+	}
+	
+	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;
+	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 = earc->v2;
+		node_end = earc->v1;
+	}
+	else
+	{
+		node_start = earc->v1;
+		node_end = earc->v2;
+	}
+	
+	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;
+	
+	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->v2->p);
+			VECCOPY(bone->tail, earc->v1->p);
+		}
+		/* or not */
+		else
+		{
+			VECCOPY(bone->head, earc->v1->p);
+			VECCOPY(bone->tail, earc->v2->p);
+		}
+	}
+	else
+	{
+		RetargetMode mode = detectArcRetargetMode(iarc);
+		
+		if (mode == RETARGET_AGGRESSIVE)
+		{
+			printf("aggresive\n");
+			retargetArctoArcAggresive(iarc);
+		}
+		else
+		{		
+			retargetArctoArcLength(iarc);
+		}
+	}
+}
+
+static void findCorrespondingArc(RigArc *start_arc, RigNode *start_node, RigArc *next_iarc)
+{
+	ReebNode *enode = start_node->link;
+	ReebArc *next_earc;
+	int symmetry_level = next_iarc->symmetry_level;
+	int symmetry_flag = next_iarc->symmetry_flag;
+	int i;
+	
+	next_iarc->link = NULL;
+		
+	for(i = 0, next_earc = enode->arcs[i]; next_earc; i++, next_earc = enode->arcs[i])
+	{
+		if (next_earc->flag == 0 && /* not already taken */
+			next_earc->symmetry_flag == symmetry_flag &&
+			next_earc->symmetry_level == symmetry_level)
+		{
+			printf("-----------------------\n");
+			printf("CORRESPONDING ARC FOUND\n");
+			RIG_printArcBones(next_iarc);
+
+			next_earc->flag = 1; // mark as taken
+			next_iarc->link = next_earc;
+			break;
+		}
+	}
+	
+	if (next_iarc->link == 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, next_earc = enode->arcs[i]; next_earc; i++, next_earc = 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;
+	RigNode *inode = start_node;
+	ReebNode *enode = start_node->link;
+	RigArc *next_iarc;
+	int i;
+		
+	retargetArctoArc(iarc);
+	
+	enode = OTHER_NODE(earc, enode);
+	inode = RIG_otherNode(iarc, inode);
+	
+	inode->link = enode;
+	
+	for(i = 0, next_iarc = inode->arcs[i]; next_iarc; i++, next_iarc = inode->arcs[i])
+	{
+		/* no back tracking */
+		if (next_iarc != iarc)
+		{
+			findCorrespondingArc(iarc, inode, next_iarc);
+			if (next_iarc->link)
+			{
+				retargetSubgraph(rigg, next_iarc, inode);
+			}
+		}
+	}
+}
+
+static void retargetGraphs(RigGraph *rigg)
+{
+	ReebGraph *reebg = rigg->link;
+	ReebArc *earc;
+	RigArc *iarc;
+	ReebNode *enode;
+	RigNode *inode;
+	
+	/* flag all ReebArcs as not taken */
+	for (earc = reebg->arcs.first; earc; earc = earc->next)
+	{
+		earc->flag = 0;
+	}
+	
+	earc = reebg->arcs.first;
+	iarc = rigg->head->arcs[0];
+	
+	iarc->link = earc;
+	earc->flag = 1;
+	
+	enode = earc->v1;
+	inode = iarc->tail;
+
+	inode->link = enode;
+
+	retargetSubgraph(rigg, iarc, inode);
+}
+
+void BIF_retargetArmature()
+{
+	Object *ob;
+	Base *base;
+	ReebGraph *reebg;
+	
+	reebg = BIF_ReebGraphFromEditMesh();
+	
+	markdownSymmetry(reebg);
+	
+	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);
+				
+				printf("Armature graph created\n");
+		
+				RIG_markdownSymmetry(rigg);
+				
+				RIG_printGraph(rigg);
+				
+				rigg->link = 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(rigg);
+			}
+		}
+	}
+
+	REEB_freeGraph(reebg);
+	
+	BIF_undo_push("Retarget Skeleton");
+	
+	exit_editmode(EM_FREEDATA|EM_FREEUNDO|EM_WAITCURSOR); // freedata, and undo
+
+	allqueue(REDRAWVIEW3D, 0);
+}
diff --git a/source/blender/src/buttons_editing.c b/source/blender/src/buttons_editing.c
index 1365baf075a..2b55fec084a 100644
--- a/source/blender/src/buttons_editing.c
+++ b/source/blender/src/buttons_editing.c
@@ -4891,6 +4891,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! */
@@ -4989,6 +4992,38 @@ static void skgen_reorder(void *option, void *arg2)
 	}
 }
 
+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", "Editing");
+	if(uiNewPanel(curarea, block, "Skeleton Retargetting", "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");
+
+	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, 							"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");
+	uiBlockEndAlign(block);
+
+	uiDefButF(block, NUM, B_DIFF, 							"Angle:",			1025, 60, 125,19, &G.scene->toolsettings->skgen_retarget_angle_weight, 0, 10, 1, 0,		"Angle Weight");
+	uiDefButF(block, NUM, B_DIFF, 							"Length:",			1150, 60, 125,19, &G.scene->toolsettings->skgen_retarget_length_weight, 0, 10, 1, 0,		"Length 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;
@@ -4996,12 +5031,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);
+	uiNewPanelTabbed("Mesh Tools More", "Editing");
 	if(uiNewPanel(curarea, block, "Skeleton 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");
 
 	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");
+	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, 							"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");
@@ -6511,8 +6548,8 @@ 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_uvautocalculation();
 			if (EM_texFaceCheck())
diff --git a/source/blender/src/editarmature.c b/source/blender/src/editarmature.c
index c166a9df762..37e6510786d 100644
--- a/source/blender/src/editarmature.c
+++ b/source/blender/src/editarmature.c
@@ -4142,542 +4142,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 ******************************************/
 
@@ -5002,8 +4467,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);
@@ -5030,35 +4493,35 @@ 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->v2->degree != 1)
 		{
 			head = arc->v2;
 			tail = arc->v1;
 			
-			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->v1->symmetry_level == 1)
 		{
 			head = arc->v1;
 			tail = arc->v2;
 			
-			arc->flags = 1; /* mark arc direction */
+			arc->flag = 1; /* mark arc direction */
 		}
-		else if (arc->v2->flags == 1)
+		else if (arc->v2->symmetry_level == 1)
 		{
 			head = arc->v2;
 			tail = arc->v1;
 			
-			arc->flags = -1; /* mark arc direction */
+			arc->flag = -1; /* mark arc direction */
 		}
 		/* otherwise, always go from low weight to high weight */
 		else
@@ -5066,7 +4529,7 @@ void generateSkeletonFromReebGraph(ReebGraph *rg)
 			head = arc->v1;
 			tail = arc->v2;
 			
-			arc->flags = 1; /* mark arc direction */
+			arc->flag = 1; /* mark arc direction */
 		}
 		
 		/* Loop over subdivision methods */	
@@ -5113,7 +4576,7 @@ void generateSkeletonFromReebGraph(ReebGraph *rg)
 			arc = 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->v1 == node && arc->flag == -1) || (arc->v2 == node && arc->flag == 1))
 			{
 				if (incomingArc == NULL)
 				{
@@ -5139,7 +4602,7 @@ void generateSkeletonFromReebGraph(ReebGraph *rg)
 				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->v1 == node && arc->flag == 1) || (arc->v2 == node && arc->flag == -1))
 				{
 					EditBone *childBone = BLI_ghash_lookup(arcBoneMap, arc);
 
@@ -5157,89 +4620,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);
-		}
-	}
-
-	verifyBuckets(rg);
+	reebg = BIF_ReebGraphFromEditMesh();
 
-	repositionNodes(rg);
-	
-	verifyBuckets(rg);
+	generateSkeletonFromReebGraph(reebg);
 
-	/* Filtering might have created degree 2 nodes, so remove them */
-	removeNormalNodes(rg);
-	
-	verifyBuckets(rg);
+	REEB_freeGraph(reebg);
 
-	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);
 }
diff --git a/source/blender/src/reeb.c b/source/blender/src/reeb.c
index 85fb5815c3e..34230e6dfc2 100644
--- a/source/blender/src/reeb.c
+++ b/source/blender/src/reeb.c
@@ -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"
 
@@ -41,6 +42,7 @@
 #include "BLI_arithb.h"
 #include "BLI_editVert.h"
 #include "BLI_edgehash.h"
+#include "BLI_ghash.h"
 
 #include "BDR_editobject.h"
 
@@ -60,6 +62,9 @@
 
 #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))
+
 /*
  * Skeleton generation algorithm based on: 
  * "Harmonic Skeleton for Realistic Character Animation"
@@ -72,10 +77,20 @@
  * SIGGRAPH 2007
  * 
  * */
+ 
+#define DEBUG_REEB
+
+typedef enum {
+	MERGE_LOWER,
+	MERGE_HIGHER,
+	MERGE_APPEND
+} MergeDirection;
 
 int mergeArcs(ReebGraph *rg, ReebArc *a0, ReebArc *a1);
 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 **********************************************/
 
@@ -227,11 +242,14 @@ void freeArc(ReebArc *arc)
 	
 	if (arc->buckets)
 		MEM_freeN(arc->buckets);
+		
+	if (arc->faces)
+		BLI_ghash_free(arc->faces, NULL, NULL);
 	
 	MEM_freeN(arc);
 }
 
-void freeGraph(ReebGraph *rg)
+void REEB_freeGraph(ReebGraph *rg)
 {
 	ReebArc *arc;
 	ReebNode *node;
@@ -292,6 +310,7 @@ void repositionNodes(ReebGraph *rg)
 
 void verifyNodeDegree(ReebGraph *rg)
 {
+#ifdef DEBUG_REEB
 	ReebNode *node = NULL;
 	ReebArc *arc = NULL;
 
@@ -310,6 +329,7 @@ void verifyNodeDegree(ReebGraph *rg)
 			printf("degree error in node %i: expected %i got %i\n", node->index, count, node->degree);
 		}
 	}
+#endif
 }
 
 void verifyBuckets(ReebGraph *rg)
@@ -345,9 +365,617 @@ void verifyBuckets(ReebGraph *rg)
 #endif
 }
 
+void verifyFaces(ReebGraph *rg)
+{
+#ifdef DEBUG_REEB
+	int total = 0;
+	ReebArc *arc = NULL;
+	for(arc = rg->arcs.first; arc; arc = arc->next)
+	{
+		total += BLI_ghash_size(arc->faces);
+	}
+	
+#endif
+}
+
+/**************************************** SYMMETRY HANDLING ******************************************/
+
+void markdownSymmetryArc(ReebArc *arc, ReebNode *node, int level);
+
+static 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], int reestablish)
+{
+	RadialArc *ring = NULL;
+	RadialArc *unit;
+	float limit = G.scene->toolsettings->skgen_symmetry_limit;
+	int symmetric = 1;
+	int count = 0;
+	int i;
+
+	/* mark topological symmetry */
+	node->symmetry_flag |= SYM_TOPOLOGICAL;
+
+	/* 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->symmetry_level == -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->symmetry_level == -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)
+	{
+		/* mark node as symmetric physically */
+		VECCOPY(node->symmetry_axis, axis);
+		node->symmetry_flag |= SYM_PHYSICAL;
+		node->symmetry_flag |= SYM_RADIAL;
+		
+		/* reestablish symmetry only if wanted */
+		if (reestablish)
+		{
+			/* 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);
+}
+
+static void setSideAxialSymmetry(ReebNode *root_node, ReebNode *end_node, ReebArc *arc)
+{
+	float vec[3];
+	
+	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;
+	}
+}
+
+void reestablishAxialSymmetry(ReebNode *node, int depth, float axis[3], int reestablish)
+{
+	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;
+	
+	/* mark topological symmetry */
+	node->symmetry_flag |= SYM_TOPOLOGICAL;
+
+	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->symmetry_level == -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(vec, node1->p, node->p);
+	Normalize(vec);
+	VecSubf(p, node->p, node2->p);
+	Normalize(p);
+	VecAddf(p, p, vec);
+	
+	
+	Crossf(vec, p, axis);
+	Crossf(nor, vec, axis);
+
+	printvecf("p", p);
+	printvecf("axis", axis);
+	printvecf("vec", vec);
+	printvecf("nor", nor);
+	
+	/* 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)
+	{
+		/* mark node as symmetric physically */
+		VECCOPY(node->symmetry_axis, nor);
+		node->symmetry_flag |= SYM_PHYSICAL;
+		node->symmetry_flag |= SYM_AXIAL;
+		
+		/* set side on arcs */
+		setSideAxialSymmetry(node, node1, arc1);
+		setSideAxialSymmetry(node, node2, arc2);
+		
+		/* reestablish symmetry only if wanted */
+		if (reestablish)
+		{
+			/* 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);
+				}
+			}
+		}
+	}
+	else
+	{
+		printf("NOT SYMMETRIC!\n");
+		printf("%f <= %f\n", VecLenf(node1->p, p), limit);
+		printvecf("axis", 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 */
+	int reestablish = 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->symmetry_level == -depth)
+		{
+			count++;
+		}
+		/* If arc is on the axis */
+		else if (connectedArc->symmetry_level == 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, reestablish);
+	}
+	else
+	{
+		reestablishRadialSymmetry(node, depth, axis, reestablish);
+	}
+
+	/* markdown secondary symetries */	
+	for (i = 0; node->arcs[i] != NULL; i++)
+	{
+		ReebArc *connectedArc = node->arcs[i];
+		
+		if (connectedArc->symmetry_level == -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->symmetry_level = 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->symmetry_level = -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->symmetry_level < 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->symmetry_level == connectedArc->symmetry_level)
+				{
+					/* 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->symmetry_level < 0)
+		{
+			/* subtree depth is store as a negative value in the flag */
+			markdownSecondarySymmetry(node, -connectedArc->symmetry_level, 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->symmetry_level = 0;
+	}
+	
+	/* mark down all nodes as not on the symmetry axis */
+	for (node = rg->nodes.first; node; node = node->next)
+	{
+		node->symmetry_level = 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->symmetry_level < 0)
+			{
+				arc->symmetry_level = 0;
+			}
+			else
+			{
+				/* mark down nodes with the lowest level symmetry axis */
+				if (arc->v1->symmetry_level == 0 || arc->v1->symmetry_level > arc->symmetry_level)
+				{
+					arc->v1->symmetry_level = arc->symmetry_level;
+				}
+				if (arc->v2->symmetry_level == 0 || arc->v2->symmetry_level > arc->symmetry_level)
+				{
+					arc->v2->symmetry_level = arc->symmetry_level;
+				}
+			}
+		}
+	}
+}
+
 /************************************** ADJACENCY LIST *************************************************/
 
-void addArcToNodeAdjacencyList(ReebNode *node, ReebArc *arc)
+static void addArcToNodeAdjacencyList(ReebNode *node, ReebArc *arc)
 {
 	ReebArc **arclist;
 
@@ -590,6 +1218,7 @@ void filterArc(ReebGraph *rg, ReebNode *newNode, ReebNode *removedNode, ReebArc
 			else
 			{
 				newNode->degree++; // incrementing degree since we're adding an arc
+				mergeArcFaces(rg, arc, srcArc);
 
 				if (merging)
 				{
@@ -762,6 +1391,181 @@ int filterExternalReebGraph(ReebGraph *rg, float threshold)
 	return value;
 }
 
+int filterSmartReebGraph(ReebGraph *rg, float threshold)
+{
+	ReebArc *arc = NULL, *nextArc = NULL;
+	int value = 0;
+	
+	BLI_sortlist(&rg->arcs, compareArcs);
+
+#ifdef DEBUG_REEB
+	{	
+		EditFace *efa;
+		for(efa=G.editMesh->faces.first; efa; efa=efa->next) {
+			efa->tmp.fp = -1;
+		}
+	}
+#endif
+
+	arc = rg->arcs.first;
+	while(arc)
+	{
+		nextArc = arc->next;
+		
+		/* need correct normals and center */
+		recalc_editnormals();
+
+		// Only test terminal arcs
+		if (arc->v1->degree == 1 || arc->v2->degree == 1)
+		{
+			GHashIterator ghi;
+			int merging = 0;
+			int total = BLI_ghash_size(arc->faces);
+			float avg_angle = 0;
+			float avg_vec[3] = {0,0,0};
+			
+			for(BLI_ghashIterator_init(&ghi, arc->faces);
+				!BLI_ghashIterator_isDone(&ghi);
+				BLI_ghashIterator_step(&ghi))
+			{
+				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->v1);
+		
+				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->v1->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
+			}
+
+
+#if 0			
+			avg_angle /= total;
+#else
+			VecMulf(avg_vec, 1.0 / total);
+			avg_angle = Inpf(avg_vec, avg_vec);
+#endif
+			
+			arc->angle = avg_angle;
+			
+#ifdef DEBUG_REEB
+			printf("angle %f total %i\n", avg_angle, total);
+#endif
+			
+			if (avg_angle > threshold)
+				merging = 1;
+			
+			if (merging)
+			{
+				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)
+				{
+					terminalNode = arc->v1;
+					middleNode = arc->v2;
+				}
+				else
+				{
+					terminalNode = arc->v2;
+					middleNode = arc->v1;
+				}
+				
+				// 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--;
+				}
+	
+				// Reset nextArc, it might have changed
+				nextArc = arc->next;
+				
+				BLI_remlink(&rg->arcs, arc);
+				freeArc(arc);
+				
+				BLI_freelinkN(&rg->nodes, removedNode);
+				value = 1;
+			}
+		}
+		
+		arc = nextArc;
+	}
+	
+	return value;
+}
+
 /************************************** WEIGHT SPREADING ***********************************************/
 
 int compareVerts( const void* a, const void* b )
@@ -858,12 +1662,12 @@ int detectCycle(ReebNode *node, ReebArc *srcArc)
 {
 	int value = 0;
 	
-	if (node->flags == 0)
+	if (node->flag == 0)
 	{
 		ReebArc ** pArc;
 
 		/* mark node as visited */
-		node->flags = 1;
+		node->flag = 1;
 
 		for(pArc = node->arcs; *pArc && value == 0; pArc++)
 		{
@@ -894,14 +1698,14 @@ int	isGraphCyclic(ReebGraph *rg)
 	/* Mark all nodes as not visited */
 	for(node = rg->nodes.first; node; node = node->next)
 	{
-		node->flags = 0;
+		node->flag = 0;
 	}
 
 	/* detectCycles in subgraphs */	
 	for(node = rg->nodes.first; node && value == 0; node = node->next)
 	{
 		/* only for nodes in subgraphs that haven't been visited yet */
-		if (node->flags == 0)
+		if (node->flag == 0)
 		{
 			value = value || detectCycle(node, NULL);
 		}		
@@ -917,9 +1721,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,6 +1740,7 @@ 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);
 		
@@ -945,11 +1749,14 @@ void exportGraph(ReebGraph *rg, int count)
 			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]);
 		}
 		
+		VecAddf(p, arc->v2->p, arc->v1->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->v2);
 	}	
 	
 	fclose(f);
-#endif
 }
 
 /***************************************** MAIN ALGORITHM **********************************************/
@@ -969,6 +1776,7 @@ ReebArc * findConnectedArc(ReebGraph *rg, ReebArc *arc, ReebNode *v)
 	return nextArc;
 }
 
+
 void removeNormalNodes(ReebGraph *rg)
 {
 	ReebArc *arc;
@@ -1041,11 +1849,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)
 {
@@ -1110,6 +1930,7 @@ int mergeConnectedArcs(ReebGraph *rg, ReebArc *a0, ReebArc *a1)
 	ReebNode *removedNode = NULL;
 	
 	mergeArcEdges(rg, a0, a1, MERGE_APPEND);
+	mergeArcFaces(rg, a0, a1);
 	
 	// Bring a0 to the combine length of both arcs
 	if (a0->v2 == a1->v1)
@@ -1146,6 +1967,7 @@ int mergeArcs(ReebGraph *rg, ReebArc *a0, ReebArc *a1)
 		if (a0->v2->weight == a1->v2->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);
 			
@@ -1162,6 +1984,7 @@ int mergeArcs(ReebGraph *rg, ReebArc *a0, ReebArc *a1)
 		else if (a0->v2->weight > a1->v2->weight) // a1->v2->weight is in the middle
 		{
 			mergeArcEdges(rg, a1, a0, MERGE_LOWER);
+			mergeArcFaces(rg, a1, a0);
 
 			// Adjust node degree
 			a0->v1->degree--;
@@ -1174,6 +1997,7 @@ int mergeArcs(ReebGraph *rg, ReebArc *a0, ReebArc *a1)
 		else // a0>n2 is in the middle
 		{
 			mergeArcEdges(rg, a0, a1, MERGE_LOWER);
+			mergeArcFaces(rg, a0, a1);
 			
 			// Adjust node degree
 			a1->v1->degree--;
@@ -1190,6 +2014,7 @@ int mergeArcs(ReebGraph *rg, ReebArc *a0, ReebArc *a1)
 		if (a0->v1->weight > a1->v1->weight) // a0->v1->weight is in the middle
 		{
 			mergeArcEdges(rg, a0, a1, MERGE_HIGHER);
+			mergeArcFaces(rg, a0, a1);
 			
 			// Adjust node degree
 			a1->v2->degree--;
@@ -1202,6 +2027,7 @@ int mergeArcs(ReebGraph *rg, ReebArc *a0, ReebArc *a1)
 		else // a1->v1->weight is in the middle
 		{
 			mergeArcEdges(rg, a1, a0, MERGE_HIGHER);
+			mergeArcFaces(rg, a1, a0);
 
 			// Adjust node degree
 			a0->v2->degree--;
@@ -1258,7 +2084,8 @@ ReebNode * addNode(ReebGraph *rg, EditVert *eve, float weight)
 	
 	node = MEM_callocN(sizeof(ReebNode), "reeb node");
 	
-	node->flags = 0; // clear flags on init
+	node->flag = 0; // clear flag on init
+	node->symmetry_level = 0;
 	node->arcs = NULL;
 	node->degree = 0;
 	node->weight = weight;
@@ -1288,7 +2115,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)
 		{
@@ -1349,7 +2178,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 +2188,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);
@@ -1398,6 +2231,17 @@ void addTriangleToGraph(ReebGraph *rg, ReebNode * n1, ReebNode * n2, ReebNode *
 	mergePaths(rg, e1, e2, e3);
 }
 
+ReebGraph * newReebGraph()
+{
+	ReebGraph *rg;
+	rg = MEM_callocN(sizeof(ReebGraph), "reeb graph");
+	
+	rg->totnodes = 0;
+	rg->emap = BLI_edgehash_new();
+	
+	return rg;
+}
+
 ReebGraph * generateReebGraph(EditMesh *em, int subdivisions)
 {
 	ReebGraph *rg;
@@ -1412,10 +2256,7 @@ ReebGraph * generateReebGraph(EditMesh *em, int subdivisions)
 	int countfaces = 0;
 #endif
  	
-	rg = MEM_callocN(sizeof(ReebGraph), "reeb graph");
-	
-	rg->totnodes = 0;
-	rg->emap = BLI_edgehash_new();
+	rg = newReebGraph();
 	
 	totvert = BLI_countlist(&em->verts);
 	totfaces = BLI_countlist(&em->faces);
@@ -1447,12 +2288,12 @@ ReebGraph * generateReebGraph(EditMesh *em, int subdivisions)
 		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);
+		addTriangleToGraph(rg, n1, n2, n3, efa);
 		
 		if (efa->v4)
 		{
 			ReebNode *n4 = (ReebNode*)efa->v4->tmp.p;
-			addTriangleToGraph(rg, n1, n3, n4);
+			addTriangleToGraph(rg, n1, n3, n4, efa);
 		}
 
 #ifdef DEBUG_REEB
@@ -1460,6 +2301,7 @@ ReebGraph * generateReebGraph(EditMesh *em, int subdivisions)
 		if (countfaces % 100 == 0)
 		{
 			printf("face %i of %i\n", countfaces, totfaces);
+			verifyFaces(rg);
 		}
 #endif
 		
@@ -1728,7 +2570,7 @@ int weightFromDistance(EditMesh *em)
 		return 0;
 	}
 
-	/* Initialize vertice flags and find at least one selected vertex */
+	/* Initialize vertice flag and find at least one selected vertex */
 	for(eve = em->verts.first; eve && vCount == 0; eve = eve->next)
 	{
 		eve->f1 = 0;
@@ -1762,7 +2604,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);
@@ -1837,17 +2679,17 @@ int weightFromDistance(EditMesh *em)
 	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;
@@ -1856,14 +2698,148 @@ void weightToVCol(EditMesh *em)
 	}
 	
 	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);
 				
-		mcol[0] = MColFromWeight(efa->v1);
-		mcol[1] = MColFromWeight(efa->v2);
-		mcol[2] = MColFromWeight(efa->v3);
+		mcol[0].a = 255;
+		mcol[1].a = 255;
+		mcol[2].a = 255;
 
 		if(efa->v4) {
-			mcol[3] = MColFromWeight(efa->v4);
+			mcol[3].a = 255;
 		}
 	}
 }
@@ -1890,6 +2866,31 @@ void initArcIterator(ReebArcIterator *iter, ReebArc *arc, ReebNode *head)
 	iter->index = iter->start - iter->stride;
 }
 
+void initArcIteratorStart(struct ReebArcIterator *iter, struct ReebArc *arc, struct ReebNode *head, int start)
+{
+	iter->arc = arc;
+	
+	if (head == arc->v1)
+	{
+		iter->start = start;
+		iter->end = arc->bcount - 1;
+		iter->stride = 1;
+	}
+	else
+	{
+		iter->start = arc->bcount - 1 - start;
+		iter->end = 0;
+		iter->stride = -1;
+	}
+	
+	iter->index = iter->start - iter->stride;
+	
+	if (start >= arc->bcount)
+	{
+		iter->start = iter->end; /* stop iterator since it's past its end */
+	}
+}
+
 void initArcIterator2(ReebArcIterator *iter, ReebArc *arc, int start, int end)
 {
 	iter->arc = arc;
@@ -1921,3 +2922,153 @@ EmbedBucket * nextBucket(ReebArcIterator *iter)
 	
 	return result;
 }
+
+EmbedBucket * nextNBucket(ReebArcIterator *iter, int n)
+{
+	EmbedBucket *result = NULL;
+	
+	iter->index += n * iter->stride;
+
+	/* check if passed end */
+	if ((iter->stride == 1 && iter->index < iter->end) ||
+		(iter->stride == -1 && iter->index > iter->end))
+	{
+		result = &(iter->arc->buckets[iter->index]);
+	}
+	else
+	{
+		/* stop iterator if passed end */
+		iter->index = iter->end; 
+	}
+	
+	return result;
+}
+
+EmbedBucket * previousBucket(struct ReebArcIterator *iter)
+{
+	EmbedBucket *result = NULL;
+	
+	if (iter->index != iter->start)
+	{
+		iter->index -= iter->stride;
+		result = &(iter->arc->buckets[iter->index]);
+	}
+	
+	return result;
+}
+
+int iteratorStopped(struct ReebArcIterator *iter)
+{
+	if (iter->index == iter->end)
+	{
+		return 1;
+	}
+	else
+	{
+		return 0;
+	}
+}
+
+struct EmbedBucket * currentBucket(struct ReebArcIterator *iter)
+{
+	EmbedBucket *result = NULL;
+	
+	if (iter->index != iter->end)
+	{
+		result = &(iter->arc->buckets[iter->index]);
+	}
+	
+	return result;
+}
+
+/************************ PUBLIC FUNCTIONS *********************************************/
+
+ReebGraph *BIF_ReebGraphFromEditMesh(void)
+{
+	EditMesh *em = G.editMesh;
+	ReebGraph *rg = NULL;
+	int i;
+	
+	if (em == NULL)
+		return NULL;
+
+	if (weightFromDistance(em) == 0)
+	{
+		error("No selected vertex\n");
+		return NULL;
+	}
+	
+	renormalizeWeight(em, 1.0f);
+
+	if (G.scene->toolsettings->skgen_options & SKGEN_HARMONIC)
+	{
+		weightToHarmonic(em);
+	}
+	
+#ifdef DEBUG_REEB
+	weightToVCol(em, 1);
+#endif
+	
+	rg = generateReebGraph(em, G.scene->toolsettings->skgen_resolution);
+
+	verifyBuckets(rg);
+	
+	verifyFaces(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);
+		}
+	}
+
+	filterSmartReebGraph(rg, 0.5);
+
+#ifdef DEBUG_REEB
+	arcToVCol(rg, em, 0);
+	//angleToVCol(em, 1);
+#endif
+
+	verifyBuckets(rg);
+
+	repositionNodes(rg);
+	
+	verifyBuckets(rg);
+
+	/* 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);
+	
+	REEB_exportGraph(rg, -1);
+	
+	return rg;
+}