diff --git a/blender_scraps/barbs.py b/blender_scraps/barbs.py
new file mode 100644
index 0000000..7f6ebd7
--- /dev/null
+++ b/blender_scraps/barbs.py
@@ -0,0 +1,122 @@
+# Hasty conversion from the Rust in prosha/src/examples.rs & Barbs
+
+# The 'main' vertices are fine (easily verified by making barbs()
+# always exit early on the limit check) .  Something is wrong with the
+# barbs.  base_incr and barb_incr seem fine on their own.  but even
+# the first iteration of barb() seems to produce garbage geometry.
+
+# Fairly sure self.sides is wrong.  I've checked it against the Rust
+# repeatedly and my best guess is that the rotation matrices are being
+# constructed differently; I don't know what the Rust code uses (it's
+# an external library) whereas mine does it from quaternions.
+
+import numpy as np
+
+import xform
+
+# Mnemonics:
+X = np.array([1.0, 0.0, 0.0])
+Y = np.array([0.0, 1.0, 0.0])
+Z = np.array([0.0, 0.0, 1.0])
+
+class Barbs(object):
+    def __init__(self):
+        # Incremental transform from each stage to the next:
+        self.base_incr = (xform.Transform().
+                          translate(0, 0, 1).
+                          rotate(Z, 0.15).
+                          rotate(X, 0.1).
+                          scale(0.95))
+        self.barb_incr = (xform.Transform().
+                          translate(0, 0, 0.5).
+                          rotate(Y, -0.2).
+                          scale(0.8))
+        # 'Base' vertices, used throughout:
+        self.base = np.array([
+            [-0.5, -0.5, 0.0],
+            [-0.5,  0.5, 0.0],
+            [ 0.5,  0.5, 0.0],
+            [ 0.5, -0.5, 0.0],
+        ])
+        # self.sides[i] gives transformation from a 'base' layer to
+        # the i'th side (0 to 3):
+        self.sides = [
+            xform.Transform().
+            rotate(Z, -i * np.pi/2).
+            rotate(Y, -np.pi/2).
+            translate(0.5, 0.0, 0.5)
+            for i in range(4)
+        ]
+        # Face & vertex accumulators:
+        self.faces = []
+        # self.faces will be a list of tuples (each one of length 4
+        # and containing indices into self.verts)
+        self.verts = []
+        # self.verts will be a list of np.array of shape (3,) but will
+        # be converted last-minute to tuples. (Why: we need to refer
+        # to prior vertices and arithmetic is much easier from an
+        # array, but Blender eventually needs tuples.)
+
+    def run(self, iters) -> (list, list):
+        # Make seed vertices, use them for 'bottom' face, and recurse:
+        self.verts.extend(self.base)
+        self.faces.append((0, 1, 2, 3))
+        self.main(iters, xform.Transform(), [0,1,2,3])
+        verts = [tuple(v) for v in self.verts]
+        faces = [tuple(f) for f in self.faces]
+        return verts, faces
+
+    def limit_check(self, xform: xform.Transform, iters) -> bool:
+        # Assume all scales are the same (for now)
+        sx,_,_ = xform.get_scale()
+        return sx < 0.005 # or iters <= 0
+
+    def main(self, iters, xform, bound):
+
+        if self.limit_check(xform, iters):
+            # Note opposite winding order
+            verts = [bound[i] for i in [3,2,1,0]]
+            self.faces.append(verts)
+            return
+
+        xform2 = xform.compose(self.base_incr)
+        g = xform2.apply_to(self.base)
+        a0 = len(self.verts)
+        self.verts.extend(g)
+
+        # TODO: Turn this to a cleaner loop?
+        self.main(iters - 1, xform2, [a0, a0 + 1, a0 + 2, a0 + 3])
+        self.barb(iters - 1, xform.compose(self.sides[0]),
+                  [bound[0], bound[1], a0 + 1, a0 + 0])
+        self.barb(iters - 1, xform.compose(self.sides[1]),
+                  [bound[1], bound[2], a0 + 2, a0 + 1])
+        self.barb(iters - 1, xform.compose(self.sides[2]),
+                  [bound[2], bound[3], a0 + 3, a0 + 2])
+        self.barb(iters - 1, xform.compose(self.sides[3]),
+                  [bound[3], bound[0], a0 + 0, a0 + 3])
+
+    def barb(self, iters, xform, bound):
+        # DEBUG: This is set True while testing until I figure out
+        # other problems
+        if True or self.limit_check(xform, iters):
+            # Note opposite winding order
+            verts = [bound[i] for i in [3,2,1,0]]
+            self.faces.append(verts)
+            return
+
+        xform2 = xform.compose(self.barb_incr)
+        g = xform2.apply_to(self.base)
+        offset = len(self.verts)
+        self.verts.extend(g)
+
+        # Connect parallel faces:
+        n = len(self.base)
+        for i, b0 in enumerate(bound):
+            j = (i + 1) % n
+            b1 = bound[j]
+            a0 = offset + i
+            a1 = offset + j
+            self.faces.append([a0, a1, b1, b0])
+
+        it2 = 1 # replace with iters-1 once fixed...
+        self.barb(it2, xform2, [offset, offset + 1, offset + 2, offset + 3])
diff --git a/blender_scraps/xform.py b/blender_scraps/xform.py
index 8729d83..5f63a65 100644
--- a/blender_scraps/xform.py
+++ b/blender_scraps/xform.py
@@ -1,8 +1,8 @@
-import numpy
+import numpy as np
 
 def quat2mat(qw, qx, qy, qz):
     s = 1
-    return numpy.array([
+    return np.array([
         [1-2*s*(qy**2+qz**2), 2*s*(qx*qy-qz*qw), 2*s*(qx*qz+qy*qw), 0],
         [2*s*(qx*qy+qz*qw), 1-2*s*(qx**2+qz**2), 2*s*(qy*qz-qx*qw), 0],
         [2*s*(qx*qz-qy*qw), 2*s*(qy*qz+qx*qw), 1-2*s*(qx**2+qy**2), 0],
@@ -16,15 +16,15 @@ def rotation_quaternion(axis, angle):
     axis -- numpy array of shape (3,), with axis to rotate around
     angle -- angle in radians by which to rotate
     """
-    qc = numpy.cos(angle / 2)
-    qs = numpy.sin(angle / 2)
-    qv = qs * numpy.array(axis)
+    qc = np.cos(angle / 2)
+    qs = np.sin(angle / 2)
+    qv = qs * np.array(axis)
     return (qc, qv[0], qv[1], qv[2])
 
 class Transform(object):
     def __init__(self, mtx=None):
         if mtx is None:
-            self.mtx = numpy.identity(4)
+            self.mtx = np.identity(4)
         else:
             self.mtx = mtx
     def _compose(self, mtx2):
@@ -44,17 +44,20 @@ class Transform(object):
         return self._compose(mtx_identity(*a, **kw))
     def apply_to(self, vs):
         # Homogeneous coords, so append a column of ones. vh is then shape (N,4):
-        vh = numpy.hstack([vs, numpy.ones((vs.shape[0], 1), dtype=vs.dtype)])
+        vh = np.hstack([vs, np.ones((vs.shape[0], 1), dtype=vs.dtype)])
         # As we have row vectors, we're doing basically (A*x)^T=(x^T)*(A^T)
         # hence transposing the matrix, while vectors are already transposed.
         return (vh @ self.mtx.T)[:,0:3]
+    def get_scale(self):
+        norms = np.linalg.norm(self.mtx, axis=0)
+        return norms[:3]
 
 def mtx_scale(sx, sy=None, sz=None):
     if sy is None:
         sy = sx
     if sz is None:
         sz = sx
-    return numpy.array([
+    return np.array([
         [sx, 0,  0,  0],
         [0, sy,  0,  0],
         [0,  0, sz,  0],
@@ -62,7 +65,7 @@ def mtx_scale(sx, sy=None, sz=None):
     ])
 
 def mtx_translate(x, y, z):
-    return numpy.array([
+    return np.array([
         [1, 0, 0, x],
         [0, 1, 0, y],
         [0, 0, 1, z],
@@ -75,10 +78,10 @@ def mtx_rotate(axis, angle):
 
 def mtx_reflect(axis):
     # axis must be norm-1
-    axis = numpy.array(axis)
-    axis = axis / numpy.linalg.norm(axis)
+    axis = np.array(axis)
+    axis = axis / np.linalg.norm(axis)
     a,b,c = axis[0], axis[1], axis[2]
-    return numpy.array([
+    return np.array([
         [1-2*a*a, -2*a*b,   -2*a*c,  0],
         [-2*a*b,  1-2*b*b,  -2*b*c,  0],
         [-2*a*c,  -2*b*c,   1-2*c*c, 0],
@@ -86,4 +89,4 @@ def mtx_reflect(axis):
     ])
 
 def mtx_identity():
-    return numpy.eye(4)
+    return np.eye(4)