Initial commit for Harmonic Skeleton generation.

This is very much a work in progress commit to allow me to work outside of home.

While it does somewhat work, I wouldn't recommend anyone to use it.
This commit is contained in:
Martin Poirier
2007-11-07 00:28:45 +00:00
parent 0de103c1cd
commit ec13425eab
12 changed files with 2343 additions and 16 deletions

168
release/scripts/reeb.py Normal file
View File

@@ -0,0 +1,168 @@
#!BPY
"""
Name: 'Reeb graph import'
Blender: 245
Group: 'Import'
Tooltip: 'Imports a reeb graph saved after skeleton generation'
"""
import Blender
def name(count):
if count == -1:
return ""
else:
return "%05" % count
def importGraph(count):
me = Blender.Mesh.New("graph%s" % name(count))
f = open("test%s.txt" % name(count), "r")
verts = []
edges = []
faces = []
i = 0
first = False
SIZE = 0.3
WITH_NODE = False
def addNode(v, s, verts, faces):
if WITH_NODE:
v1 = [v[0], v[1], v[2] + s]
i1 = len(verts)
verts.append(v1)
v2 = [v[0], v[1] + 0.959 * s, v[2] - 0.283 * s]
i2 = len(verts)
verts.append(v2)
v3 = [v[0] - 0.830 * s, v[1] - 0.479 * s, v[2] - 0.283 * s]
i3 = len(verts)
verts.append(v3)
v4 = [v[0] + 0.830 * s, v[1] - 0.479 * s, v[2] - 0.283 * s]
i4 = len(verts)
verts.append(v4)
faces.append([i1,i2,i3])
faces.append([i1,i3,i4])
faces.append([i2,i3,i4])
faces.append([i1,i2,i4])
return 4
else:
return 0
for line in f:
data = line.strip().split(" ")
if data[0] == "v1":
v = [float(x) for x in data[-3:]]
i += addNode(v, SIZE, verts, faces)
verts.append(v)
i += 1
elif data[0] == "v2":
pass
v = [float(x) for x in data[-3:]]
verts.append(v)
edges.append((i-1, i))
i += 1
i += addNode(v, SIZE, verts, faces)
elif data[0] == "b":
verts.append([float(x) for x in data[-3:]])
edges.append((i-1, i))
i += 1
me.verts.extend(verts)
me.edges.extend(edges)
me.faces.extend(faces)
scn = Blender.Scene.GetCurrent()
ob = scn.objects.new(me, "graph%s" % name(count))
#for i in range(16):
# importGraph(i)
importGraph(-1)
import Blender
def name(count):
if count == -1:
return ""
else:
return "%05" % count
def importGraph(count):
me = Blender.Mesh.New("graph%s" % name(count))
f = open("test%s.txt" % name(count), "r")
verts = []
edges = []
faces = []
i = 0
first = False
SIZE = 0.3
WITH_NODE = False
def addNode(v, s, verts, faces):
if WITH_NODE:
v1 = [v[0], v[1], v[2] + s]
i1 = len(verts)
verts.append(v1)
v2 = [v[0], v[1] + 0.959 * s, v[2] - 0.283 * s]
i2 = len(verts)
verts.append(v2)
v3 = [v[0] - 0.830 * s, v[1] - 0.479 * s, v[2] - 0.283 * s]
i3 = len(verts)
verts.append(v3)
v4 = [v[0] + 0.830 * s, v[1] - 0.479 * s, v[2] - 0.283 * s]
i4 = len(verts)
verts.append(v4)
faces.append([i1,i2,i3])
faces.append([i1,i3,i4])
faces.append([i2,i3,i4])
faces.append([i1,i2,i4])
return 4
else:
return 0
for line in f:
data = line.strip().split(" ")
if data[0] == "v1":
v = [float(x) for x in data[-3:]]
i += addNode(v, SIZE, verts, faces)
verts.append(v)
i += 1
elif data[0] == "v2":
pass
v = [float(x) for x in data[-3:]]
verts.append(v)
edges.append((i-1, i))
i += 1
i += addNode(v, SIZE, verts, faces)
elif data[0] == "b":
verts.append([float(x) for x in data[-3:]])
edges.append((i-1, i))
i += 1
me.verts.extend(verts)
me.edges.extend(edges)
me.faces.extend(faces)
scn = Blender.Scene.GetCurrent()
ob = scn.objects.new(me, "graph%s" % name(count))
#for i in range(16):
# importGraph(i)
importGraph(-1)

View File

@@ -44,7 +44,7 @@ struct ListBase;
struct MemFile;
#define BLENDER_VERSION 245
#define BLENDER_SUBVERSION 7
#define BLENDER_SUBVERSION 8
#define BLENDER_MINVERSION 240
#define BLENDER_MINSUBVERSION 0

View File

@@ -324,9 +324,10 @@ void BLI_sortlist(ListBase *listbase, int (*cmp)(void *, void *))
if (listbase->first != listbase->last)
{
for( previous = listbase->first, current = previous->next; current; previous = current, current = next )
for( previous = listbase->first, current = previous->next; current; current = next )
{
next = current->next;
previous = current->prev;
BLI_remlink(listbase, current);
@@ -335,14 +336,7 @@ void BLI_sortlist(ListBase *listbase, int (*cmp)(void *, void *))
previous = previous->prev;
}
if (previous == NULL)
{
BLI_addhead(listbase, current);
}
else
{
BLI_insertlinkafter(listbase, previous, current);
}
BLI_insertlinkafter(listbase, previous, current);
}
}
}

View File

@@ -6599,6 +6599,7 @@ static void do_versions(FileData *fd, Library *lib, Main *main)
}
}
if(main->versionfile <= 245) {
Scene *sce;
bScreen *sc;
Object *ob;
Image *ima;
@@ -6767,6 +6768,19 @@ static void do_versions(FileData *fd, Library *lib, Main *main)
}
}
if (main->versionfile < 245 || main->subversionfile < 8)
{
/* 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_threshold_angle = 45.0f;
sce->toolsettings->skgen_threshold_length = 1.3f;
sce->toolsettings->skgen_options = SKGEN_FILTER_INTERNAL|SKGEN_FILTER_EXTERNAL|SKGEN_REPOSITION|SKGEN_CUT_LENGTH|SKGEN_CUT_ANGLE;
}
}
}
if ((main->versionfile < 245) || (main->versionfile == 245 && main->subversionfile < 2)) {

View File

@@ -37,6 +37,7 @@ struct Base;
struct Bone;
struct bArmature;
struct ListBase;
struct ReebGraph;
typedef struct EditBone
{
@@ -78,6 +79,8 @@ void add_primitiveArmature(int type);
void apply_rot_armature (struct Object *ob, float mat[3][3]);
void docenter_armature (struct Object *ob, int centermode);
void generateSkeletonFromReebGraph(struct ReebGraph *rg);
void clear_armature(struct Object *ob, char mode);
void delete_armature(void);
@@ -117,8 +120,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);
@@ -140,3 +143,4 @@ void show_all_armature_bones(void);
#endif

View File

@@ -51,5 +51,8 @@ void snap_curs_to_grid(void);
void snap_curs_to_sel(void);
void snap_to_center(void);
void generateSkeleton(void);
#endif /* BSE_EDIT_H */

View File

@@ -413,6 +413,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
@@ -721,3 +723,4 @@ enum {
#endif

View File

@@ -0,0 +1,84 @@
/**
* $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];
} 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;
int bcount;
struct EmbedBucket *buckets;
} ReebArc;
struct EditMesh;
void weightToHarmonic(struct EditMesh *em);
void weightFromDistance(struct EditMesh *em);
void weightFromLoc(struct EditMesh *me, int axis);
void renormalizeWeight(struct EditMesh *em, float newmax);
void filterReebGraph(ReebGraph *rg, float threshold);
struct ReebGraph * generateReebGraph(struct EditMesh *me, int subdivisions);
#endif /*REEB_H_*/

View File

@@ -345,7 +345,7 @@ typedef struct ToolSettings {
short uvcalc_flag;
short pad2;
/* Image Paint (8 byte aligned please!) */
struct ImagePaintSettings imapaint;
@@ -363,7 +363,15 @@ 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_threshold_length;
float skgen_threshold_angle;
short skgen_options;
char pad4[6];
} ToolSettings;
/* Used by all brushes to store their properties, which can be directly set
@@ -665,6 +673,14 @@ 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_REPOSITION 4
#define SKGEN_SYMMETRY 8
#define SKGEN_CUT_LENGTH 16
#define SKGEN_CUT_ANGLE 32
#ifdef __cplusplus
}
#endif

View File

@@ -646,7 +646,6 @@ void do_common_editbuts(unsigned short event) // old name, is a mix of object an
allqueue(REDRAWACTION, 0);
break;
default:
if (G.vd==NULL)
break;
@@ -4342,6 +4341,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! */
@@ -4392,6 +4394,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");
@@ -4415,7 +4418,28 @@ static void verify_vertexgroup_name_func(void *datav, void *data2_unused)
unique_vertexgroup_name((bDeformGroup*)datav, OBACT);
}
static void editing_panel_mesh_skgen(Object *ob, Mesh *me)
{
uiBlock *block;
block= uiNewBlock(&curarea->uiblocks, "editing_panel_mesh_skgen", UI_EMBOSS, UI_HELV, curarea->win);
if(uiNewPanel(curarea, block, "Skeleton Generation", "Editing", 960, 0, 318, 204)==0) return;
uiDefBut(block, BUT, B_GEN_SKELETON, "Generate Skeleton", 1125,160,150,39, 0, 0, 0, 0, 0, "Generate Skeleton from Mesh");
uiBlockBeginAlign(block);
uiDefButS(block, NUM, B_DIFF, "Resolution:", 1125,110,150,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", 1125, 90, 58,19, &G.scene->toolsettings->skgen_options, 0, 0, 0, 0, "Filter internal small arcs from graph");
uiDefButF(block, NUM, B_DIFF, "Thresh:", 1186, 90, 89,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", 1125, 70, 58,19, &G.scene->toolsettings->skgen_options, 0, 0, 0, 0, "Filter external small arcs from graph");
uiDefButF(block, NUM, B_DIFF, "Thresh:", 1186, 70, 89,19, &G.scene->toolsettings->skgen_threshold_external,0.0, 1.0, 10, 0, "Specify the threshold ratio for filtering external arcs");
uiDefButBitS(block, TOG, SKGEN_CUT_LENGTH, B_DIFF, "Cut Length", 1125, 50, 58,19, &G.scene->toolsettings->skgen_options, 0, 0, 0, 0, "Subdivide arcs based on length");
uiDefButF(block, NUM, B_DIFF, "Thresh:", 1186, 50, 89,19, &G.scene->toolsettings->skgen_threshold_length,1.0, 2.0, 10, 0, "Specify the threshold ratio for subdivision");
uiDefButBitS(block, TOG, SKGEN_CUT_ANGLE, B_DIFF, "Cut Angle", 1125, 30, 58,19, &G.scene->toolsettings->skgen_options, 0, 0, 0, 0, "Subdivide arcs based on angle");
uiDefButF(block, NUM, B_DIFF, "Thresh:", 1186, 30, 89,19, &G.scene->toolsettings->skgen_threshold_angle,0.0, 90.0, 10, 0, "Specify the threshold angle in degrees for subdivision");
uiDefButBitS(block, TOG, SKGEN_REPOSITION, B_DIFF, "Reposition", 1125, 10,150,19, &G.scene->toolsettings->skgen_options, 0, 0, 0, 0, "Reposition nodes based on embedding instead of original vertice positions");
uiBlockEndAlign(block);
}
static void editing_panel_mesh_tools1(Object *ob, Mesh *me)
{
@@ -4456,7 +4480,7 @@ static void editing_panel_mesh_tools1(Object *ob, Mesh *me)
#endif
uiBlockEndAlign(block);
uiBlockBeginAlign(block);
uiDefButBitI(block, TOG, G_ALLEDGES, 0, "All Edges", 1125, 22,150,19, &G.f, 0, 0, 0, 0, "Displays all edges in object mode without optimization");
uiDefButBitS(block, TOG, B_MESH_X_MIRROR, B_DIFF, "X-axis mirror",1125,0,150,19, &G.scene->toolsettings->editbutflag, 0, 0, 0, 0, "While using transforms, mirrors the transformation");
@@ -5529,6 +5553,7 @@ void editing_panels()
if(G.obedit) {
editing_panel_mesh_tools(ob, ob->data);
editing_panel_mesh_tools1(ob, ob->data);
editing_panel_mesh_skgen(ob, ob->data);
editing_panel_mesh_uvautocalculation();
if (EM_texFaceCheck())
editing_panel_mesh_texface();

View 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"
@@ -103,6 +104,8 @@
#include "PIL_time.h"
#include "reeb.h" // FIX ME
#include "mydevice.h"
#include "blendef.h"
#include "nla.h"
@@ -3142,3 +3145,219 @@ void transform_armature_mirror_update(void)
}
/*****************************************************************************************************/
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;
}
void generateSkeletonFromReebGraph(ReebGraph *rg)
{
GHash *nodeEndMap = NULL;
ReebArc *arc = 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();
setcursor_space(SPACE_VIEW3D, CURSOR_EDIT);
nodeEndMap = BLI_ghash_new(BLI_ghashutil_ptrhash, BLI_ghashutil_ptrcmp);
for (arc = rg->arcs.first; arc; arc = arc->next)
{
EditBone *lastBone = NULL;
EditBone *firstBone = NULL;
EditBone *parentBone = NULL;
ReebNode *head, *tail;
int added = 0;
// printf("/***************** new arc ******************/\n");
if (arc->v1->degree >= arc->v2->degree)
{
head = arc->v1;
tail = arc->v2;
}
else
{
head = arc->v2;
tail = arc->v1;
}
if ((G.scene->toolsettings->skgen_options & SKGEN_CUT_LENGTH) &&
arcLengthRatio(arc) >= G.scene->toolsettings->skgen_threshold_length)
{
EditBone *child = NULL;
EditBone *parent = NULL;
int same = 0;
int index = 0;
int stride = 1;
float lengthLimit = 1.5f; // use value from UI
// If head is the highest node, invert stride and start index
if (head == arc->v2)
{
stride *= -1;
index = arc->bcount -1;
}
parent = add_editbone("Bone");
VECCOPY(parent->head, head->p);
firstBone = parent; /* set first bone in the chain */
while (index >= 0 && index < arc->bcount)
{
float *vec0 = NULL;
float *vec1 = arc->buckets[index].p;
/* first bucket. Previous is head */
if (index == 0 || index == arc->bcount - 1)
{
vec0 = head->p;
}
/* Previous are valid buckets */
else
{
vec0 = arc->buckets[index - 1].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);
if (isnan(f) == 0)
{
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;
parent = child; // new child is next parent
same = 1; // mark as same
}
else
{
// Next bucket
index += stride;
same = 0; // Reset same
}
}
VECCOPY(parent->tail, tail->p);
lastBone = parent; /* set last bone in the chain */
added = 1;
}
if (added == 0)
{
EditBone *bone;
bone = add_editbone("Bone");
VECCOPY(bone->head, head->p);
VECCOPY(bone->tail, tail->p);
/* set first and last bone, since there's only one */
lastBone = bone;
firstBone = bone;
}
BLI_ghash_insert(nodeEndMap, tail, lastBone);
parentBone = BLI_ghash_lookup(nodeEndMap, head);
if (parentBone != NULL)
{
firstBone->parent = parentBone;
firstBone->flag |= BONE_CONNECTED;
}
}
BLI_ghash_free(nodeEndMap, NULL, NULL);
BIF_undo_push("Generate Skeleton");
}

1797
source/blender/src/reeb.c Normal file
View File

@@ -0,0 +1,1797 @@
/**
* $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_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 "BKE_global.h"
#include "BKE_utildefines.h"
#include "BKE_customdata.h"
#include "blendef.h"
#include "ONL_opennl.h"
#include "reeb.h"
#if 0
#define DPRINT(a) a
#else
#define DPRINT(a)
#endif
/*
* 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;
DPRINT(printf("allocating %i buckets for range %f >> %f\n", arc->bcount, arc->v1->weight, arc->v2->weight));
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);
DPRINT(printf("old: %i new: %i\n", oldBCount, arc->bcount));
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);
DPRINT(printf("%i, %i, %i\n", newOffset, oldOffset, len));
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;
DPRINT(printf("copy edge [%p] -> [%p]\n", edge, newEdge));
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)
{
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));
}
}
}
}
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)
{
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);
}
}
/****************************************** 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
if (arc->v1 == newNode)
{
arc->v2->degree--;
}
else
{
arc->v1->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;
}
}
void filterInternalReebGraph(ReebGraph *rg, float threshold)
{
ReebArc *arc = NULL, *nextArc = NULL;
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);
}
arc = nextArc;
}
}
void filterExternalReebGraph(ReebGraph *rg, float threshold)
{
ReebArc *arc = NULL, *nextArc = NULL;
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)
{
printf("merging\n");
printArc(arc);
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);
}
arc = nextArc;
}
}
/************************************** 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)
{
DPRINT(printf("diff %f (%f, %f, %f)\n", eve->tmp.fp, eve->co[0], eve->co[1], eve->co[2]));
lastWeight = eve->tmp.fp;
}
else
{
DPRINT(printf("same %f (%f, %f, %f)\n", eve->tmp.fp, eve->co[0], eve->co[1], eve->co[2]));
work_needed = 1;
eve->tmp.fp = lastWeight + FLT_EPSILON * 2;
lastWeight = eve->tmp.fp;
}
}
}
MEM_freeN(verts);
}
/******************************************** 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)
{
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);
}
/***************************************** MAIN ALGORITHM **********************************************/
ReebArc * findConnectedArc(ReebGraph *rg, ReebArc *arc, ReebNode *v)
{
ReebArc *nextArc = arc->next;
for(nextArc = rg->arcs.first; nextArc; nextArc = nextArc->next)
// for( ; 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);
DPRINT(printf("arc: %f -> %f\n", arc->v1->weight, arc->v2->weight));
if (nextArc == NULL)
printf("uhm1\n");
DPRINT(printf("nextArc: %f -> %f\n", nextArc->v1->weight, nextArc->v2->weight));
// 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);
DPRINT(printf("arc: %f -> %f\n", arc->v1->weight, arc->v2->weight));
if (nextArc == NULL)
printf("uhm %p\n", arc->v2);
DPRINT(printf("nextArc: %f -> %f\n", nextArc->v1->weight, nextArc->v2->weight));
// Merge arc only if needed
if (arc->v2 == nextArc->v1)
{
mergeConnectedArcs(rg, arc, nextArc);
}
// Otherwise, mark down vert
else
{
arc->v2->degree = 3;
}
}
//DPRINT(printf("merged\n"));
}
}
}
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;
for(edge = arc->edges.first; edge && !edgeEquals(edge, e); edge = edge->next)
{ }
nextEdge = edge->nextEdge;
if (nextEdge != NULL)
{
result = nextEdge->arc;
}
if (result == arc)
{
printf("WTF");
getchar();
exit(1);
}
return result;
}
typedef enum {
MERGE_LOWER,
MERGE_HIGHER,
MERGE_APPEND
} MergeDirection;
void mergeArcEdges(ReebGraph *rg, ReebArc *aDst, ReebArc *aSrc, MergeDirection direction)
{
ReebEdge *e = NULL;
DPRINT(printf("\t\t\t\tmerge arc edges\n"));
if (direction == MERGE_APPEND)
{
for(e = aSrc->edges.first; e; e = e->next)
{
DPRINT(printf("adding edge: %i\n", e->edge->f1));
e->arc = aDst; // Edge is stolen by new arc
}
addlisttolist(&aDst->edges , &aSrc->edges);
}
else
{
DPRINT(printf("number of edges: %i\n", BLI_countlist(&aSrc->edges)));
for(e = aSrc->edges.first; e; e = e->next)
{
ReebEdge *newEdge = copyEdge(e);
DPRINT(printf("adding edge[%p]: %i\n", newEdge, e->edge->f1));
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;
}
if (newEdge->nextEdge == newEdge || e->nextEdge == e)
{
printf("ouch\n");
getchar();
}
}
}
}
// return 1 on full merge
int mergeConnectedArcs(ReebGraph *rg, ReebArc *a0, ReebArc *a1)
{
int result = 0;
ReebNode *removedNode = NULL;
DPRINT(printf("\t\t\tmerge arcs\n"));
DPRINT(printf("a0: %f -> %f\n", a0->v1->weight, a0->v2->weight));
DPRINT(printf("a1: %f -> %f\n", a1->v1->weight, a1->v2->weight));
mergeArcEdges(rg, a0, a1, MERGE_APPEND);
DPRINT(printf("append arcs"));
// Bring a0 to the combine length of both arcs
if (a0->v2 == a1->v1)
{
DPRINT(printf("new v2\n"));
removedNode = a0->v2;
a0->v2 = a1->v2;
}
else if (a0->v1 == a1->v2)
{
DPRINT(printf("new v1\n"));
removedNode = a0->v1;
a0->v1 = a1->v1;
}
else
{
printf("OMGWTFBBQ\n");
}
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;
DPRINT(printf("\t\t\tmerge arcs\n"));
DPRINT(printf("a0: %f -> %f\n", a0->v1->weight, a0->v2->weight));
DPRINT(printf("a1: %f -> %f\n", a1->v1->weight, a1->v2->weight));
// TRIANGLE POINTS DOWN
if (a0->v1->weight == a1->v1->weight) // heads are the same
{
DPRINT(printf("down\n"));
if (a0->v2->weight == a1->v2->weight) // tails also the same, arcs can be totally merge together
{
DPRINT(printf("full merge\n"));
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
{
DPRINT(printf("a0 cut in half\n"));
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
{
DPRINT(printf("a1 cut in half\n"));
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
{
DPRINT(printf("up\n"));
if (a0->v1->weight > a1->v1->weight) // a0->v1->weight is in the middle
{
DPRINT(printf("a1 cut in half\n"));
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
{
DPRINT(printf("a0 cut in half\n"));
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)
}
DPRINT(printf("MERGE RESULT:\n"));
DPRINT(printArc(a0));
if (result == 0)
{
DPRINT(printArc(a1));
}
return result;
}
void glueByMergeSort(ReebGraph *rg, ReebArc *a0, ReebArc *a1, ReebEdge *e0, ReebEdge *e1)
{
int total = 0;
//DPRINT(printf("\t\tglue: %i -> %i\n", e0->f1, e1->f1));
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)
{
DPRINT(printf("next a0\n"));
a0 = nextArcMappedToEdge(a0, e0);
}
else
{
DPRINT(printf("next a1\n"));
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;
DPRINT(printf("\tmerge paths\n"));
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->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;
DPRINT(printf("\tcreate arc for edge %i\n", eed->f1));
arc = MEM_callocN(sizeof(ReebArc), "reeb arc");
edge = MEM_callocN(sizeof(ReebEdge), "reeb edge");
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 (arc->bcount > 0)
// {
// addVertToBucket(&(arc->buckets[0]), arc->v1->co);
// addVertToBucket(&(arc->buckets[arc->bcount - 1]), arc->v2->co);
// }
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);
}
}
else
{
DPRINT(printf("\treuse arc for edge %i\n", eed->f1));
}
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;
int countfaces = 0;
int ecount = 0;
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, 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;
DPRINT(printf("adding face %i\n", countfaces));
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);
}
//exportGraph(rg, ecount);
ecount++;
countfaces++;
if (countfaces % 100 == 0)
printf("face %i of %i\n", countfaces, totfaces);
}
//exportGraph(rg, ecount);
ecount++;
BLI_listbase_from_dlist(dlist, &rg->nodes);
removeNormalNodes(rg);
exportGraph(rg, -1);
ecount++;
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;
}
}
void weightFromLoc(EditMesh *em, int axis)
{
EditVert *eve;
if (em == NULL || BLI_countlist(&em->verts) == 0 || axis < 0 || axis > 2)
return;
// Copy coordinate in weight
for(eve = em->verts.first; eve; eve = eve->next)
{
eve->tmp.fp = eve->co[axis];
}
}
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;
}
void weightToHarmonic(EditMesh *em)
{
NLboolean success;
EditVert *eve;
EditEdge *eed;
EditFace *efa;
int totvert = 0;
int index;
// 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; //minimum?0:1;
eve->f1 = 0;
//printf("extrema [%04i] (%.3f, %.3f, %.3f): %f\n", index, eve->co[0], eve->co[1], eve->co[2], w);
nlSetVariable(0, index, w);
nlLockVariable(index);
// nlRightHandSideAdd(index, w);
}
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 angle 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)
{
for(index = 0, eve = em->verts.first; eve; index++, eve = eve->next)
{
eve->tmp.fp = nlGetVariable(0, index);
//printf("w[%04i] = %f\n", index, eve->tmp.fp);
}
}
else
{
printf("whoops!\n");
}
nlDeleteContext(nlGetCurrent());
}
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;
}
void weightFromDistance(EditMesh *em)
{
EditVert *eve, *current_eve = NULL;
int totedge = 0;
int vCount = 0;
if (em == NULL || BLI_countlist(&em->verts) == 0)
{
return;
}
totedge = BLI_countlist(&em->edges);
if (totedge == 0)
{
return;
}
// Initialize vertice flags and find selected vertex
for(eve = em->verts.first; eve; eve = eve->next)
{
eve->f1 = 0;
if (current_eve == NULL && eve->f & SELECT)
{
current_eve = eve;
eve->f1 = 1;
vCount = 1;
}
}
if (current_eve != NULL)
{
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
// DPRINT(printf("adding new edges\n"));
NextEdgeForVert(NULL, NULL); // Reset next edge
for(eed = NextEdgeForVert(em, current_eve); eed; eed = NextEdgeForVert(em, current_eve))
{
if (eed->f1 == 0)
{
// DPRINT(printf("adding new edge[%i]: EID = %i\n", eIndex, eed->EID));
edges[eIndex] = eed;
eed->f1 = 1;
eIndex++;
}
}
// Find next shortest edge
// printf("finding shortest edge\n");
select_eed = NULL;
for(i = 0; i < eIndex; i++)
{
eed = edges[i];
// printf("edges[%i]: EID = %i\n", i, eed->EID);
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
{
// printf("selecting edge %i\n", i);
currentWeight = newWeight;
select_eed = eed;
}
}
}
// printf("adding edge in graph\n");
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;
//DPRINT(printf("adding vert: weight = %f\n", current_eve->tmp.fp));
}
} while (select_eed != NULL);
MEM_freeN(edges);
}
}
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);
}
}
}
void generateSkeleton(void)
{
EditMesh *em = G.editMesh;
ReebGraph *rg = NULL;
short val;
if (em == NULL)
return;
#if 0
val= pupmenu("Orientation%t|X Axis|Y Axis|Z Axis|Distance|Harmonic");
if(val>0) {
switch(val)
{
case 1:
case 2:
case 3:
weightFromLoc(em, val - 1);
break;
case 4:
weightFromDistance(em);
break;
case 5:
weightFromDistance(em);
weightToHarmonic(em);
break;
}
}
#else
weightFromDistance(em);
weightToHarmonic(em);
#endif
renormalizeWeight(em, 1.0f);
weightToVCol(em);
rg = generateReebGraph(em, G.scene->toolsettings->skgen_resolution);
verifyBuckets(rg);
// Remove arcs without embedding
filterNullReebGraph(rg);
printf("-------------- NULL\n");
verifyBuckets(rg);
if (G.scene->toolsettings->skgen_options & SKGEN_FILTER_EXTERNAL)
{
filterExternalReebGraph(rg, G.scene->toolsettings->skgen_threshold_external * G.scene->toolsettings->skgen_resolution);
}
printf("-------------- EXTERNAL\n");
verifyBuckets(rg);
if (G.scene->toolsettings->skgen_options & SKGEN_FILTER_INTERNAL)
{
filterInternalReebGraph(rg, G.scene->toolsettings->skgen_threshold_internal * G.scene->toolsettings->skgen_resolution);
}
printf("-------------- INTERNAL\n");
verifyBuckets(rg);
if (G.scene->toolsettings->skgen_options & SKGEN_REPOSITION)
{
repositionNodes(rg);
}
printf("-------------- POSITION\n");
verifyBuckets(rg);
printf("--------------\n");
// Filtering might have created degree 2 nodes, so remove them
removeNormalNodes(rg);
verifyBuckets(rg);
buildAdjacencyList(rg);
exportGraph(rg, -1);
generateSkeletonFromReebGraph(rg);
freeGraph(rg);
}