diff --git a/API.md b/API.md
index ba028db..cd03c2c 100644
--- a/API.md
+++ b/API.md
@@ -8,6 +8,12 @@ The API is written in a way to be very similar to the Minetest API.
 Underground caves have varying shapes, and the variable
 `noordstar.registered_shapes` contains all those definitions.
 
+For shapes, the following functions are available:
+
+- `noordstar_caves.register_shape(shape def)` Define a new cave shape
+- `noordstar_caves.unregister_shape(name)` Remove a defined cave shape
+- `noordstar_caves.clear_registered_shapes()` Remove all known cave shapes
+
 The shapes are defined as follows:
 
 ```lua
diff --git a/lua/engine.lua b/lua/engine.lua
index a38eded..3d63f7b 100644
--- a/lua/engine.lua
+++ b/lua/engine.lua
@@ -1,4 +1,7 @@
 
+-- Constants and magic numbers
+local mapgen_buffer = 16
+
 -- Convert 3d relative coordinates to an index on a flat array
 local function from_3d_to_flat(dx, dy, dz, nx, ny)
     return (nx * ny * dz) + (nx * dy) + dx + 1
@@ -31,6 +34,108 @@ local function iter_3d_area(minp, maxp, callback)
     end
 end
 
+-- Helper function to convert a set of coordinates to a readable string
+local function pos_to_str(pos)
+    return "(" .. pos.x .. ", " .. pos.y .. ", " .. pos.z .. " )"
+end
+
+local Flat3dArray = {}
+Flat3dArray.__index = Flat3dArray
+
+function Flat3dArray:new(minp, maxp, arr)
+    local instance = {}
+    setmetatable(instance, Flat3dArray)
+
+    local nx = maxp.x - minp.x + 1
+    local ny = maxp.y - minp.y + 1
+    local nz = maxp.z - minp.z + 1
+
+    if #arr ~= nx * ny * nz then
+        error(
+            "Input array doesn't match dimension lengths: " .. nx .. " x " ..
+            ny .. " x " .. nz .. " = " .. (nx*ny*nz) .. ", but found " .. #arr
+        )
+    end
+
+    instance.nx = nx
+    instance.ny = ny
+    instance.nz = nz
+    instance.minp = minp
+    instance.maxp = maxp
+    instance.arr = arr
+
+    return instance
+end
+
+function Flat3dArray:from_func(minp, maxp, callback)
+    local arr = {}
+
+    iter_3d_area(minp, maxp, function (i, pos)
+        arr[i] = callback(i, pos)
+    end)
+
+    return self:new(minp, maxp, arr)
+end
+
+function Flat3dArray:get_index(i)
+    local out = self.arr[i]
+    if out == nil then
+        error(
+            "Index " .. i .. " not found in array of length " .. #self.arr
+        )
+    end
+    return out
+end
+
+function Flat3dArray:get_pos(pos)
+    self:validate_pos(pos)
+
+    local dx = pos.x - self.minp.x
+    local dy = pos.y - self.minp.y
+    local dz = pos.z - self.minp.z
+
+    return self:get_index(from_3d_to_flat(dx, dy, dz, self.nx, self.ny))
+end
+
+local function is_valid_pos(pos, minp, maxp)
+    if pos.x < minp.x then
+        return false
+    elseif pos.x > maxp.x then
+        return false
+    elseif pos.y < minp.y then
+        return false
+    elseif pos.y > maxp.y then
+        return false
+    elseif pos.z < minp.z then
+        return false
+    elseif pos.z > maxp.z then
+        return false
+    else
+        return true
+    end
+end
+
+function Flat3dArray:valid_pos(pos)
+    return is_valid_pos(pos, self.minp, self.maxp)
+end
+
+function Flat3dArray:validate_pos(pos)
+    if not self:valid_pos(pos) then
+        error(
+            table.concat(
+                { "Position "
+                , pos_to_str(pos)
+                , " out of bounds from minp = "
+                , pos_to_str(self.minp)
+                , ", maxp = "
+                , pos_to_str(self.maxp)
+                },
+                ""
+            )
+        )
+    end
+end
+
 -- Get an enhanced function from the def function that warns us when the
 -- function does not behave properly
 local function enhanced_func(def)
@@ -83,7 +188,7 @@ local function get_flat_from_shape_def(def, minp, maxp)
         noise_flat_map[i] = f(pos, noise_flat_map[i])
     end)
 
-    return noise_flat_map
+    return Flat3dArray:new(minp, maxp, noise_flat_map)
 end
 
 local function get_flat_from_noise_params(minp, maxp, noise_params)
@@ -101,7 +206,7 @@ local function get_flat_from_noise_params(minp, maxp, noise_params)
         p:get_3d_map_flat(minp, buffer)
     end
 
-    return buffer
+    return Flat3dArray:new(minp, maxp, buffer)
 end
 
 -- Based on the number of cave shapes, calculate how quickly connectivity is
@@ -172,18 +277,16 @@ local function get_threshold_flat(minp, maxp)
         }
     end
 
-    local thresholds = {}
-
-    -- Fill the table
-    iter_3d_area(minp, maxp, function(i, pos)
+    -- Create the flat array
+    return Flat3dArray:from_func(minp, maxp, function(i, pos)
         local total = 0
         local count = 0
 
-        local x = connectivity[i]
-        local y = verticality[i]
+        local x = connectivity:get_pos(pos)
+        local y = verticality:get_pos(pos)
 
         for _, n in pairs(noise) do
-            local v = n.noise[i]
+            local v = n.noise:get_pos(pos)
 
             local dx = math.abs(x - n.def.connectivity_point)
             local dy = math.abs(y - n.def.verticality_point)
@@ -197,13 +300,11 @@ local function get_threshold_flat(minp, maxp)
         end
 
         if count <= 0 then
-            thresholds[i] = -1000
+            return -1000
         else
-            thresholds[i] = total / count
+            return total / count
         end
     end)
-
-    return thresholds
 end
 
 
@@ -212,7 +313,7 @@ local old_overworld_min = mcl_vars.mg_overworld_min
 
 -- If another mod doesn't override the maximum world depth, we will assume that
 -- the world depth is the following value.
-local world_depth = -800
+local world_depth = -60
 
 -- Otherwise, this variable can be changed using the following function,
 -- which will also update all the other necessary variables
@@ -246,48 +347,126 @@ local function cave_vastness(pos)
     end
 end
 
-local tpd_yet = false
+-- Get a flat array of nodes that are either in caves or not in caves
+local function get_flat_cave_bools(minp, maxp)
+    local thresholds = get_threshold_flat(minp, maxp)
+
+    return Flat3dArray:from_func(minp, maxp, function(i, pos)
+        return thresholds:get_pos(pos) >= 1 - cave_vastness(pos)
+    end)
+end
+
+local node_type =
+    { unknown = 1  -- Edge of a chunk
+    , floor = 2    -- Floor of a cave
+    , wall = 3     -- Side wall of a cave
+    , roof = 4     -- Roof of a cave
+    , content = 5  -- Air in the cave not adjacent to a wall, floor or roof
+    , stone = 6    -- Underground node not adjacent to a cave
+    }
+
+local function get_flat_cave_node_types(minp, maxp)
+    local bools = get_flat_cave_bools(minp, maxp)
+
+    return Flat3dArray:from_func(minp, maxp, function (i, pos)
+        if not bools:get_pos(pos) then
+            return node_type.stone
+        -- Floor takes precedence
+        elseif pos.y == minp.y then -- Could be floor
+            return node_type.unknown
+        elseif not bools:get_pos({ x = pos.x, y = pos.y - 1, z = pos.z }) then
+            return node_type.floor
+        -- Then roof takes precedence
+        elseif pos.y == maxp.y then -- Could be roof
+            return node_type.unknown
+        elseif not bools:get_pos({ x = pos.x, y = pos.y + 1, z = pos.z }) then
+            return node_type.roof
+        else
+            -- Check for walls
+            local left  = { x = pos.x - 1, y = pos.y, z = pos.z }
+            local right = { x = pos.x + 1, y = pos.y, z = pos.z }
+            local front = { x = pos.x, y = pos.y, z = pos.z - 1 }
+            local back  = { x = pos.x, y = pos.y, z = pos.z + 1 }
+
+            -- Check if the value is near the edge
+            local on_edge = false
+
+            if left.x < minp.x then
+                on_edge = true
+            elseif not bools:get_pos(left) then
+                return node_type.wall
+            end
+            if right.x > maxp.x then
+                on_edge = true
+            elseif not bools:get_pos(right) then
+                return node_type.wall
+            end
+            if front.z < minp.z then
+                on_edge = true
+            elseif not bools:get_pos(front) then
+                return node_type.wall
+            end
+            if back.z > maxp.z then
+                on_edge = true
+            elseif not bools:get_pos(back) then
+                return node_type.wall
+            end
+
+            if on_edge then
+                return node_type.unknown
+            else
+                return node_type.content
+            end
+        end
+    end)
+end
 
 minetest.register_on_generated(function(minp, maxp, blockseed)
-    if maxp.y < world_depth then
-        return
-    end
+    local vminp =
+        { x = minp.x - mapgen_buffer
+        , y = minp.y - mapgen_buffer
+        , z = minp.z - mapgen_buffer
+        }
+    local vmaxp =
+        { x = maxp.x + mapgen_buffer
+        , y = maxp.y + mapgen_buffer
+        , z = maxp.z + mapgen_buffer
+        }
 
     -- Get voxelmanip
     local vm = minetest.get_mapgen_object("voxelmanip")
-    local data = vm:get_data()
+    local flat_data = Flat3dArray:new(vminp, vmaxp, vm:get_data())
     
     -- Get threshold values
-    local thresholds = get_threshold_flat(minp, maxp)
-    local air = minetest.get_content_id("mcl_core:glass")
+    local node_types = get_flat_cave_node_types(vminp, vmaxp)
+    
+    local node_air   = minetest.get_content_id("air")
+    local node_floor = minetest.get_content_id("mcl_core:glass_green")
+    local node_wall  = minetest.get_content_id("mcl_core:glass_purple")
+    local node_roof  = minetest.get_content_id("mcl_core:glass_red")
+    local node_other = minetest.get_content_id("mcl_core:glass")
 
-    local count = 0
+    local nids = Flat3dArray:from_func(vminp, vmaxp, function(i, pos)
+        local nt = node_types:get_pos(pos)
 
-    iter_3d_area(minp, maxp, function(i, pos)
-        local nx = maxp.x - minp.x + 1 + 32
-        local ny = maxp.y - minp.y + 1 + 32
-
-        local dx = pos.x - minp.x
-        local dy = pos.y - minp.y
-        local dz = pos.z - minp.z
-
-
-        local vi = from_3d_to_flat(dx + 16, dy + 16, dz + 16, nx, ny)
-
-        if not data[vi] then
-            error("Vi is not in data (len " .. #data .. "): " .. vi .. " for (dx, dy, dz) = (" .. dx .. ", " .. dy .. ", " .. dz .. ") with minp = (" .. minp.x .. ", " .. minp.y .. ", " .. minp.z .. ") and maxp = (" .. maxp.x .. ", " .. maxp.y .. ", " .. maxp.z .. ") and pos = (" .. pos.x .. ", " .. pos.y .. ", " .. pos.z .. ")")
-        end
-        
-        if thresholds[i] >= (1 - cave_vastness(pos)) then
-            data[vi] = air
-
-            count = count + 1
+        if not is_valid_pos(pos, minp, maxp) then
+            return flat_data:get_pos(pos)
+        elseif nt == node_type.unknown then
+            return node_other
+        elseif nt == node_type.floor then
+            return node_floor
+        elseif nt == node_type.wall then
+            return node_wall
+        elseif nt == node_type.roof then
+            return node_roof
+        elseif nt == node_type.content then
+            return node_air
+        elseif nt == node_type.stone then
+            return flat_data:get_pos(pos)
         end
     end)
 
     -- Write all changes to the Minetest world
-    vm:set_data(data)
+    vm:set_data(nids.arr)
     vm:write_to_map()
-
-    minetest.chat_send_all("Updated " .. count .. " squares")
 end)
diff --git a/lua/shape.lua b/lua/shape.lua
index ea794cf..f28fb8d 100644
--- a/lua/shape.lua
+++ b/lua/shape.lua
@@ -46,3 +46,7 @@ function noordstar_caves.register_shape(def)
         noordstar_caves.registered_shapes[d.name] = d
     end
 end
+
+function noordstar_caves.clear_registered_shapes()
+    noordstar_caves.registered_shapes = {}
+end