Imported bsphere.c, mostly ifdef'd out for now

3 files changed, 1180 insertions, 1 deletions

diff --git a/source/blender/blenkernel/CMakeLists.txt b/source/blender/blenkernel/CMakeLists.txt
index defcef58463..e549a205325 100644
--- a/source/blender/blenkernel/CMakeLists.txt
+++ b/source/blender/blenkernel/CMakeLists.txt
@@ -81,6 +81,7 @@ set(SRC
 	intern/boids.c
 	intern/booleanops_mesh.c
 	intern/brush.c
+	intern/bsphere.c
 	intern/bullet.c
 	intern/bvhutils.c
 	intern/cdderivedmesh.c
diff --git a/source/blender/blenkernel/intern/bsphere.c b/source/blender/blenkernel/intern/bsphere.c
new file mode 100644
index 00000000000..d8a75f62172
--- /dev/null
+++ b/source/blender/blenkernel/intern/bsphere.c
@@ -0,0 +1,1178 @@
+/*
+ * $Id$
+ *
+ * ***** BEGIN GPL LICENSE BLOCK *****
+ *
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software  Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
+ *
+ * The Original Code is Copyright (C) 2011 by Nicholas Bishop
+ * All rights reserved.
+ *
+ * The Original Code is: all of this file.
+ *
+ * Contributor(s): none yet.
+ *
+ * ***** END GPL LICENSE BLOCK *****
+ */
+
+#include "MEM_guardedalloc.h"
+
+#include "DNA_mesh_types.h"
+#include "DNA_meshdata_types.h"
+#include "DNA_modifier_types.h"
+
+#include "BLI_utildefines.h"
+#include "BLI_heap.h"
+#include "BLI_ghash.h"
+#include "BLI_listbase.h"
+#include "BLI_math.h"
+
+#include "BKE_cdderivedmesh.h"
+#include "BKE_DerivedMesh.h"
+#include "BKE_global.h"
+#include "BKE_library.h"
+#include "BKE_main.h"
+#include "BKE_subsurf.h"
+
+#include <assert.h>
+
+typedef struct Tri {
+	struct Tri *next, *prev;
+	unsigned int v[3];
+	float no[3];
+	float area;
+	int marked;
+} Tri;
+
+typedef struct Edge {
+	struct Edge *next, *prev;
+	unsigned int v[2];
+	Tri *tri[2];
+} Edge;
+
+/* note: a `frame' is four vertices (contiguous within the MVert
+   array), stored simply as the index of the first vertex */
+
+static void create_frame(float frame[4][3], const float co[3], float radius,
+			 const float mat[3][3], float x_offset)
+{
+	float rx[3], ry[3], rz[3];
+	int i;
+
+	mul_v3_v3fl(ry, mat[1], radius);
+	mul_v3_v3fl(rz, mat[2], radius);
+	
+	add_v3_v3v3(frame[3], co, ry);
+	add_v3_v3v3(frame[3], frame[3], rz);
+
+	sub_v3_v3v3(frame[2], co, ry);
+	add_v3_v3v3(frame[2], frame[2], rz);
+
+	sub_v3_v3v3(frame[1], co, ry);
+	sub_v3_v3v3(frame[1], frame[1], rz);
+
+	add_v3_v3v3(frame[0], co, ry);
+	sub_v3_v3v3(frame[0], frame[0], rz);
+
+	mul_v3_v3fl(rx, mat[0], x_offset);
+	for(i = 0; i < 4; i++)
+		add_v3_v3v3(frame[i], frame[i], rx);
+}
+
+static void create_mid_frame(float frame[4][3],
+			     const float co1[3], const float co2[3],
+			     const SkinNode *n1, const SkinNode *n2,
+			     const float mat[3][3])
+{
+	create_frame(frame, co1, (n1->radius + n2->radius) / 2,
+		     mat, len_v3v3(co1, co2) / 2);
+}
+
+static void add_face(MFace *mface, int *totface,
+		     int v1, int v2, int v3, int v4)
+{
+	MFace *f;
+
+	f = &mface[*totface];
+	f->v1 = v1;
+	f->v2 = v2;
+	f->v3 = v3;
+	f->v4 = v4 == -1 ? 0 : v4;
+	if(!v4) {
+		f->v1 = v3;
+		f->v2 = v4;
+		f->v3 = v1;
+		f->v4 = v2;
+	}
+	(*totface)++;
+}
+
+static void connect_frames(MFace *mface, int *totface, int frame1, int frame2)
+{
+	int i;
+
+	for(i = 0; i < 4; i++) {
+		add_face(mface, totface,
+			 frame2 + i,
+			 frame2 + (i+1) % 4,
+			 frame1 + (i+1) % 4,
+			 frame1 + i);
+	}
+}
+
+static Tri *add_triangle(ListBase *tris, MVert *mvert, int v1, int v2, int v3)
+{
+	Tri *t;
+
+	t = MEM_callocN(sizeof(Tri), "add_triangle.tri");
+	t->v[0] = v1;
+	t->v[1] = v2;
+	t->v[2] = v3;
+	BLI_addtail(tris, t);
+
+	normal_tri_v3(t->no, mvert[t->v[0]].co, mvert[t->v[1]].co, mvert[t->v[2]].co);
+	t->area = area_tri_v3(mvert[t->v[0]].co, mvert[t->v[1]].co, mvert[t->v[2]].co);
+
+	return t;
+}
+
+static int point_tri_side(const Tri *t, MVert *mvert, int v)
+{
+	float p[3], d;
+	sub_v3_v3v3(p, mvert[v].co, mvert[t->v[0]].co);
+	d = dot_v3v3(t->no, p);
+	if(d < 0) return -1;
+	else if(d > 0) return 1;
+	else return 0;
+}
+
+static int mark_v_outside_tris(ListBase *tris, MVert *mvert, int v)
+{
+	Tri *t;
+	int outside = 0;
+
+	/* probably there's a much better way to do this */
+	for(t = tris->first; t; t = t->next) {
+		if((t->marked = point_tri_side(t, mvert, v) > 0))
+			outside = 1;
+	}
+	return outside;
+}
+
+static int edge_match(int e1_0, int e1_1, const unsigned int e2[2])
+{
+	/* XXX: maybe isn't necesseary to check both directions? */
+	return (e1_0 == e2[0] && e1_1 == e2[1]) ||
+	       (e1_0 == e2[1] && e1_1 == e2[0]);
+}
+
+/* returns true if the edge (e1, e2) is already in edges; that edge is
+   deleted here as well. if not found just returns 0 */
+static int check_for_dup(ListBase *edges, unsigned int e1, unsigned int e2)
+{
+	Edge *e, *next;
+
+	for(e = edges->first; e; e = next) {
+		next = e->next;
+
+		if(edge_match(e1, e2, e->v)) {
+			/* remove the interior edge */
+			BLI_freelinkN(edges, e);
+			return 1;
+		}
+	}
+
+	return 0;
+}
+
+static void expand_boundary_edges(ListBase *edges, Tri *t)
+{
+	Edge *new;
+	int i;
+
+	/* insert each triangle edge into the boundary list; if any of
+	   its edges are already in there, remove the edge entirely */
+	for(i = 0; i < 3; i++) {
+		if(!check_for_dup(edges, t->v[i], t->v[(i+1)%3])) {
+			new = MEM_callocN(sizeof(Edge), "Edge");
+			new->v[0] = t->v[i];
+			new->v[1] = t->v[(i+1)%3];
+			BLI_addtail(edges, new);
+		}
+	}
+}
+
+static int tri_matches_frame(const Tri *t, int frame)
+{
+	int i, j, match = 0;
+	for(i = 0; i < 3; i++) {
+		for(j = 0; j < 4; j++) {
+			if(t->v[i] == frame + j) {
+				match++;
+				if(match >= 3)
+					return 1;
+			}
+		}
+	}
+	return 0;
+}
+
+static void quad_from_tris(const Edge *e, int ndx[4])
+{
+	int i, j, opp = -1;
+
+	/* find what the second tri has that the first doesn't */
+	for(i = 0; i < 3; i++) {
+		if(e->tri[1]->v[i] != e->tri[0]->v[0] &&
+		   e->tri[1]->v[i] != e->tri[0]->v[1] &&
+		   e->tri[1]->v[i] != e->tri[0]->v[2]) {
+			opp = e->tri[1]->v[i];
+			break;
+		}
+	}
+	assert(opp != -1);
+
+	for(i = 0, j = 0; i < 3; i++, j++) {
+		ndx[j] = e->tri[0]->v[i];
+		/* when the triangle edge cuts across our quad-to-be,
+		   throw in the second triangle's vertex */
+		if((e->tri[0]->v[i] == e->v[0] || e->tri[0]->v[i] == e->v[1]) &&
+		   (e->tri[0]->v[(i+1)%3] == e->v[0] || e->tri[0]->v[(i+1)%3] == e->v[1])) {
+			j++;
+			ndx[j] = opp;
+		}
+	}
+}
+
+static void add_quad_from_tris(MFace *mface, int *totface, const Edge *e)
+{
+	int ndx[4];
+
+	quad_from_tris(e, ndx);
+
+	add_face(mface, totface, ndx[0], ndx[1], ndx[2], ndx[3]);
+}
+
+static GHash *calc_edges(ListBase *tris)
+{
+	GHash *adj;
+	ListBase *edges;
+	Edge *e;
+	Tri *t;
+	int i, e1, e2;
+
+	/* XXX: vertex range might be a little funky? so using a
+	   hash here */
+	adj = BLI_ghash_new(BLI_ghashutil_inthash,
+			    BLI_ghashutil_intcmp,
+			    "calc_edges adj");
+
+	for(t = tris->first; t; t = t->next) {
+		for(i = 0; i < 3; i++) {
+			e1 = t->v[i];
+			e2 = t->v[(i+1)%3];
+			assert(e1 != e2);
+			if(e1 > e2)
+				SWAP(int, e1, e2);
+
+			edges = BLI_ghash_lookup(adj, SET_INT_IN_POINTER(e1));
+			if(!edges) {
+				edges = MEM_callocN(sizeof(ListBase),
+						    "calc_edges ListBase");
+				BLI_ghash_insert(adj, SET_INT_IN_POINTER(e1), edges);
+			}
+
+			/* find the edge in the adjacency list */
+			for(e = edges->first; e; e = e->next) {
+				assert(e->v[0] == e1);
+				if(e->v[1] == e2)
+					break;
+			}
+
+			/* if not found, create the edge */
+			if(!e) {
+				e = MEM_callocN(sizeof(Edge), "calc_edges Edge");
+				e->v[0] = e1;
+				e->v[1] = e2;
+				BLI_addtail(edges, e);
+			}
+
+			/* should never be more than two faces
+			   attached to an edge here */
+			assert(!e->tri[0] || !e->tri[1]);
+
+			if(!e->tri[0])
+				e->tri[0] = t;
+			else if(!e->tri[1])
+				e->tri[1] = t;
+		}
+	}
+
+	return adj;
+}
+
+static int is_quad_symmetric(const MVert *mvert, const Edge *e)
+{
+	int ndx[4];
+	float a[3];
+
+	quad_from_tris(e, ndx);
+
+	copy_v3_v3(a, mvert[ndx[0]].co);
+	a[0] = -a[0];
+
+	if(len_v3v3(a, mvert[ndx[1]].co) < FLT_EPSILON) {
+		copy_v3_v3(a, mvert[ndx[2]].co);
+		a[0] = -a[0];
+		if(len_v3v3(a, mvert[ndx[3]].co) < FLT_EPSILON)
+			return 1;
+	}
+	else if(len_v3v3(a, mvert[ndx[3]].co) < FLT_EPSILON) {
+		copy_v3_v3(a, mvert[ndx[2]].co);
+		a[0] = -a[0];
+		if(len_v3v3(a, mvert[ndx[1]].co) < FLT_EPSILON)
+			return 1;
+	}
+
+	return 0;
+}
+
+static void output_hull(const MVert *mvert, MFace *mface, int *totface, ListBase *tris)
+{
+	Heap *heap;
+	GHash *adj;
+	GHashIterator *iter;
+	ListBase *edges;
+	Edge *e;
+	Tri *t;
+	float score;
+
+	heap = BLI_heap_new();
+	adj = calc_edges(tris);
+
+	/* unmark all triangles */
+	for(t = tris->first; t; t = t->next)
+		t->marked = 0;
+
+	/* build heap */
+	iter = BLI_ghashIterator_new(adj);
+	while(!BLI_ghashIterator_isDone(iter)) {
+		edges = BLI_ghashIterator_getValue(iter);
+		for(e = edges->first; e; e = e->next) {
+			/* only care if the edge is used by more than
+			   one triangle */
+			if(e->tri[0] && e->tri[1]) {
+				score = (e->tri[0]->area + e->tri[1]->area) *
+					dot_v3v3(e->tri[0]->no, e->tri[1]->no);
+
+				/* increase score if the triangles
+				   form a symmetric quad */
+				if(is_quad_symmetric(mvert, e))
+					score *= 10;
+
+				BLI_heap_insert(heap, -score, e);
+			}
+		}
+
+		BLI_ghashIterator_step(iter);
+	}
+	BLI_ghashIterator_free(iter);
+
+	while(!BLI_heap_empty(heap)) {
+		e = BLI_heap_popmin(heap);
+		
+		/* if both triangles still free, outupt as a quad */
+		if(!e->tri[0]->marked && !e->tri[1]->marked) {
+			add_quad_from_tris(mface, totface, e);
+			e->tri[0]->marked = 1;
+			e->tri[1]->marked = 1;
+		}
+	}
+
+	/* free edge list */
+	iter = BLI_ghashIterator_new(adj);
+	while(!BLI_ghashIterator_isDone(iter)) {
+		edges = BLI_ghashIterator_getValue(iter);
+		BLI_freelistN(edges);
+		MEM_freeN(edges);	
+		BLI_ghashIterator_step(iter);
+	}
+	BLI_ghashIterator_free(iter);
+	BLI_ghash_free(adj, NULL, NULL);
+	BLI_heap_free(heap, NULL);
+
+	/* write out any remaining triangles */
+	for(t = tris->first; t; t = t->next) {
+		if(!t->marked)
+			add_face(mface, totface, t->v[0], t->v[1], t->v[2], -1);
+	}
+}	
+
+/* for vertex `v', find which triangles must be deleted to extend the
+   hull; find the boundary edges of that hole so that it can be filled
+   with connections to the new vertex, and update the `tri' list to
+   delete the marked triangles */
+static void add_point(ListBase *tris, MVert *mvert, int v)
+{
+	ListBase edges = {NULL, NULL};
+	Tri *t, *next;
+	Edge *e;
+
+	for(t = tris->first; t; t = next) {
+		next = t->next;
+
+		/* check if triangle is `visible' to v */
+		if(t->marked) {
+			expand_boundary_edges(&edges, t);
+			/* remove the triangle */
+			BLI_freelinkN(tris, t);
+		}
+	}
+
+	/* fill hole boundary with triangles to new point */
+	for(e = edges.first; e; e = e->next)
+		add_triangle(tris, mvert, e->v[0], e->v[1], v);
+	BLI_freelistN(&edges);
+}
+
+static void build_hull(MFace *mface, int *totface,
+		       MVert *mvert, int *frames, int totframe)
+{
+	ListBase tris = {NULL, NULL};
+	Tri *t, *next;
+	int i, j;
+
+	/* use first frame to make initial degenerate pyramid */
+	add_triangle(&tris, mvert, frames[0], frames[0] + 1, frames[0] + 2);
+	add_triangle(&tris, mvert, frames[0] + 1, frames[0], frames[0] + 3);
+	add_triangle(&tris, mvert, frames[0] + 2, frames[0] + 1, frames[0] + 3);
+	add_triangle(&tris, mvert, frames[0], frames[0] + 2, frames[0] + 3);
+
+	for(i = 1; i < totframe; i++) {
+		for(j = 0; j < 4; j++) {
+			int v = frames[i] + j;
+
+			/* ignore point already inside hull */
+			if(mark_v_outside_tris(&tris, mvert, v)) {
+				/* expand hull and delete interior triangles */
+				add_point(&tris, mvert, v);
+			}
+		}
+	}
+
+	/* remove triangles that would fill the original frames */
+	for(t = tris.first; t; t = next) {
+		next = t->next;
+
+		for(i = 0; i < totframe; i++) {
+			if(tri_matches_frame(t, frames[i])) {
+				BLI_freelinkN(&tris, t);
+				break;
+			}
+		}
+	}
+
+	output_hull(mvert, mface, totface, &tris);
+
+	BLI_freelistN(&tris);
+}
+
+/* test: build hull on origdm */
+#if 0
+static DerivedMesh *test_hull(DerivedMesh *origdm)
+{
+	DerivedMesh *dm;
+	MVert *mvert, *origmvert;
+	MFace *mface;
+	int *frames;
+	int totvert = 0, totface = 0, totframe, i;
+
+	/* TODO */
+	totface = 2000;
+	totvert = origdm->getNumVerts(origdm);
+	dm = CDDM_new(totvert, 0, totface);
+
+	mvert = CDDM_get_verts(dm);
+	mface = CDDM_get_faces(dm);
+
+	origmvert = CDDM_get_verts(origdm);
+	for(i = 0; i < totvert; i++)
+		mvert[i] = origmvert[i];
+
+	assert(totvert % 4 == 0);
+	totframe = totvert / 4;
+
+	frames = MEM_callocN(sizeof(int) * totframe, "frames");
+	for(i = 0; i < totframe; i++)
+		frames[i] = i*4;
+
+	totface = 0;
+	build_hull(mface, &totface, mvert, frames, totframe);
+
+	CDDM_calc_edges(dm);
+	//CDDM_calc_normals(dm);
+
+	MEM_freeN(frames);
+
+	return dm;
+}
+#endif
+
+/* TODO */
+static void calc_edge_mat(float mat[3][3], const float a[3], const float b[3])
+{
+	float Z[3] = {0, 0, 1};
+	float dot;
+
+	/* x = edge direction */
+	sub_v3_v3v3(mat[0], b, a);
+	normalize_v3(mat[0]);
+
+	dot = dot_v3v3(mat[0], Z);
+	if(dot > -1 + FLT_EPSILON && dot < 1 - FLT_EPSILON) {
+		/* y = Z cross x */
+		cross_v3_v3v3(mat[1], Z, mat[0]);
+		normalize_v3(mat[1]);
+
+		/* z = x cross y */
+		cross_v3_v3v3(mat[2], mat[0], mat[1]);
+		normalize_v3(mat[2]);
+	}
+	else {
+		mat[1][0] = 1;
+		mat[1][1] = 0;
+		mat[1][2] = 0;
+		mat[2][0] = 0;
+		mat[2][1] = 1;
+		mat[2][2] = 0;
+	}
+}
+
+/* BMesh paper */
+#if 0
+static float calc_R_sub_i_squared(const BSphereNode *n)
+{
+	const float alpha = 1.5;
+	//const float alpha = 1;
+	const float R_sub_i = n->size[0] * alpha;
+	return R_sub_i * R_sub_i;
+}
+
+/* equation (2) */
+static float calc_r_squared(const float v[3], const BSphereNode *n)
+{
+	return len_squared_v3v3(v, n->co);
+}
+
+/* equation (1) */
+static float calc_f_sub_i(const float v[3], const BSphereNode *n)
+{
+	float r_squared;
+	float R_sub_i_squared;
+
+	r_squared = calc_r_squared(v, n);
+	R_sub_i_squared = calc_R_sub_i_squared(n);
+
+	if(r_squared <= R_sub_i_squared)
+		return powf(1.0f - (r_squared / R_sub_i_squared), 2);
+	else
+		return 0;
+}
+
+/* equation (3) */
+static float calc_I(const float x[3], float T, const ListBase *base)
+{
+	BSphereNode *n;
+	float a = -T;
+
+	for(n = base->first; n; n = n->next) {
+		a += calc_I(x, 0, &n->children);
+		a += calc_f_sub_i(x, n);
+	}
+
+	//printf("calc_I: %f, %f\n", -T, a);
+
+	return a;
+}
+
+/* kinda my own guess here */
+static void calc_gradient(float g[3], const float x[3], const ListBase *base)
+{
+	BSphereNode *n;
+
+	for(n = base->first; n; n = n->next) {
+		float R_sub_i_squared = calc_R_sub_i_squared(n);
+		float dist = len_v3v3(x, n->co);
+		float f = calc_f_sub_i(x, n);
+		float d = 1.0f / R_sub_i_squared;
+		float tmp[3];
+
+		/* XXX: other possibilities... */
+		if(f > 0) {
+			sub_v3_v3v3(tmp, x, n->co);
+			mul_v3_fl(tmp, -4*d * (1 - d * dist));
+
+			/* not sure if I'm doing this right; intent is
+			   to flip direction when x is `beneath' the
+			   level set */
+			if(len_v3v3(x, n->co) > n->size[0])
+				negate_v3(tmp);
+
+			add_v3_v3(g, tmp);
+		}
+
+		calc_gradient(g, x, &n->children);
+	}
+}
+
+/* equation (5) */
+static float calc_F(const float x[3], const float K[2], float T, const ListBase *base)
+{
+	float f_of_K;
+	float I_target;
+
+	/* TODO */
+	f_of_K = 1.0f / (1 + fabs(K[0]) + fabs(K[1]));
+	//f_of_K = 1;
+
+	/* TODO */
+	I_target = 0;
+
+	return (calc_I(x, T, base) - I_target) /*TODO: * f_of_K*/;
+}
+
+static float smallest_radius(const ListBase *base)
+{
+	BSphereNode *n;
+	float s = FLT_MAX, t;
+
+	for(n = base->first; n; n = n->next) {
+		if(n->size[0] < s)
+			s = n->size[0];
+		if((t = smallest_radius(&n->children)) < s)
+			s = t;
+	}
+
+	return s;
+}
+
+/* equation (7) */
+static float calc_delta_t(float F_max, float k, const ListBase *base, int print)
+{
+	float min_r_sub_i;
+	float step;
+
+	min_r_sub_i = smallest_radius(base);
+	step = min_r_sub_i / powf(2, k);
+
+	if(print) {
+		printf("min_r: %f, 2^k=%f, step=%f\n", min_r_sub_i, powf(2, k), step);
+	}
+
+	return step / F_max;
+}
+
+/* equation (6) */
+static float evolve_x(float x[3], const float K[2], float T,
+		      int k, const ListBase *base, int i, float F_max)
+{
+	float tmp[3], no[3];
+	float F, dt;
+
+	F = calc_F(x, K, T, base);
+
+	if(F < FLT_EPSILON) {
+		//printf("stability reached for ndx=%d\n", i);
+		return 0;
+	}
+
+	/* TODO ;-) */
+	if(F > F_max)
+		F_max = F;
+
+	dt = calc_delta_t(F_max, k, base, 0);
+
+	dt = F / 2;
+
+	if(i == 0)
+		;//printf("F=%.3f, dt=%.3f\n", F, calc_delta_t(F_max, k, base, 1));
+	
+	zero_v3(no);
+	calc_gradient(no, x, base);
+	normalize_v3(no);
+	negate_v3(no);
+
+	mul_v3_v3fl(tmp, no, F * dt);
+	add_v3_v3(x, tmp);
+
+	return F_max;
+}
+
+static void bsphere_evolve(MVert *mvert, int totvert, const ListBase *base,
+			   float T, int steps, float subdiv_level) 
+{
+	int i, s;
+
+	for(i = 0; i < totvert; i++) {
+		float F_max = 0;
+		/* TODO: for now just doing a constant number of steps */
+		for(s = 0; s < steps; s++) {
+			/* TODO */
+			float K[2] = {0, 0};
+
+			F_max = evolve_x(mvert[i].co, K, T, subdiv_level, base, i, F_max);
+		}
+	}
+}
+#endif
+
+typedef enum {
+	START_CAP = 1,
+	END_CAP = 2,
+	NODE_BISECT = 4,
+	CHILD_EDGE_BISECT = 8,
+	
+	IN_FRAME = 16,
+	OUT_FRAME = 32,
+} FrameFlag;
+
+typedef struct {
+	int totframe;
+	int outbase;
+	FrameFlag flag;
+	float co[0][4][3];
+} Frames;
+
+static Frames *frames_create(FrameFlag flag)
+{
+	Frames *frames;
+	int totframe;
+	totframe = (!!(flag & START_CAP) + !!(flag & END_CAP) +
+		    !!(flag & NODE_BISECT) + !!(flag & CHILD_EDGE_BISECT));
+	assert(totframe);
+	frames = MEM_callocN(sizeof(Frames) + sizeof(float)*totframe*4*3,
+			     "frames_create.frames");
+	frames->totframe = totframe;
+	frames->flag = flag;
+	assert((flag & (END_CAP|CHILD_EDGE_BISECT)) != (END_CAP|CHILD_EDGE_BISECT));
+	return frames;
+}
+
+static int frame_offset(Frames *frames, FrameFlag flag)
+{
+	int n = 0;
+	
+	assert(flag == IN_FRAME || flag == OUT_FRAME ||
+	       !(flag & (IN_FRAME|OUT_FRAME)));
+
+	if(flag == IN_FRAME)
+		return 0;
+	else if(flag == OUT_FRAME)
+		return frames->totframe-1;
+	
+	assert((flag & frames->flag) == flag);
+
+	if(flag == START_CAP) return n;
+	if(frames->flag & START_CAP) n++;
+
+	if(flag == NODE_BISECT) return n;
+	if(frames->flag & NODE_BISECT) n++;
+	
+	if(flag == END_CAP) return n;
+	if(frames->flag & END_CAP) n++;
+
+	if(flag == CHILD_EDGE_BISECT) return n;
+	if(frames->flag & CHILD_EDGE_BISECT) n++;
+
+	assert(!"bad frame offset flag");
+}
+
+static int frame_base(Frames *frames, FrameFlag flag)
+{
+	return frames->outbase + 4*frame_offset(frames, flag);
+}
+
+static void *frame_co(Frames *frames, FrameFlag flag)
+{
+	return frames->co[frame_offset(frames, flag)];
+}
+
+static int frames_totvert(Frames *frames)
+{
+	return frames->totframe*4;
+}
+
+#if 0
+static int bsphere_build_frames(BSphereNode *n, GHash *ghash,
+				float (*parent_mat)[3], FrameFlag flag);
+
+/* builds the `cap' of polygons for the end of a BSphere chain */
+static int bsphere_end_node_frames(BSphereNode *n,
+				   GHash *ghash,
+				   float parent_mat[3][3],
+				   FrameFlag flag)
+{
+	Frames *frames;
+	float rad;
+
+	assert(flag == 0 || flag == START_CAP);
+	frames = frames_create(NODE_BISECT|END_CAP|flag, n, ghash);
+
+	/* TODO */
+	rad = n->size[0];
+
+	if(flag & START_CAP)
+		create_frame(frame_co(frames, START_CAP), n, parent_mat, -rad, rad);
+
+	/* frame to bisect the node */
+	create_frame(frame_co(frames, NODE_BISECT), n, parent_mat, 0, rad);
+
+	/* frame to cap end node */
+	create_frame(frame_co(frames, END_CAP), n, parent_mat, rad, rad);
+
+	return frames_totvert(frames);
+}
+
+static int bsphere_connection_node_frames(BSphereNode *n,
+					  GHash *ghash,
+					  float parent_mat[3][3],
+					  FrameFlag flag)
+{
+	Frames *frames;
+	BSphereNode *child;
+	float childmat[3][3], nodemat[3][3], axis[3], angle, rad, child_rad;
+	int totvert;
+
+	assert(flag == 0 || flag == START_CAP);
+
+	child = n->children.first;
+
+	/* if child is not a branch node, bisect the child edge */
+	if(BLI_countlist(&child->children) <= 1)
+		flag |= CHILD_EDGE_BISECT;
+
+	frames = frames_create(NODE_BISECT|flag, n, ghash);
+	totvert = frames_totvert(frames);
+
+	/* TODO */
+	rad = n->size[0];
+	child_rad = child->size[0];
+
+	/* build matrix to give to child */
+	sub_v3_v3v3(childmat[0], child->co, n->co);
+	normalize_v3(childmat[0]);
+	angle = angle_normalized_v3v3(parent_mat[0], childmat[0]);
+	cross_v3_v3v3(axis, parent_mat[0], childmat[0]);
+	normalize_v3(axis);
+	rotate_normalized_v3_v3v3fl(childmat[1], parent_mat[1], axis, angle);
+	rotate_normalized_v3_v3v3fl(childmat[2], parent_mat[2], axis, angle);
+
+	/* build child */
+	totvert += bsphere_build_frames(child, ghash, childmat, 0);
+
+	/* frame that bisects the node */
+	angle /= 2;
+	copy_v3_v3(nodemat[0], parent_mat[0]);
+	rotate_normalized_v3_v3v3fl(nodemat[1], parent_mat[1], axis, angle);
+	rotate_normalized_v3_v3v3fl(nodemat[2], parent_mat[2], axis, angle);
+	create_frame(frame_co(frames, NODE_BISECT), n, nodemat, 0, rad);
+
+	if(flag & START_CAP)
+		create_frame(frame_co(frames, START_CAP), n, nodemat, -rad, rad);
+
+	if(flag & CHILD_EDGE_BISECT)
+		create_mid_frame(frame_co(frames, CHILD_EDGE_BISECT), n, child, childmat);
+	
+	return totvert;
+}
+
+static int bsphere_branch_node_frames(BSphereNode *n, GHash *ghash)
+{
+	BSphereNode *child;
+	int totvert = 0;
+	float rad;
+
+	rad = n->size[0];
+
+	/* build children */
+	for(child = n->children.first; child; child = child->next) {
+		float childmat[3][3];
+
+		calc_edge_mat(childmat, n->co, child->co);
+
+		totvert += bsphere_build_frames(child, ghash, childmat, 0);
+	}
+
+	/* TODO: for now, builds no frames of its own */
+
+	return totvert;
+}
+
+static int bsphere_build_frames(BSphereNode *n, GHash *ghash,
+				float parent_mat[3][3], FrameFlag flag)
+{
+	BSphereNode *child;
+
+	child = n->children.first;
+
+	if(!child)
+		return bsphere_end_node_frames(n, ghash, parent_mat, flag);
+	else if(!child->next)
+		return bsphere_connection_node_frames(n, ghash, parent_mat, flag);
+	else
+		return bsphere_branch_node_frames(n, ghash);
+}
+
+static GHash *bsphere_frames(BSphere *bs, int *totvert)
+{
+	GHash *ghash;
+	BSphereNode *n, *child;
+
+	ghash = BLI_ghash_new(BLI_ghashutil_ptrhash, BLI_ghashutil_ptrcmp,
+			      "bsphere_frames.ghash");
+	(*totvert) = 0;
+
+	for(n = bs->rootbase.first; n; n = n->next) {
+		float mat[3][3];
+
+		child = n->children.first;
+		if(!child || (child && child->next)) {
+			set_v3(mat[0], 0, 0, 1);
+			set_v3(mat[1], 1, 0, 0);
+			set_v3(mat[2], 0, 1, 0);
+		}
+		else if(!child->next) {
+			calc_edge_mat(mat, n->co, child->co);
+		}
+
+		(*totvert) += bsphere_build_frames(n, ghash, mat, START_CAP);
+	}
+
+	return ghash;
+}
+
+static void output_frames(Frames *frames, MVert *mvert, MFace *mface,
+			  int *totvert, int *totface)
+{
+	int b, i, j;
+
+	frames->outbase = (*totvert);
+	for(i = 0; i < frames->totframe; i++) {
+		for(j = 0; j < 4; j++) {
+			copy_v3_v3(mvert[(*totvert)].co, frames->co[i][j]);
+			(*totvert)++;
+		}
+
+		if(i > 0) {
+			connect_frames(mface, totface,
+				       frames->outbase + 4*(i-1),
+				       frames->outbase + 4*i);
+		}
+	}
+
+	if(frames->flag & START_CAP) {
+		b = frame_base(frames, START_CAP);
+		add_face(mface, totface, b, b+1, b+2, b+3);
+	}
+	if(frames->flag & END_CAP) {
+		b = frame_base(frames, END_CAP);
+		add_face(mface, totface, b+3, b+2, b+1, b+0);
+	}
+}
+
+static void bsphere_build_skin(BSphereNode *parent, ListBase *base,
+			       GHash *frames, MVert *mvert,
+			       MFace *mface, int *totvert, int *totface)
+{
+	BSphereNode *n, *child;
+	Frames *node_frames, *child_frames, *parent_frames;
+	int *hull_frames, totchild, i;
+
+	for(n = base->first; n; n = n->next) {
+		if((node_frames = BLI_ghash_lookup(frames, n)))
+			output_frames(node_frames, mvert, mface, totvert, totface);
+
+		bsphere_build_skin(n, &n->children, frames, mvert, mface, totvert, totface);
+
+		child = n->children.first;
+		if(child && !child->next) {
+			child_frames = BLI_ghash_lookup(frames, child);
+			if(child_frames) {
+				connect_frames(mface, totface,
+					       frame_base(node_frames, OUT_FRAME),
+					       frame_base(child_frames, IN_FRAME));
+			}
+		}
+		else if(child && child->next) {
+			totchild = BLI_countlist(&n->children);
+			hull_frames = MEM_callocN(sizeof(int) * (1+totchild),
+						  "bsphere_build_skin.hull_frames");
+			i = 0;
+			for(child = n->children.first; child; child = child->next) {
+				child_frames = BLI_ghash_lookup(frames, child);
+				assert(child_frames);
+				hull_frames[i++] = frame_base(child_frames, IN_FRAME);
+			}
+			parent_frames = BLI_ghash_lookup(frames, parent);
+			if(parent_frames)
+				hull_frames[i] = frame_base(parent_frames, OUT_FRAME);
+			else
+				totchild--;
+			build_hull(mface, totface, mvert, hull_frames, totchild+1);
+			MEM_freeN(hull_frames);
+		}
+	}
+}
+
+static DerivedMesh *bsphere_base_skin(BSphere *bs)
+{
+	DerivedMesh *dm;
+	GHash *ghash;
+	MVert *mvert;
+	MFace *mface;
+	int totvert, totface = 0;
+
+	/* no nodes, nothing to do */
+	if(!bs->rootbase.first)
+		return NULL;
+
+	/* could probably do all this in one fell [recursive] swoop
+	   using some tricksies and clevers, but multiple passes is
+	   easier :) */
+	ghash = bsphere_frames(bs, &totvert);
+
+	totface = totvert * 1.5;
+	/* XXX: sizes OK? */
+	dm = CDDM_new(totvert, 0, totface);
+
+	mvert = CDDM_get_verts(dm);
+	mface = CDDM_get_faces(dm);
+
+	totvert = totface = 0;
+	bsphere_build_skin(NULL, &bs->rootbase, ghash, mvert, mface, &totvert, &totface);
+	BLI_ghash_free(ghash, NULL, (void*)MEM_freeN);
+
+	CDDM_lower_num_verts(dm, totvert);
+	CDDM_lower_num_faces(dm, totface);
+
+	CDDM_calc_edges(dm);
+	CDDM_calc_normals(dm);
+
+	return dm;
+}
+
+/* generate random points (just for testing) */
+DerivedMesh *bsphere_random_points(BSphere *UNUSED(bs))
+{
+	DerivedMesh *dm;
+	MVert *mvert;
+	const int totvert = 1;//50000;
+	int i;
+
+	dm = CDDM_new(totvert, 0, 0);
+
+	mvert = CDDM_get_verts(dm);
+	
+	srand(1);
+	for(i = 0; i < totvert; i++) {
+		mvert[i].co[0] = rand() / (float)RAND_MAX - 0.5;
+		mvert[i].co[1] = rand() / (float)RAND_MAX - 0.5;
+		mvert[i].co[2] = rand() / (float)RAND_MAX - 0.5;
+	}
+	
+	return dm;
+}
+
+DerivedMesh *bsphere_test_surface(BSphere *bs, DerivedMesh *input)
+{
+	DerivedMesh *dm;
+	MVert *mvert;
+
+	if(!input)
+		return NULL;
+
+	dm = CDDM_copy(input);
+	mvert = CDDM_get_verts(dm);
+	bsphere_evolve(mvert, input->getNumVerts(input), &bs->rootbase,
+		       bs->skin_threshold, bs->evolve, 2);
+
+	return dm;
+}
+
+DerivedMesh *bsphere_test_gradient(BSphere *bs, DerivedMesh *input)
+{
+	DerivedMesh *dm;
+	MVert *mvert;
+	MEdge *medge;
+	int totvert;
+	int i;
+
+	if(!input)
+		return NULL;
+
+	totvert = input->getNumVerts(input);
+	dm = CDDM_new(totvert*2, totvert, 0);
+
+	mvert = CDDM_get_verts(dm);
+	medge = CDDM_get_edges(dm);
+
+	memcpy(mvert, CDDM_get_verts(input), sizeof(MVert) * totvert);
+	
+	for(i = 0; i < totvert; i++) {
+		const float *b = mvert[i].co;
+		float *g = mvert[totvert+i].co;
+		float K[2] = {0, 0};
+
+		zero_v3(g);
+		calc_gradient(g, b, &bs->rootbase);
+		normalize_v3(g);
+		
+		mul_v3_fl(g, -calc_F(b, K, bs->skin_threshold, &bs->rootbase));
+		add_v3_v3(g, b);
+
+		medge[i].v1 = i;
+		medge[i].v2 = totvert+i;
+	}
+	
+	return dm;
+}
+
+DerivedMesh *bsphere_get_skin_dm(BSphere *bs)
+{
+	DerivedMesh *dm, *subdm;
+
+	dm = bsphere_base_skin(bs);
+
+	/* subdivide once */
+	if(dm && bs->subdiv_level > 0) {
+		SubsurfModifierData smd;
+
+		memset(&smd, 0, sizeof(smd));
+		smd.subdivType = ME_CC_SUBSURF;
+		smd.levels = bs->subdiv_level;
+		subdm = subsurf_make_derived_from_derived(dm, &smd, 0, 0, 0, 0);
+		dm->release(dm);
+		/* for convenience, convert back to cddm (could use ccgdm
+		   later for better performance though? (TODO)) */
+		dm = CDDM_copy(subdm);
+		subdm->release(subdm);
+
+		bsphere_evolve(CDDM_get_verts(dm), dm->getNumVerts(dm), &bs->rootbase,
+			       bs->skin_threshold, bs->evolve, bs->subdiv_level);
+	}
+
+	return dm;
+}
+#endif
diff --git a/source/blender/blenkernel/intern/customdata.c b/source/blender/blenkernel/intern/customdata.c
index 1f18c78d070..b9d6badfbe9 100644
--- a/source/blender/blenkernel/intern/customdata.c
+++ b/source/blender/blenkernel/intern/customdata.c
@@ -811,7 +811,7 @@ static void layerDefault_skin(void *data, int count)
 	int i;
 	
 	for(i = 0; i < count; i++) {
-		n[i].radius = 1;
+		n[i].radius = 0.25;
 	}
 }