diff --git a/src/neighborhood.py b/src/neighborhood.py
index cd539eb..fa9c5f0 100644
--- a/src/neighborhood.py
+++ b/src/neighborhood.py
@@ -2,36 +2,7 @@
 
 
 """
-import enum
-
-
-class Cells(enum.Enum):
-    """
-    Allows for specification of which cells should be considered in a 2D or 3D matrix.
-
-    For example, to specify the CENTER and NORTH cells, we consider CENTER | NORTH. It
-    does not make sense to use FORWARD and BACKWARD in a 2D matrix, as that specifies
-    looking up and down a bitplane for further cells.
-
-    For higher level dimensions, forgo use of this enumeration entirely, as described in
-    the Neighborhood class.
-    """
-
-    # 2D & 3D
-    CENTER      = 1 << 0
-    NORTH       = 1 << 1
-    NORTHEAST   = 1 << 2
-    EAST        = 1 << 3
-    SOUTHEAST   = 1 << 4
-    SOUTH       = 1 << 5
-    SOUTHWEST   = 1 << 6
-    WEST        = 1 << 7
-    NORTHWEST   = 1 << 8
-    NORTH       = 1 << 9
-
-    # 3D
-    FORWARD     = 1 << 10
-    BACKWARD    = 1 << 11
+import itertools
 
 
 class Neighborhood:
@@ -42,13 +13,136 @@ class Neighborhood:
     the basic Moore neighborhood comprises of the 8 cells surrounding the center, but what if we wanted
     these 8 and include the cell north of north? The following enables this:
 
-    ...
-
+    m_neighborhood = Neighborhood.moore(2)
+    m_neighborhood.extend({(-2, 0): True})
 
     This allows indexing at levels beyond 3D, which the Cells enumeration does not allow, though visualization
     at this point isn't possible.
     """
-    def __init__(self, grid):
 
-        pass
+    class NeighborhoodKey:
+        """
+        Allows proper sorting of neighborhoods.
 
+        Lists should be returned in order, where cell's with smaller indices (in most significant axis first)
+        are listed before cell's with larger ones. For example, in a 3D grid, the neighbors corresponding to:
+
+            offsets = (-1, -1, -1), (-1, 1, 0), (-1, 0, -1), and (1, 0, -1)
+
+        are returned in the following order:
+
+            offsets = (-1, -1, -1), (-1, 0, -1), (1, 0, -1), (-1, 1, 0)
+
+        since the z-axis is most significant, followed by the y-axis, and lastly the x-axis.
+        """
+        def __init__(self, obj, *args):
+            self.obj = obj
+        def __lt__(self, other):
+             return self.compare(self.obj, other.obj) < 0
+        def __gt__(self, other):
+            return self.compare(self.obj, other.obj) > 0
+        def __eq__(self, other):
+            return self.compare(self.obj, other.obj) == 0
+        def __le__(self, other):
+            return self.compare(self.obj, other.obj) <= 0
+        def __ge__(self, other):
+            return self.compare(self.obj, other.obj) >= 0
+        def __ne__(self, other):
+            return self.compare(self.obj, other.obj) != 0
+        def compare(self, other):
+            for i in reversed(range(len(a))):
+                if a[i] < b[i]:
+                    return -1
+                elif a[i] > b[i]:
+                    return 1
+            return 0
+
+
+    def __init__(self):
+        """
+        Sets up an empty neighborhood.
+
+        Initially, no cells are included in a neighborhood. All neighborhoods must be extended.
+        Note the offsets have a tuple as a key representing the position being offsetted by, and as a value,
+        the current state the given cell at the offset is checked to be.
+        """
+        self.offsets = {}
+
+
+    def neighbors(self, bit, grid):
+        """
+        Returns all bits in the given neighborhood.
+
+        The returned cells are grouped with the value the cell is checked to be (a 2-tuple (Bit, value) pair).
+        These are sorted based on the NeighborhoodKey comparison class defined above.
+        """
+        bits = []
+        for k in sorted(self.offsets.keys()):
+            bits.append((k, self.offsets[k]))
+
+        return bits
+
+
+    def extend(self, offsets, strict=False):
+        """
+        Adds new offsets to the instance member offsets.
+
+        We complain if the strict flag is set to True and an offset has already been declared with a different value.
+        """
+        if not strict:
+            self.offsets.update(offsets)
+        else:
+            for k in offsets.keys():
+                value = self.offsets.get(k, None)
+                if value is None:
+                    self.offsets[k] = offsets[k]
+                elif value != offsets[k]:
+                    raise KeyError
+
+
+    @classmethod
+    def moore(cls, dimen, value=True):
+        """
+        Returns a neighborhood corresponding to the Moore neighborhood.
+
+        The Moore neighborhood consists of all adjacent cells. In 2D, these correspond to the 8 touching cells
+        N, NE, E, SE, S, SW, S, and NW. In 3D, this corresponds to all cells in the "backward" and "forward"
+        layer that adjoin the nine cells in the "center" layer. This concept can be extended to N dimensions.
+
+        Note the center cell is excluded, so the total number of offsets are 3^N - 1.
+        """
+        offsets = {}
+        variants = ([-1, 0, 1],) * dimen
+        for current in itertools.product(*variants):
+            if any(current):
+                offsets[current] = value
+
+        m_neighborhood = cls()
+        m_neighborhood.extend(offsets)
+
+        return m_neighborhood
+
+
+    @classmethod
+    def neumann(cls, dimen, value=True):
+        """
+        Returns a neighborhood corresponding to the Von Neumann neighborhood.
+
+        The Von Neumann neighborhood consists of adjacent cells that directly share a face with the current cell.
+        In 2D, these correspond to the 4 touching cells N, S, E, W. In 3D, we include the "backward" and "forward"
+        cell. This concept can be extended to N dimensions.
+
+        Note the center cell is excluded, so the total number of offsets are 2N.
+        """
+        offsets = {}
+        variant = [0] * dimen
+        for i in range(len(variant)):
+            for j in [-1, 1]:
+                variant[i] = j
+                offsets[tuple(variant)] = value
+            variant[i] = 0
+
+        n_neighborhood = cls()
+        n_neighborhood.extend(offsets)
+
+        return n_neighborhood