diff --git a/CHANGELOG.md b/CHANGELOG.md
index 2967f1b..5e31cae 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,8 @@
 # Change log
 
+## 0.1.0
+ * Feature: 3-dimensional prime lattices.
+
 ## 0.0.3
  * Bug fix: Connect nodes that are one unit apart within epsilon tolerance. [#12](https://github.com/xenharmonic-devs/ji-lattice/issues/12)
  * Bug fix: Don't create "diagonal" edges even if the nodes are one unit apart. [#13](https://github.com/xenharmonic-devs/ji-lattice/issues/13)
diff --git a/src/__tests__/lattice-3d.spec.ts b/src/__tests__/lattice-3d.spec.ts
new file mode 100644
index 0000000..3616e23
--- /dev/null
+++ b/src/__tests__/lattice-3d.spec.ts
@@ -0,0 +1,64 @@
+import {describe, it, expect} from 'vitest';
+import {Fraction, LOG_PRIMES, toMonzo} from 'xen-dev-utils';
+import {WGP9, primeSphere, spanLattice3D} from '../lattice-3d';
+
+describe('Wilson-Grady-Pakkanen lattice', () => {
+  it('works for a 7-limit box', () => {
+    const monzos: number[][] = [];
+    for (let i = 0; i < 3; ++i) {
+      const fifths = new Fraction(3).pow(i)!.geoMod(2);
+      monzos.push(toMonzo(fifths));
+      monzos.push(toMonzo(fifths.mul(5).geoMod(2)));
+      monzos.push(toMonzo(fifths.mul(7).geoMod(2)));
+    }
+    monzos[0].push(1);
+    const {vertices, edges} = spanLattice3D(monzos, WGP9());
+    expect(vertices).toEqual([
+      {index: 0, x: 0, y: 0, z: 0},
+      {index: 1, x: 0, y: -40, z: 0},
+      {index: 2, x: 0, y: 0, z: 40},
+      {index: 3, x: 40, y: 0, z: 0},
+      {index: 4, x: 40, y: -40, z: 0},
+      {index: 5, x: 40, y: 0, z: 40},
+      {index: 6, x: 80, y: 0, z: 0},
+      {index: 7, x: 80, y: -40, z: 0},
+      {index: 8, x: 80, y: 0, z: 40},
+    ]);
+    expect(edges).toEqual([
+      {x1: 0, y1: 0, z1: 0, x2: 0, y2: -40, z2: 0, type: 'primary'},
+      {x1: 0, y1: 0, z1: 0, x2: 0, y2: 0, z2: 40, type: 'primary'},
+      {x1: 0, y1: 0, z1: 0, x2: 40, y2: 0, z2: 0, type: 'primary'},
+      {x1: 0, y1: -40, z1: 0, x2: 40, y2: -40, z2: 0, type: 'primary'},
+      {x1: 0, y1: 0, z1: 40, x2: 40, y2: 0, z2: 40, type: 'primary'},
+      {x1: 40, y1: 0, z1: 0, x2: 40, y2: -40, z2: 0, type: 'primary'},
+      {x1: 40, y1: 0, z1: 0, x2: 40, y2: 0, z2: 40, type: 'primary'},
+      {x1: 40, y1: 0, z1: 0, x2: 80, y2: 0, z2: 0, type: 'primary'},
+      {x1: 40, y1: -40, z1: 0, x2: 80, y2: -40, z2: 0, type: 'primary'},
+      {x1: 40, y1: 0, z1: 40, x2: 80, y2: 0, z2: 40, type: 'primary'},
+      {x1: 80, y1: 0, z1: 0, x2: 80, y2: -40, z2: 0, type: 'primary'},
+      {x1: 80, y1: 0, z1: 0, x2: 80, y2: 0, z2: 40, type: 'primary'},
+    ]);
+  });
+});
+
+describe('Prime sphere coordinates', () => {
+  it('produces coordinates for the 11-limit', () => {
+    const {horizontalCoordinates, verticalCoordinates, depthwiseCoordinates} =
+      primeSphere(LOG_PRIMES.slice(0, 5));
+    const coords: string[] = [];
+    for (let i = 0; i < 5; ++i) {
+      coords.push(
+        `${horizontalCoordinates[i].toFixed(3)}, ${verticalCoordinates[
+          i
+        ].toFixed(3)}, ${depthwiseCoordinates[i].toFixed(3)}`
+      );
+    }
+    expect(coords).toEqual([
+      '0.000, 0.000, 0.000',
+      '1.861, -0.509, 0.000',
+      '1.437, 0.900, 0.000',
+      '0.647, -0.211, -0.912',
+      '1.968, -0.086, 0.237',
+    ]);
+  });
+});
diff --git a/src/index.ts b/src/index.ts
index b664f11..c5999fc 100644
--- a/src/index.ts
+++ b/src/index.ts
@@ -1,20 +1,7 @@
 import {LOG_PRIMES, dot, mmod, monzosEqual, sub} from 'xen-dev-utils';
-
-// Small radius of tolerance to accept near unit distances between fractional coordinates as edges.
-const EPSILON = 1e-6;
-
-/**
- * The type of an edge connecting two vertices or a gridline.
- *
- * `"primary"`: Prime-wise connection between two vertices.
- *
- * `"custom"`: User-defined connection between two vertices.
- *
- * `"auxiliary"`: Connection where at least one vertex is auxiliary.
- *
- * `"gridline"`: Line extending across the screen.
- */
-export type EdgeType = 'primary' | 'custom' | 'auxiliary' | 'gridline';
+import {EdgeType} from './types';
+import {connect, project, unproject} from './utils';
+export * from './lattice-3d';
 
 /**
  * A vertex of a 2D graph.
@@ -56,12 +43,6 @@ export type Edge = {
   type: EdgeType;
 };
 
-type Connection = {
-  index1: number;
-  index2: number;
-  type: EdgeType;
-};
-
 /**
  * Options for {@link spanLattice}.
  */
@@ -204,178 +185,6 @@ export function mergeEdges(edges: Edge[]) {
   return result;
 }
 
-/**
- * Calculate the taxicab norm / Manhattan distance between two integral vectors.
- * Restrict movement to whole steps for fractional vectors.
- * Has a tolerance for small errors.
- * @param a Prime exponents of a musical interval.
- * @param b Prime exponents of a musical interval.
- * @returns Integer representing the number of "moves" required to reach `b`from `a`. `NaN` if no legal moves exist.
- */
-function taxicabDistance(
-  a: number[],
-  b: number[],
-  tolerance = EPSILON
-): number {
-  if (a.length > b.length) {
-    return taxicabDistance(b, a);
-  }
-  let result = 0;
-  for (let i = 0; i < a.length; ++i) {
-    const distance = Math.abs(a[i] - b[i]);
-    const move = Math.round(distance);
-    if (Math.abs(distance - move) <= tolerance) {
-      result += move;
-    } else {
-      return NaN;
-    }
-  }
-  for (let i = a.length; i < b.length; ++i) {
-    const distance = Math.abs(b[i]);
-    const move = Math.round(distance);
-    if (Math.abs(distance - move) <= tolerance) {
-      result += move;
-    } else {
-      return NaN;
-    }
-  }
-  return result;
-}
-
-/**
- * Connect monzos that are pre-processed to ignore equaves.
- * @param monzos Array of arrays of prime exponents of musical intervals (usually without prime 2).
- * @param maxDistance Maximum taxicab distance to connect.
- * @returns An array of connections and an array of auxillary nodes.
- */
-function connect(monzos: number[][], maxDistance: number) {
-  if (maxDistance > 2) {
-    throw new Error('Only up to max distance = 2 implemented.');
-  }
-
-  const connections: Connection[] = [];
-  const connectingMonzos: number[][] = [];
-
-  if (maxDistance > 1) {
-    for (let i = 0; i < monzos.length; ++i) {
-      for (let j = i + 1; j < monzos.length; ++j) {
-        const distance = taxicabDistance(monzos[i], monzos[j]);
-        if (distance > 1 && distance <= maxDistance) {
-          const len = Math.max(monzos[i].length, monzos[j].length);
-          gapSearch: for (let k = 0; k < len; ++k) {
-            const gap = (monzos[i][k] ?? 0) - (monzos[j][k] ?? 0);
-            if (Math.abs(gap) === 2) {
-              const monzo = [...monzos[j]];
-              for (let l = monzos[j].length; l < monzos[i].length; ++l) {
-                monzo[l] = 0;
-              }
-              monzo[k] += gap / 2;
-              for (const existing of monzos.concat(connectingMonzos)) {
-                if (monzosEqual(monzo, existing)) {
-                  break gapSearch;
-                }
-              }
-              connectingMonzos.push(monzo);
-              break;
-            } else if (Math.abs(gap) === 1) {
-              for (let l = k + 1; l < len; ++l) {
-                const otherGap = (monzos[i][l] ?? 0) - (monzos[j][l] ?? 0);
-                const monzo = [...monzos[j]];
-                for (let m = monzos[j].length; m < monzos[i].length; ++m) {
-                  monzo[m] = 0;
-                }
-                const otherWay = [...monzo];
-                monzo[k] += gap;
-                otherWay[l] += otherGap;
-                let monzoUnique = true;
-                let otherUnique = true;
-                for (const existing of monzos.concat(connectingMonzos)) {
-                  if (monzosEqual(monzo, existing)) {
-                    monzoUnique = false;
-                  }
-                  if (monzosEqual(otherWay, existing)) {
-                    otherUnique = false;
-                  }
-                }
-                if (monzoUnique) {
-                  connectingMonzos.push(monzo);
-                }
-                if (otherUnique) {
-                  connectingMonzos.push(otherWay);
-                }
-                break gapSearch;
-              }
-            }
-          }
-        }
-      }
-    }
-  }
-  if (maxDistance >= 1) {
-    const primaryLength = monzos.length;
-    monzos = monzos.concat(connectingMonzos);
-    for (let i = 0; i < monzos.length; ++i) {
-      for (let j = i + 1; j < monzos.length; ++j) {
-        const distance = taxicabDistance(monzos[i], monzos[j]);
-        if (distance === 1) {
-          connections.push({
-            index1: i,
-            index2: j,
-            type:
-              i < primaryLength && j < primaryLength ? 'primary' : 'auxiliary',
-          });
-        }
-      }
-    }
-  }
-  return {
-    connections,
-    connectingMonzos,
-  };
-}
-
-function project(monzos: number[][], options: LatticeOptions) {
-  const {horizontalCoordinates, verticalCoordinates} = options;
-  const projected = monzos.map(m => [...m]);
-  const limit = Math.max(
-    horizontalCoordinates.length,
-    verticalCoordinates.length
-  );
-  for (const m of projected) {
-    m.length = Math.min(limit, m.length);
-  }
-  for (let i = limit - 1; i >= 0; --i) {
-    if (horizontalCoordinates[i] || verticalCoordinates[i]) {
-      continue;
-    }
-    for (const m of projected) {
-      m.splice(i, 1);
-    }
-  }
-  return projected;
-}
-
-function unproject(monzos: number[][], options: LatticeOptions) {
-  if (!monzos.length) {
-    return [];
-  }
-  const unprojected = monzos.map(m => [...m]);
-  const {horizontalCoordinates, verticalCoordinates} = options;
-  const limit = Math.max(
-    horizontalCoordinates.length,
-    verticalCoordinates.length
-  );
-  for (let i = 0; i < limit; ++i) {
-    if (horizontalCoordinates[i] || verticalCoordinates[i]) {
-      continue;
-    }
-    for (const u of unprojected) {
-      u.splice(i, 0, 0);
-    }
-  }
-  return unprojected;
-}
-
 /**
  * Compute vertices and edges for a 2D graph representing the lattice of a musical scale in just intonation.
  * @param monzos Prime exponents of the musical intervals in the scale.
@@ -386,11 +195,17 @@ export function spanLattice(monzos: number[][], options: LatticeOptions) {
   const {horizontalCoordinates, verticalCoordinates} = options;
   const maxDistance = options.maxDistance ?? 1;
 
-  let projected = project(monzos, options);
+  let projected = project(monzos, [
+    options.horizontalCoordinates,
+    options.verticalCoordinates,
+  ]);
 
   const {connections, connectingMonzos} = connect(projected, maxDistance);
 
-  const unprojected = unproject(connectingMonzos, options);
+  const unprojected = unproject(connectingMonzos, [
+    options.horizontalCoordinates,
+    options.verticalCoordinates,
+  ]);
 
   const vertices: Vertex[] = [];
   let edges: Edge[] = [];
@@ -424,7 +239,10 @@ export function spanLattice(monzos: number[][], options: LatticeOptions) {
 
   if (options.edgeMonzos) {
     projected = projected.concat(connectingMonzos);
-    let ems = project(options.edgeMonzos, options);
+    let ems = project(options.edgeMonzos, [
+      options.horizontalCoordinates,
+      options.verticalCoordinates,
+    ]);
     ems = ems.concat(ems.map(em => em.map(e => -e)));
     for (let i = 0; i < projected.length; ++i) {
       for (let j = i + 1; j < projected.length; ++j) {
diff --git a/src/lattice-3d.ts b/src/lattice-3d.ts
new file mode 100644
index 0000000..bd3377d
--- /dev/null
+++ b/src/lattice-3d.ts
@@ -0,0 +1,354 @@
+import {dot, mmod, monzosEqual, sub} from 'xen-dev-utils';
+import {EdgeType} from './types';
+import {connect, project, unproject} from './utils';
+
+const EPSILON = 1e-6;
+
+/**
+ * A vertex of a 3D graph.
+ */
+export type Vertex3D = {
+  /** Horizontal coordinate. */
+  x: number;
+  /** Vertical coordinate. */
+  y: number;
+  /** Depthwise coordinate. */
+  z: number;
+  /** Index to input array. */
+  index?: number;
+};
+
+/**
+ * An edge connecting two vertices of a 3D graph.
+ */
+export type Edge3D = {
+  /** First horizontal coordinate. */
+  x1: number;
+  /** First vertical coordinate. */
+  y1: number;
+  /** First depthwise coordinate. */
+  z1: number;
+  /** Second horizontal coordinate. */
+  x2: number;
+  /** Second vertical coordinate. */
+  y2: number;
+  /** Second depthwise coordinate. */
+  z2: number;
+  /** Type of connection. */
+  type: EdgeType;
+};
+
+/**
+ * Options for {@link spanLattice3D}.
+ */
+export type LatticeOptions3D = {
+  /** Mapping for prime x-coordinates. */
+  horizontalCoordinates: number[];
+  /** Mapping for prime y-coordinates. */
+  verticalCoordinates: number[];
+  /** Mapping for prime z-coordinates. */
+  depthwiseCoordinates: number[];
+  /** Maximum prime-wise distance for connecting two inputs. */
+  maxDistance?: number;
+  /** Prime-count vectors of connections in addition the the primes. */
+  edgeMonzos?: number[][];
+  /** Flag to merge short edges into a long ones wherever possible. */
+  mergeEdges?: boolean;
+};
+
+// Coordinates are based on SVG so positive y-direction points down.
+// Based on Kraig Grady's coordinate system https://anaphoria.com/wilsontreasure.html
+// Coordinates for prime 2 and third dimension by Lumi Pakkanen.
+// X-coordinates for every prime up to 23.
+const WGP_X = [23, 40, 0, 0, -14, -8, -5, 0, 20];
+// Y-coordinates for every prime up to 23.
+const WGP_Y = [-45, 0, -40, 0, -18, -4, -32, -25, -3];
+// Z-coordinates for every prime up to 23.
+const WGP_Z = [19, 0, 0, 40, 13, 2, 5, 9, 15];
+
+/**
+ * Combine edges that share an endpoint and slope into longer ones.
+ * @param edges Large number of short edges to merge.
+ * @returns Smaller number of long edges.
+ */
+export function mergeEdges3D(edges: Edge3D[]) {
+  // Choose a canonical orientation.
+  const oriented: Edge3D[] = [];
+  for (const edge of edges) {
+    if (
+      edge.x2 < edge.x1 ||
+      (edge.x2 === edge.x1 && edge.y2 < edge.y1) ||
+      (edge.y2 === edge.y1 && edge.z2 < edge.z1)
+    ) {
+      oriented.push({
+        x1: edge.x2,
+        y1: edge.y2,
+        z1: edge.z2,
+        x2: edge.x1,
+        y2: edge.y1,
+        z2: edge.z1,
+        type: edge.type,
+      });
+    } else {
+      oriented.push(edge);
+    }
+  }
+  oriented.sort((a, b) => a.x1 - b.x1 || a.y1 - b.y1 || a.z1 - b.z1);
+  const result: Edge3D[] = [];
+  const spent = new Set<number>();
+  for (let i = 0; i < oriented.length; ++i) {
+    if (spent.has(i)) {
+      continue;
+    }
+    // eslint-disable-next-line prefer-const
+    let {x1, y1, x2, y2, z1, z2, type} = oriented[i];
+    const dx = x2 - x1;
+    const dy = y2 - y1;
+    const dz = z2 - z1;
+    for (let j = i + 1; j < oriented.length; ++j) {
+      const e = oriented[j];
+      if (e.x1 === x2 && e.y1 === y2 && e.z1 === e.z2 && e.type === type) {
+        const dex = e.x2 - e.x1;
+        const dey = e.y2 - e.y1;
+        const dez = e.z2 - e.z1;
+        if (dex * dy === dx * dey && dex * dz === dz * dez) {
+          x2 = e.x2;
+          y2 = e.y2;
+          z2 = e.z2;
+          spent.add(j);
+        }
+      }
+    }
+    result.push({x1, x2, y1, y2, z1, z2, type});
+  }
+  return result;
+}
+
+/**
+ * Compute vertices and edges for a 2D graph representing the lattice of a musical scale in just intonation.
+ * @param monzos Prime exponents of the musical intervals in the scale.
+ * @param options Options for connecting vertices in the graph.
+ * @returns Vertices and edges of the graph.
+ */
+export function spanLattice3D(monzos: number[][], options: LatticeOptions3D) {
+  const {horizontalCoordinates, verticalCoordinates, depthwiseCoordinates} =
+    options;
+  const maxDistance = options.maxDistance ?? 1;
+
+  const coordss = [
+    verticalCoordinates,
+    horizontalCoordinates,
+    depthwiseCoordinates,
+  ];
+  let projected = project(monzos, coordss);
+
+  const {connections, connectingMonzos} = connect(projected, maxDistance);
+
+  const unprojected = unproject(connectingMonzos, coordss);
+
+  const vertices: Vertex3D[] = [];
+  let edges: Edge3D[] = [];
+
+  for (let index = 0; index < monzos.length; ++index) {
+    vertices.push({
+      index,
+      x: dot(monzos[index], horizontalCoordinates),
+      y: dot(monzos[index], verticalCoordinates),
+      z: dot(monzos[index], depthwiseCoordinates),
+    });
+  }
+
+  for (const monzo of unprojected) {
+    vertices.push({
+      index: undefined,
+      x: dot(monzo, horizontalCoordinates),
+      y: dot(monzo, verticalCoordinates),
+      z: dot(monzo, depthwiseCoordinates),
+    });
+  }
+
+  for (const connection of connections) {
+    const {index1, index2, type} = connection;
+    edges.push({
+      x1: vertices[index1].x,
+      y1: vertices[index1].y,
+      z1: vertices[index1].z,
+      x2: vertices[index2].x,
+      y2: vertices[index2].y,
+      z2: vertices[index2].z,
+      type,
+    });
+  }
+
+  if (options.edgeMonzos) {
+    projected = projected.concat(connectingMonzos);
+    let ems = project(options.edgeMonzos, coordss);
+    ems = ems.concat(ems.map(em => em.map(e => -e)));
+    for (let i = 0; i < projected.length; ++i) {
+      for (let j = i + 1; j < projected.length; ++j) {
+        const diff = sub(projected[i], projected[j]);
+        for (const em of ems) {
+          if (monzosEqual(diff, em)) {
+            edges.push({
+              x1: vertices[i].x,
+              y1: vertices[i].y,
+              z1: vertices[i].z,
+              x2: vertices[j].x,
+              y2: vertices[j].y,
+              z2: vertices[j].z,
+              type:
+                i < monzos.length && j < monzos.length ? 'custom' : 'auxiliary',
+            });
+          }
+        }
+      }
+    }
+  }
+
+  if (options.mergeEdges) {
+    edges = mergeEdges3D(edges);
+  }
+
+  return {
+    vertices,
+    edges,
+  };
+}
+
+function allUnique(vector: number[]) {
+  for (let i = 0; i < vector.length; ++i) {
+    for (let j = i + 1; j < vector.length; ++j) {
+      if (vector[i] === vector[j]) {
+        return false;
+      }
+    }
+  }
+  return true;
+}
+
+/**
+ * Calculate rounded logarithms modulo divisions. All forced to be unique.
+ * @param logs Logarithms of primes.
+ * @param divisions Number of divisions of the first prime.
+ * @param searchResolution Resolution for GPV search. Set to 0 to disable (default).
+ * @returns Array of steps for each prime modulo the number of divisions.
+ */
+export function modVal(
+  logs: number[],
+  divisions: number,
+  searchResolution = 0
+) {
+  if (logs.length > divisions) {
+    throw new Error(`Too many logarithms to fit into ${divisions} notes.`);
+  }
+  // Try to find a GPV.
+  for (let i = 0; i < searchResolution; ++i) {
+    const offset = (0.5 * i) / searchResolution;
+    const normalizer = (divisions + offset) / logs[0];
+    const modval = logs.map(l => mmod(Math.round(l * normalizer), divisions));
+    if (allUnique(modval)) {
+      return modval;
+    }
+    if (i) {
+      const normalizer = (divisions - offset) / logs[0];
+      const modval = logs.map(l => mmod(Math.round(l * normalizer), divisions));
+      if (allUnique(modval)) {
+        return modval;
+      }
+    }
+  }
+  // No GPV with unique entries found. Use force.
+  const normalizer = divisions / logs[0];
+  const val = logs.map(l => Math.round(l * normalizer), divisions);
+  for (let i = 1; i < val.length; ++i) {
+    const reserved = new Set<number>();
+    for (let j = 0; j < i; ++j) {
+      reserved.add(mmod(val[j], divisions));
+    }
+    if (!reserved.has(mmod(val[i], divisions))) {
+      continue;
+    }
+    const s = Math.sign(logs[i] * divisions - val[i]);
+    for (let j = 0; 2 * j <= divisions; ++j) {
+      if (!reserved.has(mmod(val[i] + j * s, divisions))) {
+        val[i] += j * s;
+        break;
+      }
+      if (!reserved.has(mmod(val[i] - j * s, divisions))) {
+        val[i] -= j * s;
+        break;
+      }
+    }
+  }
+  return val.map(v => mmod(v, divisions));
+}
+
+/**
+ * Get Wilson-Grady-Pakkanen coordinates for the first 9 primes.
+ * @param equaveIndex Index of the prime to use as the interval of equivalence.
+ * @returns An array of horizontal coordinates for each prime and the same for vertical and depthwise coordinates.
+ */
+export function WGP9(equaveIndex = 0): LatticeOptions3D {
+  const horizontalCoordinates = [...WGP_X];
+  const verticalCoordinates = [...WGP_Y];
+  const depthwiseCoordinates = [...WGP_Z];
+  horizontalCoordinates[equaveIndex] = 0;
+  verticalCoordinates[equaveIndex] = 0;
+  depthwiseCoordinates[equaveIndex] = 0;
+  return {
+    horizontalCoordinates,
+    verticalCoordinates,
+    depthwiseCoordinates,
+  };
+}
+
+/**
+ * Compute coordinates based on sizes of primes that lie on the surface of a sphere offset on the x-axis.
+ * @param logs Logarithms of (formal) primes with the prime of equivalence first.
+ * @param searchResolution Search resolution for optimizing orthogonality of the resulting set.
+ * @returns An array of horizontal coordinates for each prime and the same for vertical and depthwise coordinates.
+ */
+export function primeSphere(logs: number[], searchResolution = 1024) {
+  const dp = (2 * Math.PI) / searchResolution;
+  const horizontalCoordinates: number[] = [];
+  const verticalCoordinates: number[] = [];
+  const depthwiseCoordinates: number[] = [];
+  const dt = (2 * Math.PI) / logs[0];
+  for (const log of logs) {
+    const theta = log * dt;
+    const x = 1 - Math.cos(theta);
+    const u = Math.sin(theta);
+    let y = u;
+    let z = 0;
+    if (horizontalCoordinates.length > 1) {
+      // Find the most orthogonal rotation around the x-axis
+      let bestError = Infinity;
+      for (let j = 0; j < 1024; ++j) {
+        const phi = dp * j;
+        const yc = Math.cos(phi) * u;
+        const zc = Math.sin(phi) * u;
+        let error = 0;
+        for (let i = 0; i < horizontalCoordinates.length; ++i) {
+          error +=
+            (x * horizontalCoordinates[i] +
+              yc * verticalCoordinates[i] +
+              zc * depthwiseCoordinates[i]) **
+            2;
+        }
+        if (error + EPSILON < bestError) {
+          bestError = error;
+          y = yc;
+          z = zc;
+        }
+      }
+    }
+    horizontalCoordinates.push(x);
+    verticalCoordinates.push(y);
+    depthwiseCoordinates.push(z);
+  }
+  return {
+    horizontalCoordinates,
+    verticalCoordinates,
+    depthwiseCoordinates,
+  };
+}
diff --git a/src/types.ts b/src/types.ts
new file mode 100644
index 0000000..9ca8edb
--- /dev/null
+++ b/src/types.ts
@@ -0,0 +1,12 @@
+/**
+ * The type of an edge connecting two vertices or a gridline.
+ *
+ * `"primary"`: Prime-wise connection between two vertices.
+ *
+ * `"custom"`: User-defined connection between two vertices.
+ *
+ * `"auxiliary"`: Connection where at least one vertex is auxiliary.
+ *
+ * `"gridline"`: Line extending across the screen.
+ */
+export type EdgeType = 'primary' | 'custom' | 'auxiliary' | 'gridline';
diff --git a/src/utils.ts b/src/utils.ts
new file mode 100644
index 0000000..7dac00a
--- /dev/null
+++ b/src/utils.ts
@@ -0,0 +1,175 @@
+import {monzosEqual} from 'xen-dev-utils';
+import {type EdgeType} from './types';
+
+// Small radius of tolerance to accept near unit distances between fractional coordinates as edges.
+const EPSILON = 1e-6;
+
+type Connection = {
+  index1: number;
+  index2: number;
+  type: EdgeType;
+};
+
+/**
+ * Calculate the taxicab norm / Manhattan distance between two integral vectors.
+ * Restrict movement to whole steps for fractional vectors.
+ * Has a tolerance for small errors.
+ * @param a Prime exponents of a musical interval.
+ * @param b Prime exponents of a musical interval.
+ * @returns Integer representing the number of "moves" required to reach `b`from `a`. `NaN` if no legal moves exist.
+ */
+function taxicabDistance(
+  a: number[],
+  b: number[],
+  tolerance = EPSILON
+): number {
+  if (a.length > b.length) {
+    return taxicabDistance(b, a);
+  }
+  let result = 0;
+  for (let i = 0; i < a.length; ++i) {
+    const distance = Math.abs(a[i] - b[i]);
+    const move = Math.round(distance);
+    if (Math.abs(distance - move) <= tolerance) {
+      result += move;
+    } else {
+      return NaN;
+    }
+  }
+  for (let i = a.length; i < b.length; ++i) {
+    const distance = Math.abs(b[i]);
+    const move = Math.round(distance);
+    if (Math.abs(distance - move) <= tolerance) {
+      result += move;
+    } else {
+      return NaN;
+    }
+  }
+  return result;
+}
+
+/**
+ * Connect monzos that are pre-processed to ignore equaves.
+ * @param monzos Array of arrays of prime exponents of musical intervals (usually without prime 2).
+ * @param maxDistance Maximum taxicab distance to connect.
+ * @returns An array of connections and an array of auxillary nodes.
+ */
+export function connect(monzos: number[][], maxDistance: number) {
+  if (maxDistance > 2) {
+    throw new Error('Only up to max distance = 2 implemented.');
+  }
+
+  const connections: Connection[] = [];
+  const connectingMonzos: number[][] = [];
+
+  if (maxDistance > 1) {
+    for (let i = 0; i < monzos.length; ++i) {
+      for (let j = i + 1; j < monzos.length; ++j) {
+        const distance = taxicabDistance(monzos[i], monzos[j]);
+        if (distance > 1 && distance <= maxDistance) {
+          const len = Math.max(monzos[i].length, monzos[j].length);
+          gapSearch: for (let k = 0; k < len; ++k) {
+            const gap = (monzos[i][k] ?? 0) - (monzos[j][k] ?? 0);
+            if (Math.abs(gap) === 2) {
+              const monzo = [...monzos[j]];
+              for (let l = monzos[j].length; l < monzos[i].length; ++l) {
+                monzo[l] = 0;
+              }
+              monzo[k] += gap / 2;
+              for (const existing of monzos.concat(connectingMonzos)) {
+                if (monzosEqual(monzo, existing)) {
+                  break gapSearch;
+                }
+              }
+              connectingMonzos.push(monzo);
+              break;
+            } else if (Math.abs(gap) === 1) {
+              for (let l = k + 1; l < len; ++l) {
+                const otherGap = (monzos[i][l] ?? 0) - (monzos[j][l] ?? 0);
+                const monzo = [...monzos[j]];
+                for (let m = monzos[j].length; m < monzos[i].length; ++m) {
+                  monzo[m] = 0;
+                }
+                const otherWay = [...monzo];
+                monzo[k] += gap;
+                otherWay[l] += otherGap;
+                let monzoUnique = true;
+                let otherUnique = true;
+                for (const existing of monzos.concat(connectingMonzos)) {
+                  if (monzosEqual(monzo, existing)) {
+                    monzoUnique = false;
+                  }
+                  if (monzosEqual(otherWay, existing)) {
+                    otherUnique = false;
+                  }
+                }
+                if (monzoUnique) {
+                  connectingMonzos.push(monzo);
+                }
+                if (otherUnique) {
+                  connectingMonzos.push(otherWay);
+                }
+                break gapSearch;
+              }
+            }
+          }
+        }
+      }
+    }
+  }
+  if (maxDistance >= 1) {
+    const primaryLength = monzos.length;
+    monzos = monzos.concat(connectingMonzos);
+    for (let i = 0; i < monzos.length; ++i) {
+      for (let j = i + 1; j < monzos.length; ++j) {
+        const distance = taxicabDistance(monzos[i], monzos[j]);
+        if (distance === 1) {
+          connections.push({
+            index1: i,
+            index2: j,
+            type:
+              i < primaryLength && j < primaryLength ? 'primary' : 'auxiliary',
+          });
+        }
+      }
+    }
+  }
+  return {
+    connections,
+    connectingMonzos,
+  };
+}
+
+export function project(monzos: number[][], coordss: number[][]) {
+  const projected = monzos.map(m => [...m]);
+  const limit = Math.max(...coordss.map(coords => coords.length));
+  for (const m of projected) {
+    m.length = Math.min(limit, m.length);
+  }
+  for (let i = limit - 1; i >= 0; --i) {
+    if (coordss.map(coords => coords[i]).some(Boolean)) {
+      continue;
+    }
+    for (const m of projected) {
+      m.splice(i, 1);
+    }
+  }
+  return projected;
+}
+
+export function unproject(monzos: number[][], coordss: number[][]) {
+  if (!monzos.length) {
+    return [];
+  }
+  const unprojected = monzos.map(m => [...m]);
+  const limit = Math.max(...coordss.map(coords => coords.length));
+  for (let i = 0; i < limit; ++i) {
+    if (coordss.map(coords => coords[i]).some(Boolean)) {
+      continue;
+    }
+    for (const u of unprojected) {
+      u.splice(i, 0, 0);
+    }
+  }
+  return unprojected;
+}