diff --git a/examples/growth2.py b/examples/growth2.py
new file mode 100644
index 0000000..6b5b257
--- /dev/null
+++ b/examples/growth2.py
@@ -0,0 +1,145 @@
+import axi
+import heapq
+import layers
+import random
+
+from collections import defaultdict
+from math import pi, sin, cos, hypot, floor
+from shapely.geometry import LineString
+
+W, H = axi.A3_SIZE
+
+def make_layer():
+    x = layers.Noise(8).add(layers.Constant(0.6)).clamp()
+    x = x.translate(random.random() * 1000, random.random() * 1000)
+    x = x.scale(0.25, 0.25)
+    x = x.power(1.5)
+    # x = x.subtract(layers.Distance(W / 2, H / 2, min(W, H) / 2, 4))
+    return x
+
+class Grid(object):
+    def __init__(self, r):
+        self.r = r
+        self.size = r / 2 ** 0.5
+        self.points = {}
+        self.lines = {}
+
+    def normalize(self, x, y):
+        i = int(floor(x / self.size))
+        j = int(floor(y / self.size))
+        return (i, j)
+
+    def nearby(self, x, y):
+        points = []
+        lines = []
+        i, j = self.normalize(x, y)
+        for p in range(i - 2, i + 3):
+            for q in range(j - 2, j + 3):
+                if (p, q) in self.points:
+                    points.append(self.points[(p, q)])
+                if (p, q) in self.lines:
+                    lines.append(self.lines[(p, q)])
+        return points, lines
+
+    def insert(self, x, y, line=None):
+        points, lines = self.nearby(x, y)
+        for bx, by in points:
+            if hypot(x - bx, y - by) < self.r:
+                return False
+        i, j = self.normalize(x, y)
+        if line:
+            for other in lines:
+                if line.crosses(other):
+                    return False
+            self.lines[(i, j)] = line
+        self.points[(i, j)] = (x, y)
+        return True
+
+    def remove(self, x, y):
+        i, j = self.normalize(x, y)
+        self.points.pop((i, j))
+        self.lines.pop((i, j))
+
+def new_angle(a, d):
+    if d < 0.1:
+        return random.random() * 2 * pi
+    else:
+        return random.gauss(a, pi / 12)
+
+def poisson_disc(layer, x1, y1, x2, y2, r, n):
+    grid = Grid(r)
+    active = []
+    g = 0
+    while len(active) < 1:
+    # for i in range(1):
+        x = x1 + random.random() * (x2 - x1)
+        y = y1 + random.random() * (y2 - y1)
+        score = layer.get(x, y)
+        if score < 0.9:
+            continue
+        # x = (x1 + x2) / 2.0
+        # y = (y1 + y2) / 2.0
+        a = random.random() * 2 * pi
+        if grid.insert(x, y):
+            print(x, y)
+            heapq.heappush(active, (-score, x, y, a, 0, 0, g))
+            g += 1
+    pairs = []
+    while active:
+        ascore, ax, ay, aa, ai, ad, ag = active[0]
+        for i in range(n):
+            a = new_angle(aa, ad)
+            d = random.random() * r + r
+            x = ax + cos(a) * d
+            y = ay + sin(a) * d
+            if x < x1 or y < y1 or x > x2 or y > y2:
+                continue
+            pair = ((ax, ay), (x, y))
+            line = LineString(pair)
+            if not grid.insert(x, y, line):
+                continue
+            score = layer.get(x, y)
+            # if score < 0.25:
+            #     continue
+            if random.random() < 0.75 and random.random() ** 3 > score:
+                heapq.heappop(active)
+                break
+            pairs.append(pair)
+            heapq.heappush(active, (-score, x, y, a, ai + 1, ad + d, ag))
+            break
+        else:
+            heapq.heappop(active)
+    return grid.points.values(), pairs
+
+def make_path(pairs):
+    lookup = defaultdict(list)
+    for parent, child in pairs:
+        lookup[parent].append(child)
+    root = pairs[0][0]
+    path = []
+    stack = []
+    stack.append(root)
+    while stack:
+        point = stack[-1]
+        path.append(point)
+        if not lookup[point]:
+            stack.pop()
+            continue
+        child = lookup[point].pop()
+        stack.append(child)
+    return path
+
+def main():
+    layer = make_layer()
+    layer.save('layer.png', 0, 0, W, H, 50)
+    points, pairs = poisson_disc(layer, 0, 0, W, H, 0.05, 8)
+    path = make_path(pairs)
+    d = axi.Drawing([path])
+    # d = d.rotate_and_scale_to_fit(W, H, step=90)
+    d = d.scale_to_fit(W, H)
+    d.dump('growth.axi')
+    d.render(bounds=(0, 0, W, H)).write_to_png('growth.png')
+    # axi.draw(d)
+
+if __name__ == '__main__':
+    main()
diff --git a/examples/layers.py b/examples/layers.py
new file mode 100644
index 0000000..372ce32
--- /dev/null
+++ b/examples/layers.py
@@ -0,0 +1,143 @@
+# from alpha_shape import alpha_shape
+from math import hypot
+from PIL import Image
+import noise
+
+class Layer(object):
+    def translate(self, x, y):
+        return Translate(self, x, y)
+    def scale(self, x, y):
+        return Scale(self, x, y)
+    def power(self, power):
+        return Power(self, power)
+    def add(self, other):
+        return Add(self, other)
+    def subtract(self, other):
+        return Subtract(self, other)
+    def multiply(self, other):
+        return Multiply(self, other)
+    def threshold(self, threshold):
+        return Threshold(self, threshold)
+    def clamp(self, lo=0, hi=1):
+        return Clamp(self, lo, hi)
+    def normalize(self, lo, hi, new_lo, new_hi):
+        return Normalize(self, lo, hi, new_lo, new_hi)
+    def filter_points(self, points, lo, hi):
+        return [(x, y) for x, y in points if lo <= self.get(x, y) < hi]
+    def alpha_shape(self, points, lo, hi, alpha):
+        points = self.filter_points(points, lo, hi)
+        return alpha_shape(points, alpha)
+    def save(self, path, x1, y1, x2, y2, scale=1, lo=0, hi=1):
+        w = int(round((x2 - x1) * scale))
+        h = int(round((y2 - y1) * scale))
+        data = bytearray(w * h)
+        for y in range(h):
+            for x in range(w):
+                sx = x1 + (x2 - x1) * x / (w - 1)
+                sy = y1 + (y2 - y1) * y / (h - 1)
+                v = (self.get(sx, sy) - lo) / (hi - lo)
+                v = max(0, min(255, int(v * 255)))
+                data[y*w+x] = v
+        # for y in range(y1, y2):
+        #     for x in range(x1, x2):
+        #         v = (self.get(x, y) - lo) / (hi - lo)
+        #         v = max(0, min(255, int(v * 255)))
+        #         data[y*w+x] = v
+        im = Image.frombytes('L', (w, h), bytes(data))
+        im.save(path, 'png')
+
+class Constant(Layer):
+    def __init__(self, value):
+        self.value = value
+    def get(self, x, y):
+        return self.value
+
+class Noise(Layer):
+    def __init__(self, octaves=1):
+        self.octaves = octaves
+    def get(self, x, y):
+        return noise.snoise2(x, y, self.octaves)
+
+class Translate(Layer):
+    def __init__(self, layer, x, y):
+        self.layer = layer
+        self.x = x
+        self.y = y
+    def get(self, x, y):
+        return self.layer.get(self.x + x, self.y + y)
+
+class Scale(Layer):
+    def __init__(self, layer, x, y):
+        self.layer = layer
+        self.x = x
+        self.y = y
+    def get(self, x, y):
+        return self.layer.get(self.x * x, self.y * y)
+
+class Power(Layer):
+    def __init__(self, layer, power):
+        self.layer = layer
+        self.power = power
+    def get(self, x, y):
+        return self.layer.get(x, y) ** self.power
+
+class Add(Layer):
+    def __init__(self, a, b):
+        self.a = a
+        self.b = b
+    def get(self, x, y):
+        return self.a.get(x, y) + self.b.get(x, y)
+
+class Subtract(Layer):
+    def __init__(self, a, b):
+        self.a = a
+        self.b = b
+    def get(self, x, y):
+        return self.a.get(x, y) - self.b.get(x, y)
+
+class Multiply(Layer):
+    def __init__(self, a, b):
+        self.a = a
+        self.b = b
+    def get(self, x, y):
+        return self.a.get(x, y) * self.b.get(x, y)
+
+class Threshold(Layer):
+    def __init__(self, layer, threshold):
+        self.layer = layer
+        self.threshold = threshold
+    def get(self, x, y):
+        return 0 if self.layer.get(x, y) < self.threshold else 1
+
+class Clamp(Layer):
+    def __init__(self, layer, lo=0, hi=1):
+        self.layer = layer
+        self.lo = lo
+        self.hi = hi
+    def get(self, x, y):
+        v = self.layer.get(x, y)
+        v = min(v, self.hi)
+        v = max(v, self.lo)
+        return v
+
+class Normalize(Layer):
+    def __init__(self, layer, lo, hi, new_lo, new_hi):
+        self.layer = layer
+        self.lo = lo
+        self.hi = hi
+        self.new_lo = new_lo
+        self.new_hi = new_hi
+    def get(self, x, y):
+        v = self.layer.get(x, y)
+        p = (v - self.lo) / (self.hi - self.lo)
+        v = self.new_lo + p * (self.new_hi - self.new_lo)
+        return v
+
+class Distance(Layer):
+    def __init__(self, x, y, maximum, gamma=1):
+        self.x = x
+        self.y = y
+        self.maximum = maximum
+        self.gamma = gamma
+    def get(self, x, y):
+        return (hypot(x - self.x, y - self.y) / self.maximum) ** self.gamma