diff --git a/examples/generate_2d.py b/examples/generate_2d.py
new file mode 100644
index 00000000..ff9979ac
--- /dev/null
+++ b/examples/generate_2d.py
@@ -0,0 +1,7 @@
+from sdf import *
+import numpy as np
+
+f = rectangle(1) - circle(0.5)
+points = f.generate(workers=1)
+
+print(np.array(points).shape)
diff --git a/examples/show_2d.py b/examples/show_2d.py
new file mode 100644
index 00000000..0ed0e1d2
--- /dev/null
+++ b/examples/show_2d.py
@@ -0,0 +1,6 @@
+from sdf import *
+from matplotlib import pyplot as plt
+
+f = polygon([[0,0], [1,0], [0,1]]).translate([-3.5, -0.5]) | rectangle(1) |    circle(0.5).translate([3, 0])
+
+f.show()
diff --git a/examples/show_3d.py b/examples/show_3d.py
new file mode 100644
index 00000000..3c3a65e6
--- /dev/null
+++ b/examples/show_3d.py
@@ -0,0 +1,11 @@
+from sdf import *
+  
+t = torus(2, 0.5)
+f = union(
+    t.orient(X),
+    t.orient(Y),
+    t.orient(Z),
+    k=0.1,
+)
+
+f.show()
\ No newline at end of file
diff --git a/sdf/d2.py b/sdf/d2.py
index f16a9d42..d0fe7edd 100644
--- a/sdf/d2.py
+++ b/sdf/d2.py
@@ -2,7 +2,7 @@
 import numpy as np
 import operator
 
-from . import dn, d3, ease
+from . import dn, d3, ease, mesh
 
 # Constants
 
@@ -20,6 +20,7 @@
 class SDF2:
     def __init__(self, f):
         self.f = f
+        self.dim = 2
     def __call__(self, p):
         return self.f(p).reshape((-1, 1))
     def __getattr__(self, name):
@@ -36,6 +37,10 @@ def __sub__(self, other):
     def k(self, k=None):
         self._k = k
         return self
+    def generate(self, *args, **kwargs):
+        return mesh.generate(self, *args, **kwargs)
+    def show(self, *args, **kwargs):
+        return mesh.show(self, *args, **kwargs)
 
 def sdf2(f):
     def wrapper(*args, **kwargs):
diff --git a/sdf/d3.py b/sdf/d3.py
index 237c1cd6..74170526 100644
--- a/sdf/d3.py
+++ b/sdf/d3.py
@@ -21,6 +21,7 @@
 class SDF3:
     def __init__(self, f):
         self.f = f
+        self.dim = 3
     def __call__(self, p):
         return self.f(p).reshape((-1, 1))
     def __getattr__(self, name):
@@ -41,6 +42,8 @@ def generate(self, *args, **kwargs):
         return mesh.generate(self, *args, **kwargs)
     def save(self, path, *args, **kwargs):
         return mesh.save(path, self, *args, **kwargs)
+    def show(self, *args, **kwargs):
+        return mesh.show(self, *args, **kwargs)
     def show_slice(self, *args, **kwargs):
         return mesh.show_slice(self, *args, **kwargs)
 
diff --git a/sdf/mesh.py b/sdf/mesh.py
index 8895d23c..68447a7d 100644
--- a/sdf/mesh.py
+++ b/sdf/mesh.py
@@ -16,7 +16,11 @@
 def _marching_cubes(volume, level=0):
     verts, faces, _, _ = measure.marching_cubes(volume, level)
     return verts[faces].reshape((-1, 3))
-
+    
+def _marching_squares(area, level=0):
+    verts = measure.find_contours(area, level)
+    return np.row_stack(verts)
+    
 def _cartesian_product(*arrays):
     la = len(arrays)
     dtype = np.result_type(*arrays)
@@ -26,59 +30,118 @@ def _cartesian_product(*arrays):
     return arr.reshape(-1, la)
 
 def _skip(sdf, job):
-    X, Y, Z = job
-    x0, x1 = X[0], X[-1]
-    y0, y1 = Y[0], Y[-1]
-    z0, z1 = Z[0], Z[-1]
-    x = (x0 + x1) / 2
-    y = (y0 + y1) / 2
-    z = (z0 + z1) / 2
-    r = abs(sdf(np.array([(x, y, z)])).reshape(-1)[0])
-    d = np.linalg.norm(np.array((x-x0, y-y0, z-z0)))
-    if r <= d:
-        return False
-    corners = np.array(list(itertools.product((x0, x1), (y0, y1), (z0, z1))))
-    values = sdf(corners).reshape(-1)
-    same = np.all(values > 0) if values[0] > 0 else np.all(values < 0)
-    return same
+    if sdf.dim == 2:
+        X, Y = job
+        x0, x1 = X[0], X[-1]
+        y0, y1 = Y[0], Y[-1]
+        x = (x0 + x1) / 2
+        y = (y0 + y1) / 2
+        r = abs(sdf(np.array([(x, y)])).reshape(-1)[0])
+        d = np.linalg.norm(np.array((x-x0, y-y0)))
+        if r <= d:
+            return False
+        corners = np.array(list(itertools.product((x0, x1), (y0, y1))))
+        values = sdf(corners).reshape(-1)
+        same = np.all(values > 0) if values[0] > 0 else np.all(values < 0)
+        return same
+    elif sdf.dim == 3:
+        X, Y, Z = job
+        x0, x1 = X[0], X[-1]
+        y0, y1 = Y[0], Y[-1]
+        z0, z1 = Z[0], Z[-1]
+        x = (x0 + x1) / 2
+        y = (y0 + y1) / 2
+        z = (z0 + z1) / 2
+        r = abs(sdf(np.array([(x, y, z)])).reshape(-1)[0])
+        d = np.linalg.norm(np.array((x-x0, y-y0, z-z0)))
+        if r <= d:
+            return False
+        corners = np.array(list(itertools.product((x0, x1), (y0, y1), (z0, z1))))
+        values = sdf(corners).reshape(-1)
+        same = np.all(values > 0) if values[0] > 0 else np.all(values < 0)
+        return same
+    else:
+        raise ValueError('sdf must be an instance of SDF2 or SDF3')
+    
 
 def _worker(sdf, job, sparse):
-    X, Y, Z = job
-    if sparse and _skip(sdf, job):
-        return None
-        # return _debug_triangles(X, Y, Z)
-    P = _cartesian_product(X, Y, Z)
-    volume = sdf(P).reshape((len(X), len(Y), len(Z)))
-    try:
-        points = _marching_cubes(volume)
-    except Exception:
-        return []
-        # return _debug_triangles(X, Y, Z)
-    scale = np.array([X[1] - X[0], Y[1] - Y[0], Z[1] - Z[0]])
-    offset = np.array([X[0], Y[0], Z[0]])
-    return points * scale + offset
-
-def _estimate_bounds(sdf):
-    # TODO: raise exception if bound estimation fails
-    s = 16
-    x0 = y0 = z0 = -1e9
-    x1 = y1 = z1 = 1e9
-    prev = None
-    for i in range(32):
-        X = np.linspace(x0, x1, s)
-        Y = np.linspace(y0, y1, s)
-        Z = np.linspace(z0, z1, s)
-        d = np.array([X[1] - X[0], Y[1] - Y[0], Z[1] - Z[0]])
-        threshold = np.linalg.norm(d) / 2
-        if threshold == prev:
-            break
-        prev = threshold
+    if sdf.dim == 2: # isinstance check causes circuilar import
+        X, Y = job
+        if sparse and _skip(sdf, job):
+            return None
+            # return _debug_triangles(X, Y)
+        P = _cartesian_product(X, Y)
+        area = sdf(P).reshape((len(X), len(Y)))
+        try:
+            points = _marching_squares(area)
+        except Exception:
+            return []
+            # return _debug_triangles(X, Y, Z)
+        scale = np.array([X[1] - X[0], Y[1] - Y[0]])
+        offset = np.array([X[0], Y[0]])
+        return points * scale + offset
+    elif sdf.dim == 3:
+        X, Y, Z = job
+        if sparse and _skip(sdf, job):
+            return None
+            # return _debug_triangles(X, Y, Z)
         P = _cartesian_product(X, Y, Z)
         volume = sdf(P).reshape((len(X), len(Y), len(Z)))
-        where = np.argwhere(np.abs(volume) <= threshold)
-        x1, y1, z1 = (x0, y0, z0) + where.max(axis=0) * d + d / 2
-        x0, y0, z0 = (x0, y0, z0) + where.min(axis=0) * d - d / 2
-    return ((x0, y0, z0), (x1, y1, z1))
+        try:
+            points = _marching_cubes(volume)
+        except Exception:
+            return []
+            # return _debug_triangles(X, Y, Z)
+        scale = np.array([X[1] - X[0], Y[1] - Y[0], Z[1] - Z[0]])
+        offset = np.array([X[0], Y[0], Z[0]])
+        return points * scale + offset
+    else:
+        raise ValueError('sdf must be an instance of SDF2 or SDF3')
+    
+def _estimate_bounds(sdf):
+    # TODO: raise exception if bound estimation fails
+        
+    if sdf.dim == 2:
+        s = 16
+        x0 = y0 = -1e9
+        x1 = y1  = 1e9
+        prev = None
+        for i in range(32):
+            X = np.linspace(x0, x1, s)
+            Y = np.linspace(y0, y1, s)
+            d = np.array([X[1] - X[0], Y[1] - Y[0]])
+            threshold = np.linalg.norm(d) / 2
+            if threshold == prev:
+                break
+            prev = threshold
+            P = _cartesian_product(X, Y)
+            area = sdf(P).reshape((len(X), len(Y)))
+            where = np.argwhere(np.abs(area) <= threshold)
+            x1, y1 = (x0, y0) + where.max(axis=0) * d + d / 2
+            x0, y0 = (x0, y0) + where.min(axis=0) * d - d / 2
+        return ((x0, y0), (x1, y1))
+    elif sdf.dim == 3:
+        s = 16
+        x0 = y0 = z0 = -1e9
+        x1 = y1 = z1 = 1e9
+        prev = None
+        for i in range(32):
+            X = np.linspace(x0, x1, s)
+            Y = np.linspace(y0, y1, s)
+            Z = np.linspace(z0, z1, s)
+            d = np.array([X[1] - X[0], Y[1] - Y[0], Z[1] - Z[0]])
+            threshold = np.linalg.norm(d) / 2
+            if threshold == prev:
+                break
+            prev = threshold
+            P = _cartesian_product(X, Y, Z)
+            volume = sdf(P).reshape((len(X), len(Y), len(Z)))
+            where = np.argwhere(np.abs(volume) <= threshold)
+            x1, y1, z1 = (x0, y0, z0) + where.max(axis=0) * d + d / 2
+            x0, y0, z0 = (x0, y0, z0) + where.min(axis=0) * d - d / 2
+        return ((x0, y0, z0), (x1, y1, z1))
+    else:
+        raise ValueError('sdf must be an instance of SDF2 or SDF3')
 
 def generate(
         sdf,
@@ -87,6 +150,135 @@ def generate(
         verbose=True, sparse=True):
 
     start = time.time()
+    
+    if sdf.dim == 2:
+        if bounds is None:
+            bounds = _estimate_bounds(sdf)
+        (x0, y0), (x1, y1) = bounds
+
+        if step is None and samples is not None:
+            volume = (x1 - x0) * (y1 - y0)
+            step = (volume / samples) ** (1 / 2)
+
+        try:
+            dx, dy = step
+        except TypeError:
+            dx = dy = step
+
+        if verbose:
+            print('min %g, %g' % (x0, y0))
+            print('max %g, %g' % (x1, y1))
+            print('step %g, %g' % (dx, dy))
+
+        X = np.arange(x0, x1, dx)
+        Y = np.arange(y0, y1, dy)
+
+        s = batch_size
+        Xs = [X[i:i+s+1] for i in range(0, len(X), s)]
+        Ys = [Y[i:i+s+1] for i in range(0, len(Y), s)]
+
+        batches = list(itertools.product(Xs, Ys))
+        num_batches = len(batches)
+        num_samples = sum(len(xs) * len(ys) for xs, ys in batches)
+
+        if verbose:
+            print('%d samples in %d batches with %d workers' %
+                (num_samples, num_batches, workers))
+
+        points = []
+        skipped = empty = nonempty = 0
+        bar = progress.Bar(num_batches, enabled=verbose)
+        pool = ThreadPool(workers)
+        f = partial(_worker, sdf, sparse=sparse)
+        for result in pool.imap(f, batches):
+            bar.increment(1)
+            if result is None:
+                skipped += 1
+            elif len(result) == 0:
+                empty += 1
+            else:
+                nonempty += 1
+                points.extend(result)
+        bar.done()
+
+        if verbose:
+            print('%d skipped, %d empty, %d nonempty' % (skipped, empty, nonempty))
+            seconds = time.time() - start
+            print('%d points in %g seconds' % (len(points), seconds))
+
+        return points
+        
+    elif sdf.dim == 3:
+        if bounds is None:
+            bounds = _estimate_bounds(sdf)
+        (x0, y0, z0), (x1, y1, z1) = bounds
+
+        if step is None and samples is not None:
+            volume = (x1 - x0) * (y1 - y0) * (z1 - z0)
+            step = (volume / samples) ** (1 / 3)
+
+        try:
+            dx, dy, dz = step
+        except TypeError:
+            dx = dy = dz = step
+
+        if verbose:
+            print('min %g, %g, %g' % (x0, y0, z0))
+            print('max %g, %g, %g' % (x1, y1, z1))
+            print('step %g, %g, %g' % (dx, dy, dz))
+
+        X = np.arange(x0, x1, dx)
+        Y = np.arange(y0, y1, dy)
+        Z = np.arange(z0, z1, dz)
+
+        s = batch_size
+        Xs = [X[i:i+s+1] for i in range(0, len(X), s)]
+        Ys = [Y[i:i+s+1] for i in range(0, len(Y), s)]
+        Zs = [Z[i:i+s+1] for i in range(0, len(Z), s)]
+
+        batches = list(itertools.product(Xs, Ys, Zs))
+        num_batches = len(batches)
+        num_samples = sum(len(xs) * len(ys) * len(zs)
+            for xs, ys, zs in batches)
+
+        if verbose:
+            print('%d samples in %d batches with %d workers' %
+                (num_samples, num_batches, workers))
+
+        points = []
+        skipped = empty = nonempty = 0
+        bar = progress.Bar(num_batches, enabled=verbose)
+        pool = ThreadPool(workers)
+        f = partial(_worker, sdf, sparse=sparse)
+        for result in pool.imap(f, batches):
+            bar.increment(1)
+            if result is None:
+                skipped += 1
+            elif len(result) == 0:
+                empty += 1
+            else:
+                nonempty += 1
+                points.extend(result)
+        bar.done()
+
+        if verbose:
+            print('%d skipped, %d empty, %d nonempty' % (skipped, empty, nonempty))
+            triangles = len(points) // 3
+            seconds = time.time() - start
+            print('%d triangles in %g seconds' % (triangles, seconds))
+
+        return points
+        
+    else:
+        raise ValueError('sdf must be an instance of SDF2 or SDF3')
+        
+def generate_2d(
+        sdf,
+        step=None, bounds=None, samples=SAMPLES,
+        workers=WORKERS, batch_size=BATCH_SIZE,
+        verbose=True, sparse=True):
+
+    start = time.time()
 
     if bounds is None:
         bounds = _estimate_bounds(sdf)
@@ -148,6 +340,7 @@ def generate(
 
     return points
 
+
 def save(path, *args, **kwargs):
     points = generate(*args, **kwargs)
     if path.lower().endswith('.stl'):
@@ -241,3 +434,29 @@ def show_slice(*args, **kwargs):
     plt.ylabel(axes[1])
     plt.colorbar(im)
     plt.show()
+    
+def show(sdf, *args, **kwargs):
+    if sdf.dim == 2:
+        from matplotlib import pyplot as plt
+        
+        verts = np.array(sdf.generate(*args, **kwargs))
+        
+        plt.plot(verts[:,0], verts[:,1], 'o')
+        plt.show()
+    elif sdf.dim == 3:
+        import vtkplotlib as vpl
+        import stl #numpy-stl
+
+        points = sdf.generate(*args, **kwargs)
+        vertices, cells = np.unique(points, axis=0, return_inverse=True)
+        faces = cells.reshape((-1, 3))
+
+        shape = stl.mesh.Mesh(np.empty(faces.shape[0], dtype=stl.mesh.Mesh.dtype))
+        for i, f in enumerate(faces):
+            for j in range(3):
+                shape.vectors[i][j] = vertices[f[j],:]
+                
+        vpl.mesh_plot(shape)
+        vpl.show()
+    else:
+        raise ValueError('sdf must be an instance of SDF2 or SDF3')