Merge from Harmonic Skeleton branch
authorMartin Poirier <theeth@yahoo.com>
Mon, 10 Dec 2007 21:14:19 +0000 (21:14 +0000)
committerMartin Poirier <theeth@yahoo.com>
Mon, 10 Dec 2007 21:14:19 +0000 (21:14 +0000)
This code adds a basic and simple skeleton generator.

Examples and links are in the wiki, docs will come eventually: http://wiki.blender.org/index.php/User:Theeth/skeletor

In a nutshell, select a vertex at the top of the head and press "Generate Skeleton". UI Panel is in the Editing buttons in Edit Mode, tooltips and semi-useful.

projectfiles_vc7/blender/blender.sln
projectfiles_vc7/blender/src/BL_src.vcproj
source/blender/blenlib/BLI_arithb.h
source/blender/blenloader/intern/readfile.c
source/blender/include/BIF_editarmature.h
source/blender/include/butspace.h
source/blender/include/reeb.h [new file with mode: 0644]
source/blender/makesdna/DNA_scene_types.h
source/blender/src/buttons_editing.c
source/blender/src/editarmature.c
source/blender/src/reeb.c [new file with mode: 0644]

index 0888cf7971a7d57487ff0dd0b5a84a787ee62707..60c2d8f2d4c821e30cbdfe8399a2e19243d84316 100644 (file)
@@ -248,8 +248,6 @@ Global
                Debug = Debug
                Release = Release
        EndGlobalSection
-       GlobalSection(ProjectDependencies) = postSolution
-       EndGlobalSection
        GlobalSection(ProjectConfiguration) = postSolution
                {F78B7FC9-DE32-465E-9F26-BB0B6B7A2EAF}.3D Plugin Debug.ActiveCfg = Blender Debug|Win32
                {F78B7FC9-DE32-465E-9F26-BB0B6B7A2EAF}.3D Plugin Release.ActiveCfg = Blender Release|Win32
index 6cc8f27f3988bdd350b1f32b5a649c3519bdcab7..9c03b21f1a82c689a6fbd595368fe9621f1ee821 100644 (file)
                        <File
                                RelativePath="..\..\..\source\blender\src\prvicons.c">
                        </File>
+                       <File
+                               RelativePath="..\..\..\source\blender\src\reeb.c">
+                       </File>
                        <File
                                RelativePath="..\..\..\source\blender\src\renderwin.c">
                        </File>
                        <File
                                RelativePath="..\..\..\source\blender\include\playanim_ext.h">
                        </File>
+                       <File
+                               RelativePath="..\..\..\source\blender\include\reeb.h">
+                       </File>
                        <File
                                RelativePath="..\..\..\source\blender\include\transform.h">
                        </File>
index 11a1fd277da2065fff5a20f3c111311d0a15c4db..1e86f8fe5b37b6150b3ac985431eb17052bdb54d 100644 (file)
@@ -54,6 +54,10 @@ extern "C" {
 #define M_SQRT1_2      0.70710678118654752440
 #endif
 
+#ifdef WIN32
+    #define isnan(n) _isnan(n)
+#endif
+
 #define MAT4_UNITY {{ 1.0, 0.0, 0.0, 0.0},\
                                        { 0.0, 1.0, 0.0, 0.0},\
                                        { 0.0, 0.0, 1.0, 0.0},\
index 68bf1d1815c75b57942c89dfdc16e91c556576c6..4b8418a78b908fd9a727d31fc9f6fcf2e356e15d 100644 (file)
@@ -6652,7 +6652,7 @@ static void do_versions(FileData *fd, Library *lib, Main *main)
                                else
                                        wrld->ao_samp_method = WO_AOSAMP_HAMMERSLEY;
                                
-                               wrld->ao_adapt_thresh = 0.005;
+                               wrld->ao_adapt_thresh = 0.005f;
                        }
                        
                        for(la=main->lamp.first; la; la= la->id.next) {
@@ -6666,6 +6666,7 @@ static void do_versions(FileData *fd, Library *lib, Main *main)
                }
        }
        if(main->versionfile <= 245) {
+               Scene *sce;
                bScreen *sc;
                Object *ob;
                Image *ima;
@@ -6838,6 +6839,27 @@ static void do_versions(FileData *fd, Library *lib, Main *main)
                for(part=main->particle.first; part; part=part->id.next)
                        if(part->ren_child_nbr==0)
                                part->ren_child_nbr= part->child_nbr;
+               if (main->versionfile < 245 || main->subversionfile < 11)
+               {
+                       /* initialize skeleton generation toolsettings */
+                       for(sce=main->scene.first; sce; sce = sce->id.next)
+                       {
+                               sce->toolsettings->skgen_resolution = 50;
+                               sce->toolsettings->skgen_threshold_internal     = 0.01f;
+                               sce->toolsettings->skgen_threshold_external     = 0.01f;
+                               sce->toolsettings->skgen_angle_limit                    = 45.0f;
+                               sce->toolsettings->skgen_length_ratio                   = 1.3f;
+                               sce->toolsettings->skgen_length_limit                   = 1.5f;
+                               sce->toolsettings->skgen_correlation_limit              = 0.98f;
+                               sce->toolsettings->skgen_symmetry_limit                 = 0.1f;
+                               sce->toolsettings->skgen_postpro = SKGEN_SMOOTH;
+                               sce->toolsettings->skgen_postpro_passes = 1;
+                               sce->toolsettings->skgen_options = SKGEN_FILTER_INTERNAL|SKGEN_FILTER_EXTERNAL|SKGEN_SUB_CORRELATION;
+                               sce->toolsettings->skgen_subdivisions[0] = SKGEN_SUB_CORRELATION;
+                               sce->toolsettings->skgen_subdivisions[1] = SKGEN_SUB_LENGTH;
+                               sce->toolsettings->skgen_subdivisions[2] = SKGEN_SUB_ANGLE;
+                       }
+               }
        }
 
        if ((main->versionfile < 245) || (main->versionfile == 245 && main->subversionfile < 2)) {
index 903c663e5f3783f040c7fb3fca79539478f4d3c1..31a5cc2ec4c4612c283f9871627e8568bf777331 100644 (file)
@@ -101,6 +101,8 @@ void        make_trans_bones (char mode);
 
 int            do_pose_selectbuffer(struct Base *base, unsigned int *buffer, short hits);
 
+void generateSkeleton(void);
+
 void   mouse_armature(void);
 void   remake_editArmature(void);
 void   selectconnected_armature(void);
@@ -117,8 +119,8 @@ void        hide_selected_pose_bones(void);
 void   hide_unselected_pose_bones(void);
 void   show_all_pose_bones(void);
 
-int            bone_looper(Object *ob, struct Bone *bone, void *data,
-                               int (*bone_func)(Object *, struct Bone *, void *));
+int            bone_looper(struct Object *ob, struct Bone *bone, void *data,
+                               int (*bone_func)(struct Object *, struct Bone *, void *));
 
 void   undo_push_armature(char *name);
 void   armature_bone_rename(struct bArmature *arm, char *oldname, char *newname);
index debac14f63e996505342c3af4d336e7835092310..5029c88fba753162a167a6f85df750d678e9b69a 100644 (file)
@@ -420,6 +420,8 @@ void curvemap_buttons(struct uiBlock *block, struct CurveMapping *cumap, char la
 #define B_SETTFACE_RND         2082
 #define B_SETMCOL_RND          2083
 
+#define B_GEN_SKELETON         2090
+
 /* *********************** */
 #define B_VGROUPBUTS           2100
 
diff --git a/source/blender/include/reeb.h b/source/blender/include/reeb.h
new file mode 100644 (file)
index 0000000..c8352ae
--- /dev/null
@@ -0,0 +1,127 @@
+/**
+ * $Id: 
+ *
+ * ***** BEGIN GPL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version. The Blender
+ * Foundation also sells licenses for use in proprietary software under
+ * the Blender License.  See http://www.blender.org/BL/ for information
+ * about this.
+ *
+ * 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 *****
+ */
+#ifndef REEB_H_
+#define REEB_H_
+
+#include "DNA_listBase.h"
+
+struct EdgeHash;
+struct ReebArc;
+struct ReebEdge;
+struct ReebNode;
+
+typedef struct ReebGraph {
+       ListBase arcs;
+       ListBase nodes;
+       int totnodes;
+       struct EdgeHash *emap;
+} ReebGraph;
+
+typedef struct EmbedBucket {
+       float val;
+       int       nv;
+       float p[3];
+} EmbedBucket;
+
+typedef struct ReebNode {
+       struct ReebNode *next, *prev;
+       struct ReebArc **arcs;
+       int index;
+       int degree;
+       float weight;
+       float p[3];
+       int flags;
+} ReebNode;
+
+typedef struct ReebEdge {
+       struct ReebEdge *next, *prev;
+       struct ReebArc  *arc;
+       struct ReebNode *v1, *v2;
+       struct ReebEdge *nextEdge;
+} ReebEdge;
+
+typedef struct ReebArc {
+       struct ReebArc *next, *prev;
+       ListBase edges;
+       struct ReebNode *v1, *v2;
+       struct EmbedBucket *buckets;
+       int     bcount;
+       int flags;
+} ReebArc;
+
+typedef struct ReebArcIterator {
+       struct ReebArc  *arc;
+       int index;
+       int start;
+       int end;
+       int stride;
+} ReebArcIterator;
+
+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 renormalizeWeight(struct EditMesh *em, float newmax);
+
+ReebGraph * generateReebGraph(struct EditMesh *me, int subdivisions);
+void freeGraph(ReebGraph *rg);
+void exportGraph(ReebGraph *rg, int count);
+
+#define OTHER_NODE(arc, node) ((arc->v1 == node) ? arc->v2 : arc->v1)
+
+void initArcIterator(struct ReebArcIterator *iter, struct ReebArc *arc, struct ReebNode *head);
+void initArcIterator2(struct ReebArcIterator *iter, struct ReebArc *arc, int start, int end);
+struct EmbedBucket * nextBucket(struct ReebArcIterator *iter);
+
+/* Filtering */
+void filterNullReebGraph(ReebGraph *rg);
+int filterExternalReebGraph(ReebGraph *rg, float threshold);
+int filterInternalReebGraph(ReebGraph *rg, float threshold);
+
+/* Post-Build processing */
+void repositionNodes(ReebGraph *rg);
+void postprocessGraph(ReebGraph *rg, char mode);
+void removeNormalNodes(ReebGraph *rg);
+
+/* Graph processing */
+void buildAdjacencyList(ReebGraph *rg);
+
+void sortNodes(ReebGraph *rg);
+void sortArcs(ReebGraph *rg);
+
+int subtreeDepth(ReebNode *node, ReebArc *rootArc);
+int countConnectedArcs(ReebGraph *rg, ReebNode *node);
+int hasAdjacencyList(ReebGraph *rg); 
+int    isGraphCyclic(ReebGraph *rg);
+
+/* Sanity check */
+void verifyBuckets(ReebGraph *rg);
+
+#endif /*REEB_H_*/
index f9418d1fc85fa2eff139083c779a727f448ce7df..2dfe0ac43b6843a01637209a5c120e5e730b9a01 100644 (file)
@@ -364,7 +364,7 @@ typedef struct ToolSettings {
        short uvcalc_flag;
 
        short pad2;
-       
+
        /* Image Paint (8 byte aligned please!) */
        struct ImagePaintSettings imapaint;
 
@@ -385,7 +385,21 @@ typedef struct ToolSettings {
        /* Multires */
        char multires_subdiv_type;
        
-       char  pad4[2];
+       /* Skeleton generation */
+       short skgen_resolution;
+       float skgen_threshold_internal;
+       float skgen_threshold_external;
+       float skgen_length_ratio;
+       float skgen_length_limit;
+       float skgen_angle_limit;
+       float skgen_correlation_limit;
+       float skgen_symmetry_limit;
+       short skgen_options;
+       char  skgen_postpro;
+       char  skgen_postpro_passes;
+       char  skgen_subdivisions[3];
+       
+       char pad3[5];
 } ToolSettings;
 
 /* Used by all brushes to store their properties, which can be directly set
@@ -723,6 +737,24 @@ typedef struct Scene {
 #define RETOPO_LINE 2
 #define RETOPO_ELLIPSE 4
 
+/* toolsettings->skgen_options */
+#define SKGEN_FILTER_INTERNAL  1
+#define SKGEN_FILTER_EXTERNAL  2
+#define        SKGEN_SYMMETRY                  4
+#define        SKGEN_CUT_LENGTH                8
+#define        SKGEN_CUT_ANGLE                 16
+#define        SKGEN_CUT_CORRELATION   32
+
+#define        SKGEN_SUB_LENGTH                0
+#define        SKGEN_SUB_ANGLE                 1
+#define        SKGEN_SUB_CORRELATION   2
+#define        SKGEN_SUB_TOTAL                 3
+
+/* toolsettings->skgen_postpro */
+#define SKGEN_SMOOTH                   0
+#define SKGEN_AVERAGE                  1
+#define SKGEN_SHARPEN                  2
+
 #ifdef __cplusplus
 }
 #endif
index ebe7986cd75492018cbb226acd363092713abfa9..a6ba07bfb34ee0e8f78abb564b0cfe726d202b8b 100644 (file)
@@ -4621,6 +4621,9 @@ void do_meshbuts(unsigned short event)
        case B_JOINTRIA:
                join_triangles();
                break;
+       case B_GEN_SKELETON:
+               generateSkeleton();
+               break;
        }
 
        /* WATCH IT: previous events only in editmode! */
@@ -4671,6 +4674,7 @@ static void editing_panel_mesh_tools(Object *ob, Mesh *me)
        uiDefButS(block, NUM, B_DIFF, "Turns:",                 210,55,115,19, &G.scene->toolsettings->turn,1.0,360.0, 0, 0, "Specifies the number of revolutions the screw turns");
        uiDefButBitS(block, TOG, B_KEEPORIG, B_DIFF, "Keep Original",10,35,200,19, &G.scene->toolsettings->editbutflag, 0, 0, 0, 0, "Keeps a copy of the original vertices and faces after executing tools");
        uiDefButBitS(block, TOG, B_CLOCKWISE, B_DIFF, "Clockwise",      210,35,115,19, &G.scene->toolsettings->editbutflag, 0, 0, 0, 0, "Specifies the direction for 'Screw' and 'Spin'");
+       uiBlockEndAlign(block);
 
        uiBlockBeginAlign(block);
        uiDefBut(block, BUT,B_EXTREP, "Extrude Dup",    10,10,150,19, 0, 0, 0, 0, 0, "Creates copies of the selected vertices in a straight line away from the current viewport");
@@ -4694,7 +4698,81 @@ static void verify_vertexgroup_name_func(void *datav, void *data2_unused)
        unique_vertexgroup_name((bDeformGroup*)datav, OBACT);
 }
 
+static void skgen_reorder(void *option, void *arg2)
+{
+       char tmp;
+       switch ((int)option)
+       {
+               case 0:
+                       tmp = G.scene->toolsettings->skgen_subdivisions[0];
+                       G.scene->toolsettings->skgen_subdivisions[0] = G.scene->toolsettings->skgen_subdivisions[1];
+                       G.scene->toolsettings->skgen_subdivisions[1] = tmp;
+                       break;
+               case 1:
+                       tmp = G.scene->toolsettings->skgen_subdivisions[2];
+                       G.scene->toolsettings->skgen_subdivisions[2] = G.scene->toolsettings->skgen_subdivisions[1];
+                       G.scene->toolsettings->skgen_subdivisions[1] = tmp;
+                       break;
+               case 2:
+                       tmp = G.scene->toolsettings->skgen_subdivisions[0];
+                       G.scene->toolsettings->skgen_subdivisions[0] = G.scene->toolsettings->skgen_subdivisions[2];
+                       G.scene->toolsettings->skgen_subdivisions[2] = G.scene->toolsettings->skgen_subdivisions[1];
+                       G.scene->toolsettings->skgen_subdivisions[1] = tmp;
+                       break;
+       }
+}
+
+static void editing_panel_mesh_skgen(Object *ob, Mesh *me)
+{
+       uiBlock *block;
+       uiBut *but;
+       int i;
+
+       block= uiNewBlock(&curarea->uiblocks, "editing_panel_mesh_skgen", UI_EMBOSS, UI_HELV, curarea->win);
+       if(uiNewPanel(curarea, block, "Skeleton Generator", "Editing", 960, 0, 318, 204)==0) return;
+       
+       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");
+       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");
 
+       for(i = 0; i < SKGEN_SUB_TOTAL; i++)
+       {
+               int y = 90 - 20 * i;
+               
+               but = uiDefIconBut(block, BUT, B_MODIFIER_RECALC, VICON_MOVE_DOWN,              1025, y, 16, 19, NULL, 0.0, 0.0, 0.0, 0.0, "Change the order the subdivisions algorithm are applied");
+               uiButSetFunc(but, skgen_reorder, (void *)i, NULL);
+               
+               switch(G.scene->toolsettings->skgen_subdivisions[i])
+               {
+                       case SKGEN_SUB_LENGTH:
+                               uiDefButBitS(block, TOG, SKGEN_CUT_LENGTH, B_DIFF,              "Length",               1041, y, 67,19, &G.scene->toolsettings->skgen_options, 0, 0, 0, 0,                              "Subdivide arcs in bones of equal length");
+                               uiDefButF(block, NUM, B_DIFF,                                                   "T:",                   1111, y, 82,19, &G.scene->toolsettings->skgen_length_ratio,1.0, 4.0, 10, 0,             "Specify the ratio limit between straight arc and embeddings to trigger equal subdivisions");
+                               uiDefButF(block, NUM, B_DIFF,                                                   "L:",                   1193, y, 82,19, &G.scene->toolsettings->skgen_length_limit,0.1,50.0, 10, 0,             "Maximum length of the bones when subdividing");
+                               break;
+                       case SKGEN_SUB_ANGLE:
+                               uiDefButBitS(block, TOG, SKGEN_CUT_ANGLE, B_DIFF,               "Angle",                1041, y, 67,19, &G.scene->toolsettings->skgen_options, 0, 0, 0, 0,                                      "Subdivide arcs based on angle");
+                               uiDefButF(block, NUM, B_DIFF,                                                   "T:",                   1111, y,164,19, &G.scene->toolsettings->skgen_angle_limit,0.0, 90.0, 10, 0,                     "Specify the threshold angle in degrees for subdivision");
+                               break;
+                       case SKGEN_SUB_CORRELATION:
+                               uiDefButBitS(block, TOG, SKGEN_CUT_CORRELATION, B_DIFF, "Correlation",  1041, y, 67,19, &G.scene->toolsettings->skgen_options, 0, 0, 0, 0,                                      "Subdivide arcs based on correlation");
+                               uiDefButF(block, NUM, B_DIFF,                                                   "T:",                   1111, y,164,19, &G.scene->toolsettings->skgen_correlation_limit,0.0, 1.0, 0.01, 0,      "Specify the threshold correlation for subdivision");
+                               break;
+               }
+       }
+
+       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_tools1(Object *ob, Mesh *me)
 {
@@ -5951,6 +6029,8 @@ void editing_panels()
                if(G.obedit) {
                        editing_panel_mesh_tools(ob, ob->data);
                        editing_panel_mesh_tools1(ob, ob->data);
+                       uiNewPanelTabbed("Mesh Tools 1", "Editing");
+                       editing_panel_mesh_skgen(ob, ob->data);
                        editing_panel_mesh_uvautocalculation();
                        if (EM_texFaceCheck())
                                editing_panel_mesh_texface();
index 66ee7a83f3a83c4b804039cbfca713baf108b001..a77c965768dea88c946e8d9f2080a103ffcf08e9 100644 (file)
@@ -62,6 +62,7 @@
 #include "BLI_blenlib.h"
 #include "BLI_arithb.h"
 #include "BLI_editVert.h"
+#include "BLI_ghash.h"
 
 #include "BKE_action.h"
 #include "BKE_armature.h"
 
 #include "PIL_time.h"
 
+#include "reeb.h" // FIX ME
+
 #include "mydevice.h"
 #include "blendef.h"
 #include "nla.h"
@@ -332,7 +335,7 @@ void editbones_to_armature (ListBase *list, Object *ob)
        fix_bonelist_roll (&arm->bonebase, list);
        
        /* so all users of this armature should get rebuilt */
-       for(obt= G.main->object.first; obt; obt= obt->id.next) {
+       for (obt= G.main->object.first; obt; obt= obt->id.next) {
                if(obt->data==arm)
                        armature_rebuild_pose(obt, arm);
        }
@@ -3278,5 +3281,1108 @@ 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 = 2;
+               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 ******************************************/
+
+EditBone * subdivideByAngle(ReebArc *arc, ReebNode *head, ReebNode *tail)
+{
+       EditBone *lastBone = NULL;
+       if (G.scene->toolsettings->skgen_options & SKGEN_CUT_ANGLE)
+       {
+               ReebArcIterator iter;
+               EmbedBucket *current = NULL;
+               EmbedBucket *previous = NULL;
+               EditBone *child = NULL;
+               EditBone *parent = NULL;
+               EditBone *root = NULL;
+               float angleLimit = (float)cos(G.scene->toolsettings->skgen_angle_limit * M_PI / 180.0f);
+               
+               parent = add_editbone("Bone");
+               parent->flag |= BONE_SELECTED|BONE_TIPSEL|BONE_ROOTSEL;
+               VECCOPY(parent->head, head->p);
+               
+               root = parent;
+               
+               for (initArcIterator(&iter, arc, head), previous = nextBucket(&iter), current = nextBucket(&iter);
+                       current;
+                       previous = current, current = nextBucket(&iter))
+               {
+                       float vec1[3], vec2[3];
+                       float len1, len2;
+
+                       VecSubf(vec1, previous->p, parent->head);
+                       VecSubf(vec2, current->p, previous->p);
+
+                       len1 = Normalize(vec1);
+                       len2 = Normalize(vec2);
+
+                       if (len1 > 0.0f && len2 > 0.0f && Inpf(vec1, vec2) < angleLimit)
+                       {
+                               VECCOPY(parent->tail, previous->p);
+
+                               child = add_editbone("Bone");
+                               VECCOPY(child->head, parent->tail);
+                               child->parent = parent;
+                               child->flag |= BONE_CONNECTED|BONE_SELECTED|BONE_TIPSEL|BONE_ROOTSEL;
+                               
+                               parent = child; /* new child is next parent */
+                       }
+               }
+               VECCOPY(parent->tail, tail->p);
+               
+               /* If the bone wasn't subdivided, delete it and return NULL
+                * to let subsequent subdivision methods do their thing. 
+                * */
+               if (parent == root)
+               {
+                       delete_bone(parent);
+                       parent = NULL;
+               }
+               
+               lastBone = parent; /* set last bone in the chain */
+       }
+       
+       return lastBone;
+}
+
+float calcCorrelation(ReebArc *arc, int start, int end, float v0[3], float n[3])
+{
+       int len = 2 + abs(end - start);
+       
+       if (len > 2)
+       {
+               ReebArcIterator iter;
+               EmbedBucket *bucket = NULL;
+               float avg_t = 0.0f;
+               float s_t = 0.0f;
+               float s_xyz = 0.0f;
+               
+               /* First pass, calculate average */
+               for (initArcIterator2(&iter, arc, start, end), bucket = nextBucket(&iter);
+                       bucket;
+                       bucket = nextBucket(&iter))
+               {
+                       float v[3];
+                       
+                       VecSubf(v, bucket->p, v0);
+                       avg_t += Inpf(v, n);
+               }
+               
+               avg_t /= Inpf(n, n);
+               avg_t += 1.0f; /* adding start (0) and end (1) values */
+               avg_t /= len;
+               
+               /* Second pass, calculate s_xyz and s_t */
+               for (initArcIterator2(&iter, arc, start, end), bucket = nextBucket(&iter);
+                       bucket;
+                       bucket = nextBucket(&iter))
+               {
+                       float v[3], d[3];
+                       float dt;
+                       
+                       VecSubf(v, bucket->p, v0);
+                       Projf(d, v, n);
+                       VecSubf(v, v, d);
+                       
+                       dt = VecLength(d) - avg_t;
+                       
+                       s_t += dt * dt;
+                       s_xyz += Inpf(v, v);
+               }
+               
+               /* adding start(0) and end(1) values to s_t */
+               s_t += (avg_t * avg_t) + (1 - avg_t) * (1 - avg_t);
+               
+               return 1.0f - s_xyz / s_t; 
+       }
+       else
+       {
+               return 1.0f;
+       }
+}
+
+EditBone * subdivideByCorrelation(ReebArc *arc, ReebNode *head, ReebNode *tail)
+{
+       ReebArcIterator iter;
+       float n[3];
+       float CORRELATION_THRESHOLD = G.scene->toolsettings->skgen_correlation_limit;
+       EditBone *lastBone = NULL;
+       
+       /* init iterator to get start and end from head */
+       initArcIterator(&iter, arc, head);
+       
+       /* Calculate overall */
+       VecSubf(n, arc->buckets[iter.end].p, head->p);
+       
+       if (G.scene->toolsettings->skgen_options & SKGEN_CUT_CORRELATION && 
+               calcCorrelation(arc, iter.start, iter.end, head->p, n) < CORRELATION_THRESHOLD)
+       {
+               EmbedBucket *bucket = NULL;
+               EmbedBucket *previous = NULL;
+               EditBone *child = NULL;
+               EditBone *parent = NULL;
+               int boneStart = iter.start;
+
+               parent = add_editbone("Bone");
+               parent->flag = BONE_SELECTED|BONE_TIPSEL|BONE_ROOTSEL;
+               VECCOPY(parent->head, head->p);
+               
+               for (previous = nextBucket(&iter), bucket = nextBucket(&iter);
+                       bucket;
+                       previous = bucket, bucket = nextBucket(&iter))
+               {
+                       /* Calculate normal */
+                       VecSubf(n, bucket->p, parent->head);
+
+                       if (calcCorrelation(arc, boneStart, iter.index, parent->head, n) < CORRELATION_THRESHOLD)
+                       {
+                               VECCOPY(parent->tail, previous->p);
+
+                               child = add_editbone("Bone");
+                               VECCOPY(child->head, parent->tail);
+                               child->parent = parent;
+                               child->flag |= BONE_CONNECTED|BONE_SELECTED|BONE_TIPSEL|BONE_ROOTSEL;
+                               
+                               parent = child; // new child is next parent
+                               boneStart = iter.index; // start from end
+                       }
+               }
+
+               VECCOPY(parent->tail, tail->p);
+               
+               lastBone = parent; /* set last bone in the chain */
+       }
+       
+       return lastBone;
+}
+
+float arcLengthRatio(ReebArc *arc)
+{
+       float arcLength = 0.0f;
+       float embedLength = 0.0f;
+       int i;
+       
+       arcLength = VecLenf(arc->v1->p, arc->v2->p);
+       
+       if (arc->bcount > 0)
+       {
+               /* Add the embedding */
+               for ( i = 1; i < arc->bcount; i++)
+               {
+                       embedLength += VecLenf(arc->buckets[i - 1].p, arc->buckets[i].p);
+               }
+               /* Add head and tail -> embedding vectors */
+               embedLength += VecLenf(arc->v1->p, arc->buckets[0].p);
+               embedLength += VecLenf(arc->v2->p, arc->buckets[arc->bcount - 1].p);
+       }
+       else
+       {
+               embedLength = arcLength;
+       }
+       
+       return embedLength / arcLength; 
+}
+
+EditBone * subdivideByLength(ReebArc *arc, ReebNode *head, ReebNode *tail)
+{
+       EditBone *lastBone = NULL;
+       if ((G.scene->toolsettings->skgen_options & SKGEN_CUT_LENGTH) &&
+               arcLengthRatio(arc) >= G.scene->toolsettings->skgen_length_ratio)
+       {
+               ReebArcIterator iter;
+               EmbedBucket *bucket = NULL;
+               EmbedBucket *previous = NULL;
+               EditBone *child = NULL;
+               EditBone *parent = NULL;
+               float lengthLimit = G.scene->toolsettings->skgen_length_limit;
+               int same = 0;
+               
+               parent = add_editbone("Bone");
+               parent->flag |= BONE_SELECTED|BONE_TIPSEL|BONE_ROOTSEL;
+               VECCOPY(parent->head, head->p);
+
+               initArcIterator(&iter, arc, head);
+
+               bucket = nextBucket(&iter);
+               
+               while (bucket != NULL)
+               {
+                       float *vec0 = NULL;
+                       float *vec1 = bucket->p;
+
+                       /* first bucket. Previous is head */
+                       if (previous == NULL)
+                       {
+                               vec0 = head->p;
+                       }
+                       /* Previous is a valid bucket */
+                       else
+                       {
+                               vec0 = previous->p;
+                       }
+                       
+                       /* If lengthLimit hits the current segment */
+                       if (VecLenf(vec1, parent->head) > lengthLimit)
+                       {
+                               if (same == 0)
+                               {
+                                       float dv[3], off[3];
+                                       float a, b, c, f;
+                                       
+                                       /* Solve quadratic distance equation */
+                                       VecSubf(dv, vec1, vec0);
+                                       a = Inpf(dv, dv);
+                                       
+                                       VecSubf(off, vec0, parent->head);
+                                       b = 2 * Inpf(dv, off);
+                                       
+                                       c = Inpf(off, off) - (lengthLimit * lengthLimit);
+                                       
+                                       f = (-b + (float)sqrt(b * b - 4 * a * c)) / (2 * a);
+                                       
+                                       //printf("a %f, b %f, c %f, f %f\n", a, b, c, f);
+                                       
+                                       if (isnan(f) == 0 && f < 1.0f)
+                                       {
+                                               VECCOPY(parent->tail, dv);
+                                               VecMulf(parent->tail, f);
+                                               VecAddf(parent->tail, parent->tail, vec0);
+                                       }
+                                       else
+                                       {
+                                               VECCOPY(parent->tail, vec1);
+                                       }
+                               }
+                               else
+                               {
+                                       float dv[3];
+                                       
+                                       VecSubf(dv, vec1, vec0);
+                                       Normalize(dv);
+                                        
+                                       VECCOPY(parent->tail, dv);
+                                       VecMulf(parent->tail, lengthLimit);
+                                       VecAddf(parent->tail, parent->tail, parent->head);
+                               }
+                               
+                               child = add_editbone("Bone");
+                               VECCOPY(child->head, parent->tail);
+                               child->parent = parent;
+                               child->flag |= BONE_CONNECTED|BONE_SELECTED|BONE_TIPSEL|BONE_ROOTSEL;
+                               
+                               parent = child; // new child is next parent
+                               
+                               same = 1; // mark as same
+                       }
+                       else
+                       {
+                               previous = bucket;
+                               bucket = nextBucket(&iter);
+                               same = 0; // Reset same
+                       }
+               }
+               VECCOPY(parent->tail, tail->p);
+               
+               lastBone = parent; /* set last bone in the chain */
+       }
+       
+       return lastBone;
+}
+
+/***************************************** MAIN ALGORITHM ********************************************/
+
+void generateSkeletonFromReebGraph(ReebGraph *rg)
+{
+       GHash *arcBoneMap = NULL;
+       ReebArc *arc = NULL;
+       ReebNode *node = NULL;
+       Object *src = NULL;
+       Object *dst = NULL;
+       
+       src = BASACT->object;
+       
+       if (G.obedit != NULL)
+       {
+               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);
+       G.obedit= BASACT->object;
+       
+       /* Copy orientation from source */
+       VECCOPY(dst->loc, src->obmat[3]);
+       Mat4ToEul(src->obmat, dst->rot);
+       Mat4ToSize(src->obmat, dst->size);
+       
+       where_is_object(G.obedit);
+       
+       make_editArmature();
+
+       arcBoneMap = BLI_ghash_new(BLI_ghashutil_ptrhash, BLI_ghashutil_ptrcmp);
+       
+       markdownSymmetry(rg);
+       
+       for (arc = rg->arcs.first; arc; arc = arc->next) 
+       {
+               EditBone *lastBone = NULL;
+               ReebNode *head, *tail;
+               int i;
+
+               /* 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)
+                * 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)
+                */
+
+               /* if arc is a symmetry axis, internal bones go up the tree */          
+               if (arc->flags == 1 && arc->v2->degree != 1)
+               {
+                       head = arc->v2;
+                       tail = arc->v1;
+                       
+                       arc->flags = -1; /* mark arc direction */
+               }
+               /* Bones point AWAY from the symmetry axis */
+               else if (arc->v1->flags == 1)
+               {
+                       head = arc->v1;
+                       tail = arc->v2;
+                       
+                       arc->flags = 1; /* mark arc direction */
+               }
+               else if (arc->v2->flags == 1)
+               {
+                       head = arc->v2;
+                       tail = arc->v1;
+                       
+                       arc->flags = -1; /* mark arc direction */
+               }
+               /* otherwise, always go from low weight to high weight */
+               else
+               {
+                       head = arc->v1;
+                       tail = arc->v2;
+                       
+                       arc->flags = 1; /* mark arc direction */
+               }
+               
+               /* Loop over subdivision methods */     
+               for (i = 0; lastBone == NULL && i < SKGEN_SUB_TOTAL; i++)
+               {
+                       switch(G.scene->toolsettings->skgen_subdivisions[i])
+                       {
+                               case SKGEN_SUB_LENGTH:
+                                       lastBone = subdivideByLength(arc, head, tail);
+                                       break;
+                               case SKGEN_SUB_ANGLE:
+                                       lastBone = subdivideByAngle(arc, head, tail);
+                                       break;
+                               case SKGEN_SUB_CORRELATION:
+                                       lastBone = subdivideByCorrelation(arc, head, tail);
+                                       break;
+                       }
+               }
+       
+               if (lastBone == NULL)
+               {
+                       EditBone        *bone;
+                       bone = add_editbone("Bone");
+                       bone->flag |= BONE_SELECTED|BONE_TIPSEL|BONE_ROOTSEL;
+                       
+                       VECCOPY(bone->head, head->p);
+                       VECCOPY(bone->tail, tail->p);
+                       
+                       /* set first and last bone, since there's only one */
+                       lastBone = bone;
+               }
+               
+               BLI_ghash_insert(arcBoneMap, arc, lastBone);
+       }
+
+       /* Second pass, setup parent relationship between arcs */
+       for (node = rg->nodes.first; node; node = node->next)
+       {
+               ReebArc *incomingArc = NULL;
+               int i;
+
+               for (i = 0; node->arcs[i] != NULL; i++)
+               {
+                       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 (incomingArc == NULL)
+                               {
+                                       incomingArc = arc;
+                                       /* loop further to make sure there's only one incoming arc */
+                               }
+                               else
+                               {
+                                       /* skip this node if more than one incomingArc */
+                                       incomingArc = NULL;
+                                       break; /* No need to look further, we are skipping already */
+                               }
+                       }
+               }
+
+               if (incomingArc != NULL)
+               {
+                       EditBone *parentBone = BLI_ghash_lookup(arcBoneMap, incomingArc);
+
+                       /* Look for outgoing arcs and parent their bones */
+                       for (i = 0; node->arcs[i] != NULL; i++)
+                       {
+                               arc = node->arcs[i];
+
+                               /* if arc is outgoing from the node */
+                               if ((arc->v1 == node && arc->flags == 1) || (arc->v2 == node && arc->flags == -1))
+                               {
+                                       EditBone *childBone = BLI_ghash_lookup(arcBoneMap, arc);
+
+                                       /* find the root bone */
+                                       while(childBone->parent != NULL)
+                                       {
+                                               childBone = childBone->parent;
+                                       }
+
+                                       childBone->parent = parentBone;
+                                       childBone->flag |= BONE_CONNECTED;
+                               }
+                       }
+               }
+       }
+       
+       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;
+       
+       if (em == NULL)
+               return;
+
+       setcursor_space(SPACE_VIEW3D, CURSOR_WAIT);
+
+       if (weightFromDistance(em) == 0)
+       {
+               error("No selected vertex\n");
+               return;
+       }
+       
+       weightToHarmonic(em);
+               
+       renormalizeWeight(em, 1.0f);
+
+//#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);
+
+       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);
+       
+       generateSkeletonFromReebGraph(rg);
 
+       freeGraph(rg);
+}
diff --git a/source/blender/src/reeb.c b/source/blender/src/reeb.c
new file mode 100644 (file)
index 0000000..85fb581
--- /dev/null
@@ -0,0 +1,1923 @@
+/**
+ * $Id: 
+ *
+ * ***** BEGIN GPL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version. The Blender
+ * Foundation also sells licenses for use in proprietary software under
+ * the Blender License.  See http://www.blender.org/BL/ for information
+ * about this.
+ *
+ * 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 *****
+ */
+#include <math.h>
+#include <string.h> // for memcpy
+#include <stdio.h>
+#include <stdlib.h> // for qsort
+
+#include "DNA_listBase.h"
+#include "DNA_scene_types.h"
+#include "DNA_space_types.h"
+#include "DNA_meshdata_types.h"
+
+#include "MEM_guardedalloc.h"
+
+#include "BLI_blenlib.h"
+#include "BLI_arithb.h"
+#include "BLI_editVert.h"
+#include "BLI_edgehash.h"
+
+#include "BDR_editobject.h"
+
+#include "BIF_editmesh.h"
+#include "BIF_editarmature.h"
+#include "BIF_interface.h"
+#include "BIF_toolbox.h"
+#include "BIF_graphics.h"
+
+#include "BKE_global.h"
+#include "BKE_utildefines.h"
+#include "BKE_customdata.h"
+
+#include "blendef.h"
+
+#include "ONL_opennl.h"
+
+#include "reeb.h"
+
+/*
+ * Skeleton generation algorithm based on: 
+ * "Harmonic Skeleton for Realistic Character Animation"
+ * Gregoire Aujay, Franck Hetroy, Francis Lazarus and Christine Depraz
+ * SIGGRAPH 2007
+ * 
+ * Reeb graph generation algorithm based on: 
+ * "Robust On-line Computation of Reeb Graphs: Simplicity and Speed"
+ * Valerio Pascucci, Giorgio Scorzelli, Peer-Timo Bremer and Ajith Mascarenhas
+ * SIGGRAPH 2007
+ * 
+ * */
+
+int mergeArcs(ReebGraph *rg, ReebArc *a0, ReebArc *a1);
+int mergeConnectedArcs(ReebGraph *rg, ReebArc *a0, ReebArc *a1);
+EditEdge * NextEdgeForVert(EditMesh *em, EditVert *v);
+
+/***************************************** BUCKET UTILS **********************************************/
+
+void addVertToBucket(EmbedBucket *b, float co[3])
+{
+       b->nv++;
+       VecLerpf(b->p, b->p, co, 1.0f / b->nv);
+}
+
+void removeVertFromBucket(EmbedBucket *b, float co[3])
+{
+       VecMulf(b->p, (float)b->nv);
+       VecSubf(b->p, b->p, co);
+       b->nv--;
+       VecMulf(b->p, 1.0f / (float)b->nv);
+}
+
+void mergeBuckets(EmbedBucket *bDst, EmbedBucket *bSrc)
+{
+       if (bDst->nv > 0 && bSrc->nv > 0)
+       {
+               bDst->nv += bSrc->nv;
+               VecLerpf(bDst->p, bDst->p, bSrc->p, (float)bSrc->nv / (float)(bDst->nv));
+       }
+       else if (bSrc->nv > 0)
+       {
+               bDst->nv = bSrc->nv;
+               VECCOPY(bDst->p, bSrc->p);
+       }
+}
+
+void mergeArcBuckets(ReebArc *aDst, ReebArc *aSrc, float start, float end)
+{
+       if (aDst->bcount > 0 && aSrc->bcount > 0)
+       {
+               int indexDst = 0, indexSrc = 0;
+               
+               start = MAX3(start, aDst->buckets[0].val, aSrc->buckets[0].val);
+               
+               while(indexDst < aDst->bcount && aDst->buckets[indexDst].val < start)
+               {
+                       indexDst++;
+               }
+
+               while(indexSrc < aSrc->bcount && aSrc->buckets[indexSrc].val < start)
+               {
+                       indexSrc++;
+               }
+               
+               for( ;  indexDst < aDst->bcount &&
+                               indexSrc < aSrc->bcount &&
+                               aDst->buckets[indexDst].val <= end &&
+                               aSrc->buckets[indexSrc].val <= end
+                               
+                        ;      indexDst++, indexSrc++)
+               {
+                       mergeBuckets(aDst->buckets + indexDst, aSrc->buckets + indexSrc);
+               }
+       }
+}
+
+void allocArcBuckets(ReebArc *arc)
+{
+       int i;
+       float start = ceil(arc->v1->weight);
+       arc->bcount = (int)(floor(arc->v2->weight) - start) + 1;
+       
+       if (arc->bcount > 0)
+       {
+               arc->buckets = MEM_callocN(sizeof(EmbedBucket) * arc->bcount, "embed bucket");
+               
+               for(i = 0; i < arc->bcount; i++)
+               {
+                       arc->buckets[i].val = start + i;
+               }
+       }
+       else
+       {
+               arc->buckets = NULL;
+       }
+       
+}
+
+void resizeArcBuckets(ReebArc *arc)
+{
+       EmbedBucket *oldBuckets = arc->buckets;
+       int oldBCount = arc->bcount;
+       
+       allocArcBuckets(arc);
+       
+       if (oldBCount != 0 && arc->bcount != 0)
+       {
+               int oldStart = (int)oldBuckets[0].val;
+               int oldEnd = (int)oldBuckets[oldBCount - 1].val;
+               int newStart = (int)arc->buckets[0].val;
+               int newEnd = (int)arc->buckets[arc->bcount - 1].val;
+               int oldOffset = 0;
+               int newOffset = 0;
+               int len;
+               
+               if (oldStart < newStart)
+               {
+                       oldOffset = newStart - oldStart;
+               }
+               else
+               {
+                       newOffset = oldStart - newStart;
+               }
+               
+               len = MIN2(oldEnd - (oldStart + oldOffset) + 1, newEnd - (newStart - newOffset) + 1);
+               
+               memcpy(arc->buckets + newOffset, oldBuckets + oldOffset, len * sizeof(EmbedBucket)); 
+       }
+
+       if (oldBuckets != NULL)
+       {
+               MEM_freeN(oldBuckets);
+       }
+}
+/***************************************** UTILS **********************************************/
+
+ReebEdge * copyEdge(ReebEdge *edge)
+{
+       ReebEdge *newEdge = NULL;
+       
+       newEdge = MEM_callocN(sizeof(ReebEdge), "reeb edge");
+       memcpy(newEdge, edge, sizeof(ReebEdge));
+       
+       newEdge->next = NULL;
+       newEdge->prev = NULL;
+       
+       return newEdge;
+}
+
+void printArc(ReebArc *arc)
+{
+       ReebEdge *edge;
+       printf("arc: (%i)%f -> (%i)%f\n", arc->v1->index, arc->v1->weight, arc->v2->index, arc->v2->weight);
+       
+       for(edge = arc->edges.first; edge ; edge = edge->next)
+       {
+               printf("\tedge (%i, %i)\n", edge->v1->index, edge->v2->index);
+       }
+}
+
+void freeArc(ReebArc *arc)
+{
+       BLI_freelistN(&arc->edges);
+       
+       if (arc->buckets)
+               MEM_freeN(arc->buckets);
+       
+       MEM_freeN(arc);
+}
+
+void freeGraph(ReebGraph *rg)
+{
+       ReebArc *arc;
+       ReebNode *node;
+       
+       // free nodes
+       for( node = rg->nodes.first; node; node = node->next )
+       {
+               // Free adjacency lists
+               if (node->arcs != NULL)
+               {
+                       MEM_freeN(node->arcs);
+               }
+       }
+       BLI_freelistN(&rg->nodes);
+       
+       // free arcs
+       arc = rg->arcs.first;
+       while( arc )
+       {
+               ReebArc *next = arc->next;
+               freeArc(arc);
+               arc = next;
+       }
+       
+       // free edge map
+       BLI_edgehash_free(rg->emap, NULL);
+       
+       MEM_freeN(rg);
+}
+
+void repositionNodes(ReebGraph *rg)
+{
+       ReebArc *arc = NULL;
+       ReebNode *node = NULL;
+       
+       // Reset node positions
+       for(node = rg->nodes.first; node; node = node->next)
+       {
+               node->p[0] = node->p[1] = node->p[2] = 0;
+       }
+       
+       for(arc = rg->arcs.first; arc; arc = arc->next)
+       {
+               if (arc->bcount > 0)
+               {
+                       float p[3];
+                       
+                       VECCOPY(p, arc->buckets[0].p);
+                       VecMulf(p, 1.0f / arc->v1->degree);
+                       VecAddf(arc->v1->p, arc->v1->p, p);
+                       
+                       VECCOPY(p, arc->buckets[arc->bcount - 1].p);
+                       VecMulf(p, 1.0f / arc->v2->degree);
+                       VecAddf(arc->v2->p, arc->v2->p, p);
+               }
+       }
+}
+
+void verifyNodeDegree(ReebGraph *rg)
+{
+       ReebNode *node = NULL;
+       ReebArc *arc = NULL;
+
+       for(node = rg->nodes.first; node; node = node->next)
+       {
+               int count = 0;
+               for(arc = rg->arcs.first; arc; arc = arc->next)
+               {
+                       if (arc->v1 == node || arc->v2 == node)
+                       {
+                               count++;
+                       }
+               }
+               if (count != node->degree)
+               {
+                       printf("degree error in node %i: expected %i got %i\n", node->index, count, node->degree);
+               }
+       }
+}
+
+void verifyBuckets(ReebGraph *rg)
+{
+#ifdef DEBUG_REEB
+       ReebArc *arc = NULL;
+       for(arc = rg->arcs.first; arc; arc = arc->next)
+       {
+               if (arc->bcount > 0)
+               {
+                       int i;
+                       for(i = 0; i < arc->bcount; i++)
+                       {
+                               if (arc->buckets[i].nv == 0)
+                               {
+                                       printArc(arc);
+                                       printf("count error in bucket %i/%i\n", i+1, arc->bcount);
+                               }
+                       }
+                       
+                       if (ceil(arc->v1->weight) < arc->buckets[0].val)
+                       {
+                               printArc(arc);
+                               printf("alloc error in first bucket: %f should be %f \n", arc->buckets[0].val, ceil(arc->v1->weight));
+                       }
+                       if (floor(arc->v2->weight) < arc->buckets[arc->bcount - 1].val)
+                       {
+                               printArc(arc);
+                               printf("alloc error in last bucket: %f should be %f \n", arc->buckets[arc->bcount - 1].val, floor(arc->v2->weight));
+                       }
+               }
+       }
+#endif
+}
+
+/************************************** ADJACENCY LIST *************************************************/
+
+void addArcToNodeAdjacencyList(ReebNode *node, ReebArc *arc)
+{
+       ReebArc **arclist;
+
+       for(arclist = node->arcs; *arclist; arclist++)
+       {       }
+       
+       *arclist = arc;
+}
+
+void buildAdjacencyList(ReebGraph *rg)
+{
+       ReebNode *node = NULL;
+       ReebArc *arc = NULL;
+
+       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(ReebArc*), "adjacency list");
+       }
+
+       for(arc = rg->arcs.first; arc; arc= arc->next)
+       {
+               addArcToNodeAdjacencyList(arc->v1, arc);
+               addArcToNodeAdjacencyList(arc->v2, arc);
+       }
+}
+
+int hasAdjacencyList(ReebGraph *rg)
+{
+       ReebNode *node;
+       
+       for(node = rg->nodes.first; node; node = node->next)
+       {
+               if (node->arcs == NULL)
+               {
+                       return 0;
+               }
+       }
+       
+       return 1;
+}
+
+int countConnectedArcs(ReebGraph *rg, ReebNode *node)
+{
+       int count = 0;
+       
+       /* use adjacency list if present */
+       if (node->arcs)
+       {
+               ReebArc **arcs;
+       
+               for(arcs = node->arcs; *arcs; arcs++)
+               {
+                       count++;
+               }
+       }
+       else
+       {
+               ReebArc *arc;
+               for(arc = rg->arcs.first; arc; arc = arc->next)
+               {
+                       if (arc->v1 == node || arc->v2 == node)
+                       {
+                               count++;
+                       }
+               }
+       }
+       
+       return count;
+}
+
+/****************************************** SMOOTHING **************************************************/
+
+void postprocessGraph(ReebGraph *rg, char mode)
+{
+       ReebArc *arc;
+       float fac1 = 0, fac2 = 1, fac3 = 0;
+
+       switch(mode)
+       {
+       case SKGEN_AVERAGE:
+               fac1 = fac2 = fac3 = 1.0f / 3.0f;
+               break;
+       case SKGEN_SMOOTH:
+               fac1 = fac3 = 0.25f;
+               fac2 = 0.5f;
+               break;
+       case SKGEN_SHARPEN:
+               fac1 = fac2 = -0.25f;
+               fac2 = 1.5f;
+               break;
+       default:
+               error("Unknown post processing mode");
+               return;
+       }
+       
+       for(arc = rg->arcs.first; arc; arc = arc->next)
+       {
+               EmbedBucket *buckets = arc->buckets;
+               int bcount = arc->bcount;
+               int index;
+
+               for(index = 1; index < bcount - 1; index++)
+               {
+                       VecLerpf(buckets[index].p, buckets[index].p, buckets[index - 1].p, fac1 / (fac1 + fac2));
+                       VecLerpf(buckets[index].p, buckets[index].p, buckets[index + 1].p, fac3 / (fac1 + fac2 + fac3));
+               }
+       }
+}
+
+/********************************************SORTING****************************************************/
+
+int compareNodesWeight(void *vnode1, void *vnode2)
+{
+       ReebNode *node1 = (ReebNode*)vnode1;
+       ReebNode *node2 = (ReebNode*)vnode2;
+       
+       if (node1->weight < node2->weight)
+       {
+               return -1;
+       }
+       if (node1->weight > node2->weight)
+       {
+               return 1;
+       }
+       else
+       {
+               return 0;
+       }
+}
+
+void sortNodes(ReebGraph *rg)
+{
+       BLI_sortlist(&rg->nodes, compareNodesWeight);
+}
+
+int compareArcsWeight(void *varc1, void *varc2)
+{
+       ReebArc *arc1 = (ReebArc*)varc1;
+       ReebArc *arc2 = (ReebArc*)varc2;
+       
+       if (arc1->v1->weight < arc2->v1->weight)
+       {
+               return -1;
+       }
+       if (arc1->v1->weight > arc2->v1->weight)
+       {
+               return 1;
+       }
+       else
+       {
+               return 0;
+       }
+}
+
+void sortArcs(ReebGraph *rg)
+{
+       BLI_sortlist(&rg->arcs, compareArcsWeight);
+}
+
+/****************************************** FILTERING **************************************************/
+
+int compareArcs(void *varc1, void *varc2)
+{
+       ReebArc *arc1 = (ReebArc*)varc1;
+       ReebArc *arc2 = (ReebArc*)varc2;
+       float len1 = arc1->v2->weight - arc1->v1->weight;
+       float len2 = arc2->v2->weight - arc2->v1->weight;
+       
+       if (len1 < len2)
+       {
+               return -1;
+       }
+       if (len1 > len2)
+       {
+               return 1;
+       }
+       else
+       {
+               return 0;
+       }
+}
+
+void filterArc(ReebGraph *rg, ReebNode *newNode, ReebNode *removedNode, ReebArc * srcArc, int merging)
+{
+       ReebArc *arc = NULL, *nextArc = NULL;
+
+       /* first pass, merge buckets for arcs that spawned the two nodes into the source arc*/
+       for(arc = rg->arcs.first; arc; arc = arc->next)
+       {
+               if (arc->v1 == srcArc->v1 && arc->v2 == srcArc->v2 && arc != srcArc)
+               {
+                       mergeArcBuckets(srcArc, arc, srcArc->v1->weight, srcArc->v2->weight);
+               }
+       }
+
+       /* second pass, replace removedNode by newNode, remove arcs that are collapsed in a loop */
+       arc = rg->arcs.first;
+       while(arc)
+       {
+               nextArc = arc->next;
+               
+               if (arc->v1 == removedNode || arc->v2 == removedNode)
+               {
+                       if (arc->v1 == removedNode)
+                       {
+                               arc->v1 = newNode;
+                       }
+                       else
+                       {
+                               arc->v2 = newNode;
+                       }
+
+                       // Remove looped arcs                   
+                       if (arc->v1 == arc->v2)
+                       {
+                               // v1 or v2 was already newNode, since we're removing an arc, decrement degree
+                               newNode->degree--;
+                               
+                               // If it's safeArc, it'll be removed later, so keep it for now
+                               if (arc != srcArc)
+                               {
+                                       BLI_remlink(&rg->arcs, arc);
+                                       freeArc(arc);
+                               }
+                       }
+                       // Remove flipped arcs
+                       else if (arc->v1->weight > arc->v2->weight)
+                       {
+                               // Decrement degree from the other node
+                               OTHER_NODE(arc, newNode)->degree--;
+                               
+                               BLI_remlink(&rg->arcs, arc);
+                               freeArc(arc);
+                       }
+                       else
+                       {
+                               newNode->degree++; // incrementing degree since we're adding an arc
+
+                               if (merging)
+                               {
+                                       // resize bucket list
+                                       resizeArcBuckets(arc);
+                                       mergeArcBuckets(arc, srcArc, arc->v1->weight, arc->v2->weight);
+                               }
+                       }
+               }
+               
+               arc = nextArc;
+       }
+}
+
+void filterNullReebGraph(ReebGraph *rg)
+{
+       ReebArc *arc = NULL, *nextArc = NULL;
+       
+       arc = rg->arcs.first;
+       while(arc)
+       {
+               nextArc = arc->next;
+               // Only collapse arcs too short to have any embed bucket
+               if (arc->bcount == 0)
+               {
+                       ReebNode *newNode = arc->v1;
+                       ReebNode *removedNode = arc->v2;
+                       float blend;
+                       
+                       blend = (float)newNode->degree / (float)(newNode->degree + removedNode->degree); // blending factors
+                       
+                       //newNode->weight = FloatLerpf(newNode->weight, removedNode->weight, blend);
+                       VecLerpf(newNode->p, newNode->p, removedNode->p, blend);
+                       
+                       filterArc(rg, newNode, removedNode, arc, 0);
+
+                       // Reset nextArc, it might have changed
+                       nextArc = arc->next;
+                       
+                       BLI_remlink(&rg->arcs, arc);
+                       freeArc(arc);
+                       
+                       BLI_freelinkN(&rg->nodes, removedNode);
+               }
+               
+               arc = nextArc;
+       }
+}
+
+int filterInternalReebGraph(ReebGraph *rg, float threshold)
+{
+       ReebArc *arc = NULL, *nextArc = NULL;
+       int value = 0;
+       
+       BLI_sortlist(&rg->arcs, compareArcs);
+
+       arc = rg->arcs.first;
+       while(arc)
+       {
+               nextArc = arc->next;
+
+               // Only collapse non-terminal arcs that are shorter than threshold
+               if ((arc->v1->degree > 1 && arc->v2->degree > 1 && arc->v2->weight - arc->v1->weight < threshold))
+               {
+                       ReebNode *newNode = NULL;
+                       ReebNode *removedNode = NULL;
+                       
+                       /* Keep the node with the highestn number of connected arcs */
+                       if (arc->v1->degree >= arc->v2->degree)
+                       {
+                               newNode = arc->v1;
+                               removedNode = arc->v2;
+                       }
+                       else
+                       {
+                               newNode = arc->v2;
+                               removedNode = arc->v1;
+                       }
+                       
+                       filterArc(rg, newNode, removedNode, arc, 1);
+
+                       // 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;
+}
+
+int filterExternalReebGraph(ReebGraph *rg, float threshold)
+{
+       ReebArc *arc = NULL, *nextArc = NULL;
+       int value = 0;
+       
+       BLI_sortlist(&rg->arcs, compareArcs);
+
+       arc = rg->arcs.first;
+       while(arc)
+       {
+               nextArc = arc->next;
+
+               // Only collapse terminal arcs that are shorter than threshold
+               if ((arc->v1->degree == 1 || arc->v2->degree == 1) && arc->v2->weight - arc->v1->weight < threshold)
+               {
+                       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 )
+{
+       EditVert *va = *(EditVert**)a;
+       EditVert *vb = *(EditVert**)b;
+       int value = 0;
+       
+       if (va->tmp.fp < vb->tmp.fp)
+       {
+               value = -1;
+       }
+       else if (va->tmp.fp > vb->tmp.fp)
+       {
+               value = 1;
+       }
+
+       return value;           
+}
+
+void spreadWeight(EditMesh *em)
+{
+       EditVert **verts, *eve;
+       float lastWeight = 0.0f;
+       int totvert = BLI_countlist(&em->verts);
+       int i;
+       int work_needed = 1;
+       
+       verts = MEM_callocN(sizeof(EditVert*) * totvert, "verts array");
+       
+       for(eve = em->verts.first, i = 0; eve; eve = eve->next, i++)
+       {
+               verts[i] = eve;
+       }
+       
+       while(work_needed == 1)
+       {
+               work_needed = 0;
+               qsort(verts, totvert, sizeof(EditVert*), compareVerts);
+               
+               for(i = 0; i < totvert; i++)
+               {
+                       eve = verts[i];
+                       
+                       if (i == 0 || (eve->tmp.fp - lastWeight) > FLT_EPSILON)
+                       {
+                               lastWeight = eve->tmp.fp;
+                       }
+                       else
+                       {
+                               work_needed = 1;
+                               eve->tmp.fp = lastWeight + FLT_EPSILON * 2;
+                               lastWeight = eve->tmp.fp;
+                       }
+               }
+       }
+       
+       MEM_freeN(verts);
+}
+/*********************************** GRAPH AS TREE FUNCTIONS *******************************************/
+
+int subtreeDepth(ReebNode *node, ReebArc *rootArc)
+{
+       int depth = 0;
+       
+       /* Base case, no arcs leading away */
+       if (node->arcs == NULL || *(node->arcs) == NULL)
+       {
+               return 0;
+       }
+       else
+       {
+               ReebArc ** pArc;
+
+               for(pArc = node->arcs; *pArc; pArc++)
+               {
+                       ReebArc *arc = *pArc;
+                       
+                       /* only arcs that go down the tree */
+                       if (arc != rootArc)
+                       {
+                               ReebNode *newNode = OTHER_NODE(arc, node);
+                               depth = MAX2(depth, subtreeDepth(newNode, arc));
+                       }
+               }
+       }
+       
+       return depth + 1;
+}
+
+/*************************************** CYCLE DETECTION ***********************************************/
+
+int detectCycle(ReebNode *node, ReebArc *srcArc)
+{
+       int value = 0;
+       
+       if (node->flags == 0)
+       {
+               ReebArc ** pArc;
+
+               /* mark node as visited */
+               node->flags = 1;
+
+               for(pArc = node->arcs; *pArc && value == 0; pArc++)
+               {
+                       ReebArc *arc = *pArc;
+                       
+                       /* don't go back on the source arc */
+                       if (arc != srcArc)
+                       {
+                               value = detectCycle(OTHER_NODE(arc, node), arc);
+                       }
+               }
+       }
+       else
+       {
+               value = 1;
+       }
+       
+       return value;
+}
+
+int    isGraphCyclic(ReebGraph *rg)
+{
+       ReebNode *node;
+       int value = 0;
+       
+       /* NEED TO CHECK IF ADJACENCY LIST EXIST */
+       
+       /* Mark all nodes as not visited */
+       for(node = rg->nodes.first; node; node = node->next)
+       {
+               node->flags = 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)
+               {
+                       value = value || detectCycle(node, NULL);
+               }               
+       }
+       
+       return value;
+}
+
+/******************************************** EXPORT ***************************************************/
+
+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)
+{
+#ifdef DEBUG_REEB
+       ReebArc *arc;
+       char filename[128];
+       FILE *f;
+       
+       if (count == -1)
+       {
+               sprintf(filename, "test.txt");
+       }
+       else
+       {
+               sprintf(filename, "test%05i.txt", count);
+       }
+       f = fopen(filename, "w");
+
+       for(arc = rg->arcs.first; arc; arc = arc->next)
+       {
+               int i;
+               
+               exportNode(f, "v1", arc->v1);
+               
+               for(i = 0; i < arc->bcount; i++)
+               {
+                       fprintf(f, "b nv:%i %f %f %f\n", arc->buckets[i].nv, arc->buckets[i].p[0], arc->buckets[i].p[1], arc->buckets[i].p[2]);
+               }
+               
+               exportNode(f, "v2", arc->v2);
+       }       
+       
+       fclose(f);
+#endif
+}
+
+/***************************************** MAIN ALGORITHM **********************************************/
+
+ReebArc * findConnectedArc(ReebGraph *rg, ReebArc *arc, ReebNode *v)
+{
+       ReebArc *nextArc = arc->next;
+       
+       for(nextArc = rg->arcs.first; nextArc; nextArc = nextArc->next)
+       {
+               if (arc != nextArc && (nextArc->v1 == v || nextArc->v2 == v))
+               {
+                       break;
+               }
+       }
+       
+       return nextArc;
+}
+
+void removeNormalNodes(ReebGraph *rg)
+{
+       ReebArc *arc;
+       
+       // Merge degree 2 nodes
+       for(arc = rg->arcs.first; arc; arc = arc->next)
+       {
+               while (arc->v1->degree == 2 || arc->v2->degree == 2)
+               {
+                       // merge at v1
+                       if (arc->v1->degree == 2)
+                       {
+                               ReebArc *nextArc = findConnectedArc(rg, arc, arc->v1);
+
+                               // Merge arc only if needed
+                               if (arc->v1 == nextArc->v2)
+                               {                               
+                                       mergeConnectedArcs(rg, arc, nextArc);
+                               }
+                               // Otherwise, mark down vert
+                               else
+                               {
+                                       arc->v1->degree = 3;
+                               }
+                       }
+                       
+                       // merge at v2
+                       if (arc->v2->degree == 2)
+                       {
+                               ReebArc *nextArc = findConnectedArc(rg, arc, arc->v2);
+                               
+                               // Merge arc only if needed
+                               if (arc->v2 == nextArc->v1)
+                               {                               
+                                       mergeConnectedArcs(rg, arc, nextArc);
+                               }
+                               // Otherwise, mark down vert
+                               else
+                               {
+                                       arc->v2->degree = 3;
+                               }
+                       }
+               }
+       }
+       
+}
+
+int edgeEquals(ReebEdge *e1, ReebEdge *e2)
+{
+       return (e1->v1 == e2->v1 && e1->v2 == e2->v2);
+}
+
+ReebArc *nextArcMappedToEdge(ReebArc *arc, ReebEdge *e)
+{
+       ReebEdge *nextEdge = NULL;
+       ReebEdge *edge = NULL;
+       ReebArc *result = NULL;
+
+       /* Find the ReebEdge in the edge list */
+       for(edge = arc->edges.first; edge && !edgeEquals(edge, e); edge = edge->next)
+       {       }
+       
+       nextEdge = edge->nextEdge;
+       
+       if (nextEdge != NULL)
+       {
+               result = nextEdge->arc;
+       }
+
+       return result;
+}
+
+typedef enum {
+       MERGE_LOWER,
+       MERGE_HIGHER,
+       MERGE_APPEND
+} MergeDirection;
+
+void mergeArcEdges(ReebGraph *rg, ReebArc *aDst, ReebArc *aSrc, MergeDirection direction)
+{
+       ReebEdge *e = NULL;
+       
+       if (direction == MERGE_APPEND)
+       {
+               for(e = aSrc->edges.first; e; e = e->next)
+               {
+                       e->arc = aDst; // Edge is stolen by new arc
+               }
+               
+               addlisttolist(&aDst->edges , &aSrc->edges);
+       }
+       else
+       {
+               for(e = aSrc->edges.first; e; e = e->next)
+               {
+                       ReebEdge *newEdge = copyEdge(e);
+
+                       newEdge->arc = aDst;
+                       
+                       BLI_addtail(&aDst->edges, newEdge);
+                       
+                       if (direction == MERGE_LOWER)
+                       {
+                               void **p = BLI_edgehash_lookup_p(rg->emap, e->v1->index, e->v2->index);
+                               
+                               newEdge->nextEdge = e;
+
+                               // if edge was the first in the list, point the edit edge to the new reeb edge instead.                                                 
+                               if (*p == e)
+                               {
+                                       *p = (void*)newEdge;
+                               }
+                               // otherwise, advance in the list until the predecessor is found then insert it there
+                               else
+                               {
+                                       ReebEdge *previous = (ReebEdge*)*p;
+                                       
+                                       while(previous->nextEdge != e)
+                                       {
+                                               previous = previous->nextEdge;
+                                       }
+                                       
+                                       previous->nextEdge = newEdge;
+                               }
+                       }
+                       else
+                       {
+                               newEdge->nextEdge = e->nextEdge;
+                               e->nextEdge = newEdge;
+                       }
+               }
+       }
+} 
+
+// return 1 on full merge
+int mergeConnectedArcs(ReebGraph *rg, ReebArc *a0, ReebArc *a1)
+{
+       int result = 0;
+       ReebNode *removedNode = NULL;
+       
+       mergeArcEdges(rg, a0, a1, MERGE_APPEND);
+       
+       // Bring a0 to the combine length of both arcs
+       if (a0->v2 == a1->v1)
+       {
+               removedNode = a0->v2;
+               a0->v2 = a1->v2;
+       }
+       else if (a0->v1 == a1->v2)
+       {
+               removedNode = a0->v1;
+               a0->v1 = a1->v1;
+       }
+       
+       resizeArcBuckets(a0);
+       // Merge a1 in a0
+       mergeArcBuckets(a0, a1, a0->v1->weight, a0->v2->weight);
+       
+       // remove a1 from graph
+       BLI_remlink(&rg->arcs, a1);
+       freeArc(a1);
+       
+       BLI_freelinkN(&rg->nodes, removedNode);
+       result = 1;
+       
+       return result;
+}
+// return 1 on full merge
+int mergeArcs(ReebGraph *rg, ReebArc *a0, ReebArc *a1)
+{
+       int result = 0;
+       // TRIANGLE POINTS DOWN
+       if (a0->v1->weight == a1->v1->weight) // heads are the same
+       {
+               if (a0->v2->weight == a1->v2->weight) // tails also the same, arcs can be totally merge together
+               {
+                       mergeArcEdges(rg, a0, a1, MERGE_APPEND);
+                       
+                       mergeArcBuckets(a0, a1, a0->v1->weight, a0->v2->weight);
+                       
+                       // Adjust node degree
+                       a1->v1->degree--;
+                       a1->v2->degree--;
+                       
+                       // remove a1 from graph
+                       BLI_remlink(&rg->arcs, a1);
+                       
+                       freeArc(a1);
+                       result = 1;
+               }
+               else if (a0->v2->weight > a1->v2->weight) // a1->v2->weight is in the middle
+               {
+                       mergeArcEdges(rg, a1, a0, MERGE_LOWER);
+
+                       // Adjust node degree
+                       a0->v1->degree--;
+                       a1->v2->degree++;
+                       
+                       mergeArcBuckets(a1, a0, a1->v1->weight, a1->v2->weight);
+                       a0->v1 = a1->v2;
+                       resizeArcBuckets(a0);
+               }
+               else // a0>n2 is in the middle
+               {
+                       mergeArcEdges(rg, a0, a1, MERGE_LOWER);
+                       
+                       // Adjust node degree
+                       a1->v1->degree--;
+                       a0->v2->degree++;
+                       
+                       mergeArcBuckets(a0, a1, a0->v1->weight, a0->v2->weight);
+                       a1->v1 = a0->v2;
+                       resizeArcBuckets(a1);
+               }
+       }
+       // TRIANGLE POINTS UP
+       else if (a0->v2->weight == a1->v2->weight) // tails are the same
+       {
+               if (a0->v1->weight > a1->v1->weight) // a0->v1->weight is in the middle
+               {
+                       mergeArcEdges(rg, a0, a1, MERGE_HIGHER);
+                       
+                       // Adjust node degree
+                       a1->v2->degree--;
+                       a0->v1->degree++;
+                       
+                       mergeArcBuckets(a0, a1, a0->v1->weight, a0->v2->weight);
+                       a1->v2 = a0->v1;
+                       resizeArcBuckets(a1);
+               }
+               else // a1->v1->weight is in the middle
+               {
+                       mergeArcEdges(rg, a1, a0, MERGE_HIGHER);
+
+                       // Adjust node degree
+                       a0->v2->degree--;
+                       a1->v1->degree++;
+
+                       mergeArcBuckets(a1, a0, a1->v1->weight, a1->v2->weight);
+                       a0->v2 = a1->v1;
+                       resizeArcBuckets(a0);
+               }
+       }
+       else
+       {
+               // Need something here (OR NOT)
+       }
+       
+       return result;
+}
+
+void glueByMergeSort(ReebGraph *rg, ReebArc *a0, ReebArc *a1, ReebEdge *e0, ReebEdge *e1)
+{
+       int total = 0;
+       while (total == 0 && a0 != a1 && a0 != NULL && a1 != NULL)
+       {
+               total = mergeArcs(rg, a0, a1);
+               
+               if (total == 0) // if it wasn't a total merge, go forward
+               {
+                       if (a0->v2->weight < a1->v2->weight)
+                       {
+                               a0 = nextArcMappedToEdge(a0, e0);
+                       }
+                       else
+                       {
+                               a1 = nextArcMappedToEdge(a1, e1);
+                       }
+               }
+       }
+}
+
+void mergePaths(ReebGraph *rg, ReebEdge *e0, ReebEdge *e1, ReebEdge *e2)
+{
+       ReebArc *a0, *a1, *a2;
+       a0 = e0->arc;
+       a1 = e1->arc;
+       a2 = e2->arc;
+       
+       glueByMergeSort(rg, a0, a1, e0, e1);
+       glueByMergeSort(rg, a0, a2, e0, e2);
+} 
+
+ReebNode * addNode(ReebGraph *rg, EditVert *eve, float weight)
+{
+       ReebNode *node = NULL;
+       
+       node = MEM_callocN(sizeof(ReebNode), "reeb node");
+       
+       node->flags = 0; // clear flags on init
+       node->arcs = NULL;
+       node->degree = 0;
+       node->weight = weight;
+       node->index = rg->totnodes;
+       VECCOPY(node->p, eve->co);      
+       
+       BLI_addtail(&rg->nodes, node);
+       rg->totnodes++;
+       
+       return node;
+}
+
+ReebEdge * createArc(ReebGraph *rg, ReebNode *node1, ReebNode *node2)
+{
+       ReebEdge *edge;
+       
+       edge = BLI_edgehash_lookup(rg->emap, node1->index, node2->index);
+       
+       // Only add existing edges that haven't been added yet
+       if (edge == NULL)
+       {
+               ReebArc *arc;
+               ReebNode *v1, *v2;
+               float len, offset;
+               int i;
+               
+               arc = MEM_callocN(sizeof(ReebArc), "reeb arc");
+               edge = MEM_callocN(sizeof(ReebEdge), "reeb edge");
+               
+               arc->flags = 0; // clear flags on init
+               
+               if (node1->weight <= node2->weight)
+               {
+                       v1 = node1;     
+                       v2 = node2;     
+               }
+               else
+               {
+                       v1 = node2;     
+                       v2 = node1;     
+               }
+               
+               arc->v1 = v1;
+               arc->v2 = v2;
+               
+               // increase node degree
+               v1->degree++;
+               v2->degree++;
+
+               BLI_edgehash_insert(rg->emap, node1->index, node2->index, edge);
+               
+               edge->arc = arc;
+               edge->nextEdge = NULL;
+               edge->v1 = v1;
+               edge->v2 = v2;
+               
+               BLI_addtail(&rg->arcs, arc);
+               BLI_addtail(&arc->edges, edge);
+               
+               /* adding buckets for embedding */
+               allocArcBuckets(arc);
+               
+               offset = arc->v1->weight;
+               len = arc->v2->weight - arc->v1->weight;
+
+#if 0
+               /* This is the actual embedding filling described in the paper
+                * the problem is that it only works with really dense meshes
+                */
+               if (arc->bcount > 0)
+               {
+                       addVertToBucket(&(arc->buckets[0]), arc->v1->co);
+                       addVertToBucket(&(arc->buckets[arc->bcount - 1]), arc->v2->co);
+               }
+#else
+               for(i = 0; i < arc->bcount; i++)
+               {
+                       float co[3];
+                       float f = (arc->buckets[i].val - offset) / len;
+                       
+                       VecLerpf(co, v1->p, v2->p, f);
+                       addVertToBucket(&(arc->buckets[i]), co);
+               }
+#endif
+
+       }
+       
+       return edge;
+}
+
+void addTriangleToGraph(ReebGraph *rg, ReebNode * n1, ReebNode * n2, ReebNode * n3)
+{
+       ReebEdge *re1, *re2, *re3;
+       ReebEdge *e1, *e2, *e3;
+       float len1, len2, len3;
+       
+       re1 = createArc(rg, n1, n2);
+       re2 = createArc(rg, n2, n3);
+       re3 = createArc(rg, n3, n1);
+       
+       len1 = (float)fabs(n1->weight - n2->weight);
+       len2 = (float)fabs(n2->weight - n3->weight);
+       len3 = (float)fabs(n3->weight - n1->weight);
+       
+       /* The rest of the algorithm assumes that e1 is the longest edge */
+       
+       if (len1 >= len2 && len1 >= len3)
+       {
+               e1 = re1;
+               e2 = re2;
+               e3 = re3;
+       }
+       else if (len2 >= len1 && len2 >= len3)
+       {
+               e1 = re2;
+               e2 = re1;
+               e3 = re3;
+       }
+       else
+       {
+               e1 = re3;
+               e2 = re2;
+               e3 = re1;
+       }
+       
+       /* And e2 is the lowest edge
+        * If e3 is lower than e2, swap them
+        */
+       if (e3->v1->weight < e2->v1->weight)
+       {
+               ReebEdge *etmp = e2;
+               e2 = e3;
+               e3 = etmp;
+       }
+       
+       
+       mergePaths(rg, e1, e2, e3);
+}
+
+ReebGraph * generateReebGraph(EditMesh *em, int subdivisions)
+{
+       ReebGraph *rg;
+       struct DynamicList * dlist;
+       EditVert *eve;
+       EditFace *efa;
+       int index;
+       int totvert;
+       int totfaces;
+       
+#ifdef DEBUG_REEB
+       int countfaces = 0;
+#endif
+       
+       rg = MEM_callocN(sizeof(ReebGraph), "reeb graph");
+       
+       rg->totnodes = 0;
+       rg->emap = BLI_edgehash_new();
+       
+       totvert = BLI_countlist(&em->verts);
+       totfaces = BLI_countlist(&em->faces);
+       
+       renormalizeWeight(em, 1.0f);
+       
+       /* Spread weight to minimize errors */
+       spreadWeight(em);
+
+       renormalizeWeight(em, (float)subdivisions);
+
+       /* Adding vertice */
+       for(index = 0, eve = em->verts.first; eve; index++, eve = eve->next)
+       {
+               eve->hash = index;
+               eve->f2 = 0;
+               eve->tmp.p = addNode(rg, eve, eve->tmp.fp);
+       }
+       
+       /* Temporarely convert node list to dynamic list, for indexed access */
+       dlist = BLI_dlist_from_listbase(&rg->nodes);
+       
+       /* Adding face, edge per edge */
+       for(efa = em->faces.first; efa; efa = efa->next)
+       {
+               ReebNode *n1, *n2, *n3;
+               
+               n1 = (ReebNode*)BLI_dlist_find_link(dlist, efa->v1->hash);
+               n2 = (ReebNode*)BLI_dlist_find_link(dlist, efa->v2->hash);
+               n3 = (ReebNode*)BLI_dlist_find_link(dlist, efa->v3->hash);
+               
+               addTriangleToGraph(rg, n1, n2, n3);
+               
+               if (efa->v4)
+               {
+                       ReebNode *n4 = (ReebNode*)efa->v4->tmp.p;
+                       addTriangleToGraph(rg, n1, n3, n4);
+               }
+
+#ifdef DEBUG_REEB
+               countfaces++;
+               if (countfaces % 100 == 0)
+               {
+                       printf("face %i of %i\n", countfaces, totfaces);
+               }
+#endif
+               
+               
+       }
+       BLI_listbase_from_dlist(dlist, &rg->nodes);
+       
+       removeNormalNodes(rg);
+       
+       return rg;
+}
+
+/***************************************** WEIGHT UTILS **********************************************/
+
+void renormalizeWeight(EditMesh *em, float newmax)
+{
+       EditVert *eve;
+       float minimum, maximum, range;
+       
+       if (em == NULL || BLI_countlist(&em->verts) == 0)
+               return;
+
+       /* First pass, determine maximum and minimum */
+       eve = em->verts.first;
+       minimum = eve->tmp.fp;
+       maximum = eve->tmp.fp;
+       for(eve = em->verts.first; eve; eve = eve->next)
+       {
+               maximum = MAX2(maximum, eve->tmp.fp);
+               minimum = MIN2(minimum, eve->tmp.fp);
+       }
+       
+       range = maximum - minimum;
+
+       /* Normalize weights */
+       for(eve = em->verts.first; eve; eve = eve->next)
+       {
+               eve->tmp.fp = (eve->tmp.fp - minimum) / range * newmax;
+       }
+}
+
+
+int weightFromLoc(EditMesh *em, int axis)
+{
+       EditVert *eve;
+       
+       if (em == NULL || BLI_countlist(&em->verts) == 0 || axis < 0 || axis > 2)
+               return 0;
+
+       /* Copy coordinate in weight */
+       for(eve = em->verts.first; eve; eve = eve->next)
+       {
+               eve->tmp.fp = eve->co[axis];
+       }
+
+       return 1;
+}
+
+static float cotan_weight(float *v1, float *v2, float *v3)
+{
+       float a[3], b[3], c[3], clen;
+
+       VecSubf(a, v2, v1);
+       VecSubf(b, v3, v1);
+       Crossf(c, a, b);
+
+       clen = VecLength(c);
+
+       if (clen == 0.0f)
+               return 0.0f;
+       
+       return Inpf(a, b)/clen;
+}
+
+int weightToHarmonic(EditMesh *em)
+{
+       NLboolean success;
+       EditVert *eve;
+       EditEdge *eed;
+       EditFace *efa;
+       int totvert = 0;
+       int index;
+       int rval;
+       
+       /* Find local extrema */
+       for(eve = em->verts.first; eve; eve = eve->next)
+       {
+               totvert++;
+       }
+
+       /* Solve with openNL */
+       
+       nlNewContext();
+
+       nlSolverParameteri(NL_NB_VARIABLES, totvert);
+
+       nlBegin(NL_SYSTEM);
+       
+       /* Find local extrema */
+       for(index = 0, eve = em->verts.first; eve; index++, eve = eve->next)
+       {
+               EditEdge *eed;
+               int maximum = 1;
+               int minimum = 1;
+               
+               eve->hash = index; /* Assign index to vertex */
+               
+               NextEdgeForVert(NULL, NULL); /* Reset next edge */
+               for(eed = NextEdgeForVert(em, eve); eed && (maximum || minimum); eed = NextEdgeForVert(em, eve))
+               {
+                       EditVert *eve2;
+                       
+                       if (eed->v1 == eve)
+                       {
+                               eve2 = eed->v2;
+                       }
+                       else
+                       {
+                               eve2 = eed->v1;
+                       }
+                       
+                       /* Adjacent vertex is bigger, not a local maximum */
+                       if (eve2->tmp.fp > eve->tmp.fp)
+                       {
+                               maximum = 0;
+                       }
+                       /* Adjacent vertex is smaller, not a local minimum */
+                       else if (eve2->tmp.fp < eve->tmp.fp)
+                       {
+                               minimum = 0;
+                       }
+               }
+               
+               if (maximum || minimum)
+               {
+                       float w = eve->tmp.fp;
+                       eve->f1 = 0;
+                       nlSetVariable(0, index, w);
+                       nlLockVariable(index);
+               }
+               else
+               {
+                       eve->f1 = 1;
+               }
+       }
+       
+       nlBegin(NL_MATRIX);
+
+       /* Zero edge weight */
+       for(eed = em->edges.first; eed; eed = eed->next)
+       {
+               eed->tmp.l = 0;
+       }
+       
+       /* Add faces count to the edge weight */
+       for(efa = em->faces.first; efa; efa = efa->next)
+       {
+               efa->e1->tmp.l++;
+               efa->e2->tmp.l++;
+               efa->e3->tmp.l++;
+       }
+
+       /* Add faces angle to the edge weight */
+       for(efa = em->faces.first; efa; efa = efa->next)
+       {
+               /* Angle opposite e1 */
+               float t1= cotan_weight(efa->v1->co, efa->v2->co, efa->v3->co) / efa->e2->tmp.l;
+               
+               /* Angle opposite e2 */
+               float t2 = cotan_weight(efa->v2->co, efa->v3->co, efa->v1->co) / efa->e3->tmp.l;
+
+               /* Angle opposite e3 */
+               float t3 = cotan_weight(efa->v3->co, efa->v1->co, efa->v2->co) / efa->e1->tmp.l;
+               
+               int i1 = efa->v1->hash;
+               int i2 = efa->v2->hash;
+               int i3 = efa->v3->hash;
+               
+               nlMatrixAdd(i1, i1, t2+t3);
+               nlMatrixAdd(i2, i2, t1+t3);
+               nlMatrixAdd(i3, i3, t1+t2);
+       
+               nlMatrixAdd(i1, i2, -t3);
+               nlMatrixAdd(i2, i1, -t3);
+       
+               nlMatrixAdd(i2, i3, -t1);
+               nlMatrixAdd(i3, i2, -t1);
+       
+               nlMatrixAdd(i3, i1, -t2);
+               nlMatrixAdd(i1, i3, -t2);
+       }
+       
+       nlEnd(NL_MATRIX);
+
+       nlEnd(NL_SYSTEM);
+
+       success = nlSolveAdvanced(NULL, NL_TRUE);
+
+       if (success)
+       {
+               rval = 1;
+               for(index = 0, eve = em->verts.first; eve; index++, eve = eve->next)
+               {
+                       eve->tmp.fp = nlGetVariable(0, index);
+               }
+       }
+       else
+       {
+               rval = 0;
+       }
+
+       nlDeleteContext(nlGetCurrent());
+
+       return rval;
+}
+
+
+EditEdge * NextEdgeForVert(EditMesh *em, EditVert *v)
+{
+       static EditEdge *e = NULL;
+       
+       /* Reset method, call with NULL mesh pointer */
+       if (em == NULL)
+       {
+               e = NULL;
+               return NULL;
+       }
+       
+       /* first pass, start at the head of the list */
+       if (e == NULL)
+       {
+               e = em->edges.first;
+       }
+       /* subsequent passes, start on the next edge */
+       else
+       {
+               e = e->next;
+       }
+
+       for( ; e ; e = e->next)
+       {
+               if (e->v1 == v || e->v2 == v)
+               {
+                       break;
+               }
+       }       
+       
+       return e;
+}
+
+int weightFromDistance(EditMesh *em)
+{
+       EditVert *eve;
+       int totedge = 0;
+       int vCount = 0;
+       
+       if (em == NULL || BLI_countlist(&em->verts) == 0)
+       {
+               return 0;
+       }
+       
+       totedge = BLI_countlist(&em->edges);
+       
+       if (totedge == 0)
+       {
+               return 0;
+       }
+
+       /* Initialize vertice flags and find at least one selected vertex */
+       for(eve = em->verts.first; eve && vCount == 0; eve = eve->next)
+       {
+               eve->f1 = 0;
+               if (eve->f & SELECT)
+               {
+                       vCount = 1;
+               }
+       }
+       
+       if (vCount == 0)
+       {
+               return 0; /* no selected vert, failure */
+       }
+       else
+       {
+               EditVert *eve, *current_eve = NULL;
+               /* Apply dijkstra spf for each selected vert */
+               for(eve = em->verts.first; eve; eve = eve->next)
+               {
+                       if (eve->f & SELECT)
+                       {
+                               current_eve = eve;
+                               eve->f1 = 1;
+                               
+                               {
+                                       EditEdge *eed = NULL;
+                                       EditEdge *select_eed = NULL;
+                                       EditEdge **edges = NULL;
+                                       float    currentWeight = 0;
+                                       int      eIndex = 0;
+                                       
+                                       edges = MEM_callocN(totedge * sizeof(EditEdge*), "Edges");
+                                       
+                                       /* Calculate edge weight and initialize edge flags */
+                                       for(eed= em->edges.first; eed; eed= eed->next)
+                                       {
+                                               eed->tmp.fp = VecLenf(eed->v1->co, eed->v2->co);
+                                               eed->f1 = 0;
+                                       }
+                                       
+                                       do {
+                                               int i;
+                                               
+                                               current_eve->f1 = 1; /* mark vertex as selected */
+                                               
+                                               /* Add all new edges connected to current_eve to the list */
+                                               NextEdgeForVert(NULL, NULL); // Reset next edge
+                                               for(eed = NextEdgeForVert(em, current_eve); eed; eed = NextEdgeForVert(em, current_eve))
+                                               { 
+                                                       if (eed->f1 == 0)
+                                                       {
+                                                               edges[eIndex] = eed;
+                                                               eed->f1 = 1;
+                                                               eIndex++;
+                                                       }
+                                               }
+                                               
+                                               /* Find next shortest edge */
+                                               select_eed = NULL;
+                                               for(i = 0; i < eIndex; i++)
+                                               {
+                                                       eed = edges[i];
+                                                       
+                                                       if (eed->f1 != 2 && (eed->v1->f1 == 0 || eed->v2->f1 == 0)) /* eed is not selected yet and leads to a new node */
+                                                       {
+                                                               float newWeight = 0;
+                                                               if (eed->v1->f1 == 1)
+                                                               {
+                                                                       newWeight = eed->v1->tmp.fp + eed->tmp.fp;
+                                                               }
+                                                               else
+                                                               {
+                                                                       newWeight = eed->v2->tmp.fp + eed->tmp.fp;
+                                                               }
+                                                               
+                                                               if (select_eed == NULL || newWeight < currentWeight) /* no selected edge or current smaller than selected */
+                                                               {
+                                                                       currentWeight = newWeight;
+                                                                       select_eed = eed;
+                                                               }
+                                                       }
+                                               }
+                                               
+                                               if (select_eed != NULL)
+                                               {
+                                                       select_eed->f1 = 2;
+                                                       
+                                                       if (select_eed->v1->f1 == 0) /* v1 is the new vertex */
+                                                       {
+                                                               current_eve = select_eed->v1;
+                                                       }
+                                                       else /* otherwise, it's v2 */
+                                                       {
+                                                               current_eve = select_eed->v2;
+                                                       }                               
+                                                       current_eve->tmp.fp = currentWeight;
+                                               }
+                                       } while (select_eed != NULL);
+                                       
+                                       MEM_freeN(edges);
+                               }
+                       }
+               }
+       }
+
+       return 1;
+}
+
+MCol MColFromWeight(EditVert *eve)
+{
+       MCol col;
+       col.a = 255;
+       col.b = (char)(eve->tmp.fp * 255);
+       col.g = 0;
+       col.r = (char)((1.0f - eve->tmp.fp) * 255);
+       return col;
+}
+
+void weightToVCol(EditMesh *em)
+{
+       EditFace *efa;
+       MCol *mcol;
+       if (!EM_vertColorCheck()) {
+               return;
+       }
+       
+       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);
+
+               if(efa->v4) {
+                       mcol[3] = MColFromWeight(efa->v4);
+               }
+       }
+}
+
+/****************************************** BUCKET ITERATOR **************************************************/
+
+void initArcIterator(ReebArcIterator *iter, ReebArc *arc, ReebNode *head)
+{
+       iter->arc = arc;
+       
+       if (head == arc->v1)
+       {
+               iter->start = 0;
+               iter->end = arc->bcount - 1;
+               iter->stride = 1;
+       }
+       else
+       {
+               iter->start = arc->bcount - 1;
+               iter->end = 0;
+               iter->stride = -1;
+       }
+       
+       iter->index = iter->start - iter->stride;
+}
+
+void initArcIterator2(ReebArcIterator *iter, ReebArc *arc, int start, int end)
+{
+       iter->arc = arc;
+       
+       iter->start = start;
+       iter->end = end;
+       
+       if (end > start)
+       {
+               iter->stride = 1;
+       }
+       else
+       {
+               iter->stride = -1;
+       }
+
+       iter->index = iter->start - iter->stride;
+}
+
+EmbedBucket * nextBucket(ReebArcIterator *iter)
+{
+       EmbedBucket *result = NULL;
+       
+       if (iter->index != iter->end)
+       {
+               iter->index += iter->stride;
+               result = &(iter->arc->buckets[iter->index]);
+       }
+       
+       return result;
+}