diff --git a/CMakeLists.txt b/CMakeLists.txt
index 4f02c57..f932dc3 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -10,8 +10,8 @@ set(CMAKE_CXX_EXTENSIONS OFF)
 enable_testing()
 
 # Options
-option(ND_WFC_BUILD_TESTS "Build tests" OFF)  # Temporarily disabled due to hash issue
-option(ND_WFC_BUILD_EXAMPLES "Build examples" OFF)  # Temporarily disabled
+option(ND_WFC_BUILD_TESTS "Build tests" OFF)  # Temporarily disabled due to gtest hash issue
+option(ND_WFC_BUILD_EXAMPLES "Build examples" ON)
 option(ND_WFC_USE_SYSTEM_LIBS "Use system libraries instead of bundled" OFF)
 
 # Add subdirectories
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index 5f02db7..9254f3b 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -1,8 +1,5 @@
 set(EXAMPLE_SOURCES
     basic_wfc_2d.cpp
-    basic_wfc_3d.cpp
-    texture_generation.cpp
-    model_generation.cpp
 )
 
 # Create executables for each example
diff --git a/include/nd-wfc/constraint.hpp b/include/nd-wfc/constraint.hpp
new file mode 100644
index 0000000..6fe74ec
--- /dev/null
+++ b/include/nd-wfc/constraint.hpp
@@ -0,0 +1,125 @@
+#pragma once
+
+#include "types.hpp"
+#include <vector>
+#include <unordered_map>
+#include <unordered_set>
+
+namespace nd_wfc {
+
+// Constraint manager for efficient adjacency rule lookups
+class ConstraintManager {
+private:
+    // Map from (tile, direction) to set of allowed adjacent tiles
+    std::unordered_map<std::uint64_t, std::unordered_set<TileId>> adjacency_rules_;
+    std::vector<TileId> tiles_;
+    std::unordered_map<TileId, Index> tile_to_index_;
+
+public:
+    ConstraintManager(const std::vector<TileId>& tiles, const std::vector<Constraint>& constraints)
+        : tiles_(tiles) {
+        
+        // Build tile index mapping
+        for (Index i = 0; i < tiles_.size(); ++i) {
+            tile_to_index_[tiles_[i]] = i;
+        }
+        
+        // Build adjacency rules from constraints
+        for (const auto& constraint : constraints) {
+            std::uint64_t key = makeKey(constraint.from_tile, constraint.direction);
+            adjacency_rules_[key].insert(constraint.to_tile);
+        }
+    }
+
+    // Check if two tiles can be adjacent in a specific direction
+    bool isAllowed(TileId from_tile, TileId to_tile, Direction direction) const {
+        std::uint64_t key = makeKey(from_tile, direction);
+        auto it = adjacency_rules_.find(key);
+        if (it == adjacency_rules_.end()) {
+            return false; // No rules defined, assume not allowed
+        }
+        return it->second.count(to_tile) > 0;
+    }
+
+    // Get all allowed adjacent tiles for a given tile and direction
+    std::vector<TileId> getAllowedTiles(TileId from_tile, Direction direction) const {
+        std::uint64_t key = makeKey(from_tile, direction);
+        auto it = adjacency_rules_.find(key);
+        if (it == adjacency_rules_.end()) {
+            return {}; // No rules defined
+        }
+        
+        std::vector<TileId> result;
+        result.reserve(it->second.size());
+        for (TileId tile : it->second) {
+            result.push_back(tile);
+        }
+        return result;
+    }
+
+    // Get all allowed adjacent tile indices (for bitset operations)
+    TileMask getAllowedTileMask(TileId from_tile, Direction direction) const {
+        TileMask mask;
+        std::uint64_t key = makeKey(from_tile, direction);
+        auto it = adjacency_rules_.find(key);
+        if (it != adjacency_rules_.end()) {
+            for (TileId tile : it->second) {
+                auto tile_it = tile_to_index_.find(tile);
+                if (tile_it != tile_to_index_.end()) {
+                    mask.set(tile_it->second);
+                }
+            }
+        }
+        return mask;
+    }
+
+    // Get tile index from tile ID
+    Index getTileIndex(TileId tile_id) const {
+        auto it = tile_to_index_.find(tile_id);
+        return (it != tile_to_index_.end()) ? it->second : static_cast<Index>(-1);
+    }
+
+    // Get tile ID from index
+    TileId getTileId(Index index) const {
+        assert(index < tiles_.size());
+        return tiles_[index];
+    }
+
+    // Get all tiles
+    const std::vector<TileId>& getTiles() const { return tiles_; }
+    Index getNumTiles() const { return tiles_.size(); }
+
+    // Validate constraints (basic validation for now)
+    bool validateConstraints() const {
+        // For now, just check that we have some constraints
+        // In a more sophisticated implementation, we could check for symmetry
+        // but for performance reasons, we'll assume the user provides correct constraints
+        return !adjacency_rules_.empty();
+    }
+
+    // Add a constraint
+    void addConstraint(const Constraint& constraint) {
+        std::uint64_t key = makeKey(constraint.from_tile, constraint.direction);
+        adjacency_rules_[key].insert(constraint.to_tile);
+    }
+
+    // Add symmetric constraint (automatically adds the reverse)
+    void addSymmetricConstraint(TileId tile1, TileId tile2, Direction direction) {
+        addConstraint({tile1, tile2, direction});
+        addConstraint({tile2, tile1, getOppositeDirection(direction)});
+    }
+
+private:
+    // Create a unique key for (tile, direction) pairs
+    std::uint64_t makeKey(TileId tile, Direction direction) const {
+        return (static_cast<std::uint64_t>(tile) << 8) | static_cast<std::uint64_t>(direction);
+    }
+
+    // Unpack a key back to (tile, direction)
+    void unpackKey(std::uint64_t key, TileId& tile, Direction& direction) const {
+        tile = static_cast<TileId>(key >> 8);
+        direction = static_cast<Direction>(key & 0xFF);
+    }
+};
+
+} // namespace nd_wfc
diff --git a/include/nd-wfc/entropy.hpp b/include/nd-wfc/entropy.hpp
new file mode 100644
index 0000000..acfcaf5
--- /dev/null
+++ b/include/nd-wfc/entropy.hpp
@@ -0,0 +1,174 @@
+#pragma once
+
+#include "types.hpp"
+#include "wave.hpp"
+#include <vector>
+#include <random>
+#include <algorithm>
+
+namespace nd_wfc {
+
+// Entropy calculator for efficient cell selection during WFC
+template<std::size_t N>
+class EntropyCalculator {
+private:
+    const Wave<N>* wave_;
+    mutable std::vector<Index> entropy_sorted_indices_;
+    mutable bool cache_dirty_;
+
+public:
+    explicit EntropyCalculator(const Wave<N>* wave) 
+        : wave_(wave), cache_dirty_(true) {
+        entropy_sorted_indices_.reserve(wave_->getTotalSize());
+    }
+
+    // Find the cell with minimum entropy (heuristic for next cell to collapse)
+    bool findMinEntropyCell(Coord<N>& coord, std::mt19937& rng) const {
+        updateCache();
+        
+        if (entropy_sorted_indices_.empty()) {
+            return false; // All cells collapsed or contradiction
+        }
+        
+        // Get the minimum entropy value
+        Index first_index = entropy_sorted_indices_[0];
+        Coord<N> first_coord = wave_->indexToCoord(first_index);
+        std::uint8_t min_entropy = wave_->getEntropy(first_coord);
+        
+        if (min_entropy == 0) {
+            return false; // Contradiction
+        }
+        
+        // Find all cells with the same minimum entropy
+        std::vector<Index> candidates;
+        for (Index idx : entropy_sorted_indices_) {
+            Coord<N> test_coord = wave_->indexToCoord(idx);
+            if (wave_->getEntropy(test_coord) == min_entropy) {
+                candidates.push_back(idx);
+            } else {
+                break; // Since sorted, no more minimum entropy cells
+            }
+        }
+        
+        // Randomly select among candidates
+        std::uniform_int_distribution<std::size_t> dist(0, candidates.size() - 1);
+        Index chosen_index = candidates[dist(rng)];
+        coord = wave_->indexToCoord(chosen_index);
+        
+        return true;
+    }
+
+    // Find all cells with minimum entropy
+    std::vector<Coord<N>> findMinEntropyCells() const {
+        updateCache();
+        
+        if (entropy_sorted_indices_.empty()) {
+            return {};
+        }
+        
+        // Get the minimum entropy value
+        Index first_index = entropy_sorted_indices_[0];
+        Coord<N> first_coord = wave_->indexToCoord(first_index);
+        std::uint8_t min_entropy = wave_->getEntropy(first_coord);
+        
+        if (min_entropy == 0) {
+            return {}; // Contradiction
+        }
+        
+        std::vector<Coord<N>> result;
+        for (Index idx : entropy_sorted_indices_) {
+            Coord<N> coord = wave_->indexToCoord(idx);
+            if (wave_->getEntropy(coord) == min_entropy) {
+                result.push_back(coord);
+            } else {
+                break; // Since sorted, no more minimum entropy cells
+            }
+        }
+        
+        return result;
+    }
+
+    // Get entropy statistics
+    struct EntropyStats {
+        std::uint8_t min_entropy;
+        std::uint8_t max_entropy;
+        float average_entropy;
+        Index uncollapsed_cells;
+    };
+
+    EntropyStats getStats() const {
+        updateCache();
+        
+        if (entropy_sorted_indices_.empty()) {
+            return {0, 0, 0.0f, 0};
+        }
+        
+        std::uint8_t min_entropy = 255;
+        std::uint8_t max_entropy = 0;
+        std::uint32_t total_entropy = 0;
+        Index uncollapsed_count = entropy_sorted_indices_.size();
+        
+        for (Index idx : entropy_sorted_indices_) {
+            Coord<N> coord = wave_->indexToCoord(idx);
+            std::uint8_t entropy = wave_->getEntropy(coord);
+            
+            min_entropy = std::min(min_entropy, entropy);
+            max_entropy = std::max(max_entropy, entropy);
+            total_entropy += entropy;
+        }
+        
+        float avg_entropy = uncollapsed_count > 0 ? 
+            static_cast<float>(total_entropy) / uncollapsed_count : 0.0f;
+        
+        return {min_entropy, max_entropy, avg_entropy, uncollapsed_count};
+    }
+
+    // Mark cache as dirty (call when wave function changes)
+    void markDirty() {
+        cache_dirty_ = true;
+    }
+
+    // Check if there are any uncollapsed cells
+    bool hasUncollapedCells() const {
+        updateCache();
+        return !entropy_sorted_indices_.empty();
+    }
+
+private:
+    // Update the internal cache of sorted indices
+    void updateCache() const {
+        if (!cache_dirty_) {
+            return;
+        }
+        
+        entropy_sorted_indices_.clear();
+        
+        // Collect all uncollapsed cell indices
+        for (Index i = 0; i < wave_->getTotalSize(); ++i) {
+            Coord<N> coord = wave_->indexToCoord(i);
+            if (!wave_->isCollapsed(coord)) {
+                entropy_sorted_indices_.push_back(i);
+            }
+        }
+        
+        // Sort by entropy (ascending order)
+        std::sort(entropy_sorted_indices_.begin(), entropy_sorted_indices_.end(),
+            [this](Index a, Index b) {
+                Coord<N> coord_a = wave_->indexToCoord(a);
+                Coord<N> coord_b = wave_->indexToCoord(b);
+                std::uint8_t entropy_a = wave_->getEntropy(coord_a);
+                std::uint8_t entropy_b = wave_->getEntropy(coord_b);
+                
+                if (entropy_a != entropy_b) {
+                    return entropy_a < entropy_b;
+                }
+                
+                // For cells with same entropy, use index as tiebreaker for determinism
+                return a < b;
+            });
+        
+        cache_dirty_ = false;
+    }
+};
+
+} // namespace nd_wfc
diff --git a/include/nd-wfc/grid.hpp b/include/nd-wfc/grid.hpp
new file mode 100644
index 0000000..a328f42
--- /dev/null
+++ b/include/nd-wfc/grid.hpp
@@ -0,0 +1,146 @@
+#pragma once
+
+#include "types.hpp"
+#include <vector>
+#include <cassert>
+
+namespace nd_wfc {
+
+// N-dimensional grid for storing tiles and managing coordinate transformations
+template<std::size_t N>
+class Grid {
+private:
+    Size<N> size_;
+    std::vector<TileId> data_;
+    Index total_size_;
+    std::array<Index, N> strides_;
+
+public:
+    explicit Grid(const Size<N>& size) : size_(size) {
+        // Calculate total size and strides for efficient coordinate mapping
+        total_size_ = 1;
+        for (std::size_t i = 0; i < N; ++i) {
+            total_size_ *= size_[i];
+        }
+        
+        // Calculate strides for row-major ordering
+        strides_[N - 1] = 1;
+        for (std::size_t i = N - 1; i > 0; --i) {
+            strides_[i - 1] = strides_[i] * size_[i];
+        }
+        
+        // Initialize with invalid tile ID
+        data_.resize(total_size_, static_cast<TileId>(-1));
+    }
+
+    // Convert n-dimensional coordinate to linear index
+    Index coordToIndex(const Coord<N>& coord) const {
+        Index index = 0;
+        for (std::size_t i = 0; i < N; ++i) {
+            assert(coord[i] < size_[i]);
+            index += coord[i] * strides_[i];
+        }
+        return index;
+    }
+
+    // Convert linear index to n-dimensional coordinate
+    Coord<N> indexToCoord(Index index) const {
+        assert(index < total_size_);
+        Coord<N> coord;
+        for (std::size_t i = 0; i < N; ++i) {
+            coord[i] = (index / strides_[i]) % size_[i];
+        }
+        return coord;
+    }
+
+    // Check if coordinate is within bounds
+    bool isValidCoord(const Coord<N>& coord) const {
+        for (std::size_t i = 0; i < N; ++i) {
+            if (coord[i] >= size_[i]) {
+                return false;
+            }
+        }
+        return true;
+    }
+
+    // Get neighbor coordinate in a specific direction
+    bool getNeighbor(const Coord<N>& coord, Direction dir, Coord<N>& neighbor) const {
+        auto offset = getDirectionOffset<N>(dir);
+        
+        for (std::size_t i = 0; i < N; ++i) {
+            // Handle signed arithmetic carefully to avoid underflow
+            if (offset[i] < 0) {
+                if (coord[i] < static_cast<Index>(-offset[i])) {
+                    return false; // Would underflow
+                }
+                neighbor[i] = coord[i] - static_cast<Index>(-offset[i]);
+            } else {
+                neighbor[i] = coord[i] + static_cast<Index>(offset[i]);
+                if (neighbor[i] >= size_[i]) {
+                    return false; // Out of bounds
+                }
+            }
+        }
+        return true;
+    }
+
+    // Access tile at coordinate
+    TileId& operator[](const Coord<N>& coord) {
+        return data_[coordToIndex(coord)];
+    }
+
+    const TileId& operator[](const Coord<N>& coord) const {
+        return data_[coordToIndex(coord)];
+    }
+
+    // Access tile by linear index
+    TileId& operator[](Index index) {
+        assert(index < total_size_);
+        return data_[index];
+    }
+
+    const TileId& operator[](Index index) const {
+        assert(index < total_size_);
+        return data_[index];
+    }
+
+    // Getters
+    const Size<N>& getSize() const { return size_; }
+    Index getTotalSize() const { return total_size_; }
+    
+    // Get all neighbor coordinates for a given position
+    std::vector<std::pair<Coord<N>, Direction>> getNeighbors(const Coord<N>& coord) const {
+        std::vector<std::pair<Coord<N>, Direction>> neighbors;
+        auto valid_dirs = getValidDirections<N>();
+        
+        for (Direction dir : valid_dirs) {
+            Coord<N> neighbor;
+            if (getNeighbor(coord, dir, neighbor)) {
+                neighbors.emplace_back(neighbor, dir);
+            }
+        }
+        return neighbors;
+    }
+
+    // Clear the grid
+    void clear() {
+        std::fill(data_.begin(), data_.end(), static_cast<TileId>(-1));
+    }
+
+    // Check if all cells are filled
+    bool isComplete() const {
+        for (const auto& tile : data_) {
+            if (tile == static_cast<TileId>(-1)) {
+                return false;
+            }
+        }
+        return true;
+    }
+};
+
+// Type aliases for common dimensions
+using Grid2D = Grid<2>;
+using Grid3D = Grid<3>;
+using Grid4D = Grid<4>;
+
+} // namespace nd_wfc
diff --git a/include/nd-wfc/propagator.hpp b/include/nd-wfc/propagator.hpp
new file mode 100644
index 0000000..8df1ecc
--- /dev/null
+++ b/include/nd-wfc/propagator.hpp
@@ -0,0 +1,185 @@
+#pragma once
+
+#include "types.hpp"
+#include "wave.hpp"
+#include "constraint.hpp"
+#include "grid.hpp"
+#include <queue>
+#include <stack>
+
+namespace nd_wfc {
+
+// Constraint propagator for efficient constraint satisfaction
+template<std::size_t N>
+class Propagator {
+private:
+    Wave<N>* wave_;
+    const ConstraintManager* constraints_;
+    Grid<N> grid_;
+    
+    // Stack for propagation (more cache-friendly than queue for this use case)
+    mutable std::stack<Coord<N>> propagation_stack_;
+
+public:
+    Propagator(Wave<N>* wave, const ConstraintManager* constraints)
+        : wave_(wave), constraints_(constraints), grid_(wave->getSize()) {}
+
+    // Propagate constraints after a cell collapse
+    bool propagate(const Coord<N>& initial_coord) {
+        propagation_stack_.push(initial_coord);
+        
+        while (!propagation_stack_.empty()) {
+            Coord<N> coord = propagation_stack_.top();
+            propagation_stack_.pop();
+            
+            if (!propagateFromCell(coord)) {
+                return false; // Contradiction found
+            }
+        }
+        
+        return true;
+    }
+
+    // Propagate constraints from a specific cell to its neighbors
+    bool propagateFromCell(const Coord<N>& coord) {
+        if (!wave_->isCollapsed(coord)) {
+            return true; // Nothing to propagate from uncollapsed cell
+        }
+        
+        TileId collapsed_tile = wave_->getCollapsedTile(coord);
+        auto valid_directions = getValidDirections<N>();
+        
+        // Check each neighbor
+        for (Direction dir : valid_directions) {
+            Coord<N> neighbor_coord;
+            if (!grid_.getNeighbor(coord, dir, neighbor_coord)) {
+                continue; // No neighbor in this direction
+            }
+            
+            if (wave_->isCollapsed(neighbor_coord)) {
+                continue; // Neighbor already collapsed
+            }
+            
+            // Get allowed tiles for this direction
+            TileMask allowed_mask = constraints_->getAllowedTileMask(collapsed_tile, dir);
+            
+            // Remove disallowed possibilities from neighbor
+            if (!constrainNeighbor(neighbor_coord, allowed_mask)) {
+                return false; // Contradiction
+            }
+        }
+        
+        return true;
+    }
+
+    // Constrain a neighbor cell based on allowed tiles
+    bool constrainNeighbor(const Coord<N>& coord, const TileMask& allowed_mask) {
+        const TileMask& current_possibilities = wave_->getPossibilities(coord);
+        TileMask new_possibilities = current_possibilities & allowed_mask;
+        
+        // Check if any possibilities were removed
+        if (new_possibilities == current_possibilities) {
+            return true; // No change needed
+        }
+        
+        // Check for contradiction
+        if (!new_possibilities.any()) {
+            return false; // No possibilities left
+        }
+        
+        // Apply the constraint by removing disallowed tiles
+        for (Index i = 0; i < constraints_->getNumTiles(); ++i) {
+            if (current_possibilities[i] && !new_possibilities[i]) {
+                TileId tile_id = constraints_->getTileId(i);
+                if (!wave_->removePossibility(coord, tile_id)) {
+                    return false; // Contradiction
+                }
+            }
+        }
+        
+        // If the cell now has only one possibility, it's effectively collapsed
+        if (new_possibilities.count() == 1) {
+            // Find the remaining tile and collapse to it
+            for (Index i = 0; i < constraints_->getNumTiles(); ++i) {
+                if (new_possibilities[i]) {
+                    TileId tile_id = constraints_->getTileId(i);
+                    wave_->collapse(coord, tile_id);
+                    propagation_stack_.push(coord); // Propagate from this newly collapsed cell
+                    break;
+                }
+            }
+        }
+        
+        return true;
+    }
+
+    // Perform full constraint propagation across the entire grid
+    bool propagateAll() {
+        // Find all collapsed cells and propagate from them
+        for (Index i = 0; i < wave_->getTotalSize(); ++i) {
+            Coord<N> coord = wave_->indexToCoord(i);
+            if (wave_->isCollapsed(coord)) {
+                if (!propagateFromCell(coord)) {
+                    return false;
+                }
+            }
+        }
+        
+        // Continue propagation until no more changes
+        while (!propagation_stack_.empty()) {
+            Coord<N> coord = propagation_stack_.top();
+            propagation_stack_.pop();
+            
+            if (!propagateFromCell(coord)) {
+                return false;
+            }
+        }
+        
+        return true;
+    }
+
+    // Check consistency of current wave function state
+    bool checkConsistency() const {
+        // Check each cell against its neighbors
+        for (Index i = 0; i < wave_->getTotalSize(); ++i) {
+            Coord<N> coord = wave_->indexToCoord(i);
+            
+            if (!wave_->isCollapsed(coord)) {
+                continue; // Skip uncollapsed cells
+            }
+            
+            TileId tile = wave_->getCollapsedTile(coord);
+            auto valid_directions = getValidDirections<N>();
+            
+            // Check each neighbor
+            for (Direction dir : valid_directions) {
+                Coord<N> neighbor_coord;
+                if (!grid_.getNeighbor(coord, dir, neighbor_coord)) {
+                    continue; // No neighbor
+                }
+                
+                if (!wave_->isCollapsed(neighbor_coord)) {
+                    // Check if any allowed tiles are still possible for the neighbor
+                    TileMask allowed_mask = constraints_->getAllowedTileMask(tile, dir);
+                    const TileMask& neighbor_possibilities = wave_->getPossibilities(neighbor_coord);
+                    
+                    if (!(neighbor_possibilities & allowed_mask).any()) {
+                        return false; // No valid possibilities for neighbor
+                    }
+                } else {
+                    // Check if the collapsed neighbor is allowed
+                    TileId neighbor_tile = wave_->getCollapsedTile(neighbor_coord);
+                    if (!constraints_->isAllowed(tile, neighbor_tile, dir)) {
+                        return false; // Invalid adjacency
+                    }
+                }
+            }
+        }
+        
+        return true;
+    }
+
+
+};
+
+} // namespace nd_wfc
diff --git a/include/nd-wfc/types.hpp b/include/nd-wfc/types.hpp
new file mode 100644
index 0000000..e54f003
--- /dev/null
+++ b/include/nd-wfc/types.hpp
@@ -0,0 +1,142 @@
+#pragma once
+
+#include <array>
+#include <cstdint>
+#include <vector>
+#include <bitset>
+
+namespace nd_wfc {
+
+// Basic types
+using Index = std::size_t;
+using TileId = std::uint32_t;
+using Weight = float;
+
+// N-dimensional coordinate and size types
+template<std::size_t N>
+using Coord = std::array<Index, N>;
+
+template<std::size_t N>
+using Size = std::array<Index, N>;
+
+// Direction enumeration for adjacency constraints
+enum class Direction : std::uint8_t {
+    // 2D directions
+    North = 0,
+    South = 1,
+    East = 2,
+    West = 3,
+    
+    // 3D directions (for future use)
+    Up = 4,
+    Down = 5,
+    
+    // 4D directions (for future use)
+    Ana = 6,
+    Kata = 7,
+    
+    Count = 8
+};
+
+// Get the opposite direction
+constexpr Direction getOppositeDirection(Direction dir) {
+    switch (dir) {
+        case Direction::North: return Direction::South;
+        case Direction::South: return Direction::North;
+        case Direction::East: return Direction::West;
+        case Direction::West: return Direction::East;
+        case Direction::Up: return Direction::Down;
+        case Direction::Down: return Direction::Up;
+        case Direction::Ana: return Direction::Kata;
+        case Direction::Kata: return Direction::Ana;
+        default: return Direction::North;
+    }
+}
+
+// Get direction offset for n-dimensional coordinates (using signed integers)
+template<std::size_t N>
+constexpr std::array<int, N> getDirectionOffset(Direction dir) {
+    std::array<int, N> offset{};
+    
+    if constexpr (N >= 2) {
+        switch (dir) {
+            case Direction::North: offset[1] = -1; break;
+            case Direction::South: offset[1] = 1; break;
+            case Direction::East: offset[0] = 1; break;
+            case Direction::West: offset[0] = -1; break;
+            default: break;
+        }
+    }
+    
+    if constexpr (N >= 3) {
+        switch (dir) {
+            case Direction::Up: offset[2] = 1; break;
+            case Direction::Down: offset[2] = -1; break;
+            default: break;
+        }
+    }
+    
+    if constexpr (N >= 4) {
+        switch (dir) {
+            case Direction::Ana: offset[3] = 1; break;
+            case Direction::Kata: offset[3] = -1; break;
+            default: break;
+        }
+    }
+    
+    return offset;
+}
+
+// Get valid directions for N dimensions
+template<std::size_t N>
+constexpr std::array<Direction, N * 2> getValidDirections() {
+    if constexpr (N == 2) {
+        return {Direction::North, Direction::South, Direction::East, Direction::West};
+    } else if constexpr (N == 3) {
+        return {Direction::North, Direction::South, Direction::East, Direction::West, 
+                Direction::Up, Direction::Down};
+    } else if constexpr (N == 4) {
+        return {Direction::North, Direction::South, Direction::East, Direction::West,
+                Direction::Up, Direction::Down, Direction::Ana, Direction::Kata};
+    } else {
+        static_assert(N <= 4, "Only up to 4 dimensions supported");
+        return {};
+    }
+}
+
+// Constraint definition
+struct Constraint {
+    TileId from_tile;
+    TileId to_tile;
+    Direction direction;
+    
+    Constraint(TileId from, TileId to, Direction dir) 
+        : from_tile(from), to_tile(to), direction(dir) {}
+};
+
+// WFC algorithm result
+enum class WfcResult : std::uint8_t {
+    Success = 0,
+    Failure = 1,
+    Contradiction = 2,
+    MaxIterationsReached = 3,
+    InvalidInput = 4
+};
+
+// WFC configuration
+struct WfcConfig {
+    std::uint32_t seed = 0;
+    Index max_iterations = 100000;
+    bool use_entropy_heuristic = true;
+    bool use_weighted_selection = false;
+    
+    WfcConfig() = default;
+};
+
+// Maximum number of tiles we support (for bitset optimization)
+constexpr std::size_t MAX_TILES = 64;
+
+// Bitset for tracking possible tiles at each cell
+using TileMask = std::bitset<MAX_TILES>;
+
+} // namespace nd_wfc
diff --git a/include/nd-wfc/wave.hpp b/include/nd-wfc/wave.hpp
new file mode 100644
index 0000000..8fbb863
--- /dev/null
+++ b/include/nd-wfc/wave.hpp
@@ -0,0 +1,276 @@
+#pragma once
+
+#include "types.hpp"
+#include "grid.hpp"
+#include <vector>
+#include <random>
+
+namespace nd_wfc {
+
+// Wave function for tracking possible states at each grid cell
+template<std::size_t N>
+class Wave {
+private:
+    Size<N> size_;
+    Index total_size_;
+    std::vector<TileMask> possibilities_;
+    std::vector<std::uint8_t> entropy_cache_;
+    std::vector<bool> is_collapsed_;
+    std::vector<TileId> tiles_;
+    Index num_tiles_;
+    bool entropy_dirty_;
+
+public:
+    Wave(const Size<N>& size, const std::vector<TileId>& tiles) 
+        : size_(size), tiles_(tiles), num_tiles_(tiles.size()), entropy_dirty_(true) {
+        
+        assert(num_tiles_ <= MAX_TILES);
+        
+        // Calculate total size
+        total_size_ = 1;
+        for (std::size_t i = 0; i < N; ++i) {
+            total_size_ *= size_[i];
+        }
+        
+        // Initialize all cells with all possibilities
+        TileMask all_possible;
+        for (Index i = 0; i < num_tiles_; ++i) {
+            all_possible.set(i);
+        }
+        
+        possibilities_.resize(total_size_, all_possible);
+        entropy_cache_.resize(total_size_, static_cast<std::uint8_t>(num_tiles_));
+        is_collapsed_.resize(total_size_, false);
+    }
+
+    // Convert coordinate to linear index
+    Index coordToIndex(const Coord<N>& coord) const {
+        Index index = 0;
+        Index stride = 1;
+        for (std::size_t i = N; i > 0; --i) {
+            assert(coord[i-1] < size_[i-1]);
+            index += coord[i-1] * stride;
+            stride *= size_[i-1];
+        }
+        return index;
+    }
+
+    // Convert linear index to coordinate
+    Coord<N> indexToCoord(Index index) const {
+        assert(index < total_size_);
+        Coord<N> coord;
+        for (std::size_t i = 0; i < N; ++i) {
+            coord[i] = index % size_[i];
+            index /= size_[i];
+        }
+        return coord;
+    }
+
+    // Get possible tiles at a coordinate
+    const TileMask& getPossibilities(const Coord<N>& coord) const {
+        return possibilities_[coordToIndex(coord)];
+    }
+
+    // Get possible tiles at an index
+    const TileMask& getPossibilities(Index index) const {
+        assert(index < total_size_);
+        return possibilities_[index];
+    }
+
+    // Check if a tile is possible at a coordinate
+    bool isPossible(const Coord<N>& coord, TileId tile_id) const {
+        Index tile_index = getTileIndex(tile_id);
+        return tile_index < num_tiles_ && possibilities_[coordToIndex(coord)][tile_index];
+    }
+
+    // Remove a tile possibility from a coordinate
+    bool removePossibility(const Coord<N>& coord, TileId tile_id) {
+        Index tile_index = getTileIndex(tile_id);
+        if (tile_index >= num_tiles_) return false;
+        
+        Index cell_index = coordToIndex(coord);
+        if (!possibilities_[cell_index][tile_index]) {
+            return false; // Already removed
+        }
+        
+        possibilities_[cell_index][tile_index] = false;
+        entropy_dirty_ = true;
+        
+        // Check for contradiction
+        return possibilities_[cell_index].any();
+    }
+
+    // Collapse a cell to a specific tile
+    bool collapse(const Coord<N>& coord, TileId tile_id) {
+        Index tile_index = getTileIndex(tile_id);
+        if (tile_index >= num_tiles_) return false;
+        
+        Index cell_index = coordToIndex(coord);
+        if (!possibilities_[cell_index][tile_index]) {
+            return false; // Tile not possible
+        }
+        
+        // Clear all possibilities except the chosen one
+        possibilities_[cell_index].reset();
+        possibilities_[cell_index][tile_index] = true;
+        is_collapsed_[cell_index] = true;
+        entropy_dirty_ = true;
+        
+        return true;
+    }
+
+    // Get entropy (number of possible tiles) at a coordinate
+    std::uint8_t getEntropy(const Coord<N>& coord) const {
+        Index cell_index = coordToIndex(coord);
+        if (entropy_dirty_) {
+            // Recalculate entropy cache
+            const_cast<Wave*>(this)->updateEntropyCache();
+        }
+        return entropy_cache_[cell_index];
+    }
+
+    // Check if a cell is collapsed
+    bool isCollapsed(const Coord<N>& coord) const {
+        return is_collapsed_[coordToIndex(coord)];
+    }
+
+    // Get the collapsed tile at a coordinate (only valid if collapsed)
+    TileId getCollapsedTile(const Coord<N>& coord) const {
+        Index cell_index = coordToIndex(coord);
+        assert(is_collapsed_[cell_index]);
+        
+        const auto& mask = possibilities_[cell_index];
+        for (Index i = 0; i < num_tiles_; ++i) {
+            if (mask[i]) {
+                return tiles_[i];
+            }
+        }
+        assert(false); // Should never reach here if properly collapsed
+        return static_cast<TileId>(-1);
+    }
+
+    // Find the cell with minimum entropy (excluding collapsed cells)
+    bool findMinEntropyCell(Coord<N>& coord, std::mt19937& rng) const {
+        if (entropy_dirty_) {
+            const_cast<Wave*>(this)->updateEntropyCache();
+        }
+        
+        std::uint8_t min_entropy = 255;
+        std::vector<Index> candidates;
+        
+        // Find minimum entropy among uncollapsed cells
+        for (Index i = 0; i < total_size_; ++i) {
+            if (!is_collapsed_[i]) {
+                std::uint8_t entropy = entropy_cache_[i];
+                if (entropy == 0) {
+                    return false; // Contradiction found
+                }
+                if (entropy < min_entropy) {
+                    min_entropy = entropy;
+                    candidates.clear();
+                    candidates.push_back(i);
+                } else if (entropy == min_entropy) {
+                    candidates.push_back(i);
+                }
+            }
+        }
+        
+        if (candidates.empty()) {
+            return false; // All cells collapsed or contradiction
+        }
+        
+        // Randomly select among candidates with same entropy
+        std::uniform_int_distribution<std::size_t> dist(0, candidates.size() - 1);
+        Index chosen_index = candidates[dist(rng)];
+        coord = indexToCoord(chosen_index);
+        
+        return true;
+    }
+
+    // Select a random tile from possibilities at a coordinate
+    TileId selectRandomTile(const Coord<N>& coord, std::mt19937& rng) const {
+        Index cell_index = coordToIndex(coord);
+        const auto& mask = possibilities_[cell_index];
+        
+        std::vector<TileId> available_tiles;
+        for (Index i = 0; i < num_tiles_; ++i) {
+            if (mask[i]) {
+                available_tiles.push_back(tiles_[i]);
+            }
+        }
+        
+        if (available_tiles.empty()) {
+            return static_cast<TileId>(-1); // No possibilities
+        }
+        
+        std::uniform_int_distribution<std::size_t> dist(0, available_tiles.size() - 1);
+        return available_tiles[dist(rng)];
+    }
+
+    // Check if the wave function is completely collapsed
+    bool isComplete() const {
+        for (Index i = 0; i < total_size_; ++i) {
+            if (!is_collapsed_[i]) {
+                return false;
+            }
+        }
+        return true;
+    }
+
+    // Check for contradictions (cells with no possibilities)
+    bool hasContradiction() const {
+        for (Index i = 0; i < total_size_; ++i) {
+            if (!is_collapsed_[i] && !possibilities_[i].any()) {
+                return true;
+            }
+        }
+        return false;
+    }
+
+    // Get the size of the grid
+    const Size<N>& getSize() const { return size_; }
+    Index getTotalSize() const { return total_size_; }
+
+    // Reset the wave function to initial state
+    void reset() {
+        TileMask all_possible;
+        for (Index i = 0; i < num_tiles_; ++i) {
+            all_possible.set(i);
+        }
+        
+        std::fill(possibilities_.begin(), possibilities_.end(), all_possible);
+        std::fill(entropy_cache_.begin(), entropy_cache_.end(), static_cast<std::uint8_t>(num_tiles_));
+        std::fill(is_collapsed_.begin(), is_collapsed_.end(), false);
+        entropy_dirty_ = true;
+    }
+
+private:
+    // Get the index of a tile ID in the tiles vector
+    Index getTileIndex(TileId tile_id) const {
+        for (Index i = 0; i < num_tiles_; ++i) {
+            if (tiles_[i] == tile_id) {
+                return i;
+            }
+        }
+        return static_cast<Index>(-1);
+    }
+
+    // Update the entropy cache
+    void updateEntropyCache() {
+        for (Index i = 0; i < total_size_; ++i) {
+            if (is_collapsed_[i]) {
+                entropy_cache_[i] = 1;
+            } else {
+                entropy_cache_[i] = static_cast<std::uint8_t>(possibilities_[i].count());
+            }
+        }
+        entropy_dirty_ = false;
+    }
+};
+
+// Type aliases for common dimensions
+using Wave2D = Wave<2>;
+using Wave3D = Wave<3>;
+using Wave4D = Wave<4>;
+
+} // namespace nd_wfc
diff --git a/include/nd-wfc/wfc.hpp b/include/nd-wfc/wfc.hpp
new file mode 100644
index 0000000..63d4c1d
--- /dev/null
+++ b/include/nd-wfc/wfc.hpp
@@ -0,0 +1,186 @@
+#pragma once
+
+#include "types.hpp"
+#include "grid.hpp"
+#include "wave.hpp"
+#include "constraint.hpp"
+#include "entropy.hpp"
+#include "propagator.hpp"
+#include <random>
+#include <chrono>
+
+namespace nd_wfc {
+
+// Main Wave Function Collapse algorithm implementation
+template<std::size_t N>
+class Wfc {
+private:
+    Size<N> size_;
+    std::vector<TileId> tiles_;
+    ConstraintManager constraints_;
+    WfcConfig config_;
+    
+    // Algorithm state
+    Wave<N> wave_;
+    Grid<N> result_grid_;
+    EntropyCalculator<N> entropy_calc_;
+    Propagator<N> propagator_;
+    
+    // Runtime state
+    std::mt19937 rng_;
+    Index iteration_count_;
+    WfcResult last_result_;
+
+public:
+    Wfc(const Size<N>& size, const std::vector<TileId>& tiles, 
+        const std::vector<Constraint>& constraints, const WfcConfig& config = WfcConfig{})
+        : size_(size), tiles_(tiles), constraints_(tiles, constraints), config_(config),
+          wave_(size, tiles), result_grid_(size), entropy_calc_(&wave_), propagator_(&wave_, &constraints_),
+          iteration_count_(0), last_result_(WfcResult::InvalidInput) {
+        
+        // Initialize random number generator
+        if (config_.seed == 0) {
+            auto now = std::chrono::high_resolution_clock::now();
+            auto duration = now.time_since_epoch();
+            auto seed = std::chrono::duration_cast<std::chrono::nanoseconds>(duration).count();
+            rng_.seed(static_cast<std::uint32_t>(seed));
+        } else {
+            rng_.seed(config_.seed);
+        }
+        
+        // Validate input
+        if (tiles_.empty()) {
+            last_result_ = WfcResult::InvalidInput;
+        } else if (!constraints_.validateConstraints()) {
+            last_result_ = WfcResult::InvalidInput;
+        } else {
+            last_result_ = WfcResult::Success; // Ready to run
+        }
+    }
+
+    // Run the WFC algorithm
+    WfcResult run() {
+        
+        // Reset state
+        wave_.reset();
+        result_grid_.clear();
+        iteration_count_ = 0;
+        entropy_calc_.markDirty();
+        
+        // Initial constraint propagation
+        if (!propagator_.propagateAll()) {
+            last_result_ = WfcResult::Contradiction;
+            return last_result_;
+        }
+        
+        // Main algorithm loop
+        while (!wave_.isComplete() && iteration_count_ < config_.max_iterations) {
+            ++iteration_count_;
+            
+            // Find cell with minimum entropy
+            Coord<N> coord;
+            if (!entropy_calc_.findMinEntropyCell(coord, rng_)) {
+                if (wave_.hasContradiction()) {
+                    last_result_ = WfcResult::Contradiction;
+                } else {
+                    last_result_ = WfcResult::Success; // All cells collapsed
+                }
+                break;
+            }
+            
+            // Select a random tile from possibilities
+            TileId selected_tile = wave_.selectRandomTile(coord, rng_);
+            if (selected_tile == static_cast<TileId>(-1)) {
+                last_result_ = WfcResult::Contradiction;
+                break;
+            }
+            
+            // Collapse the cell
+            if (!wave_.collapse(coord, selected_tile)) {
+                last_result_ = WfcResult::Contradiction;
+                break;
+            }
+            
+            // Propagate constraints
+            if (!propagator_.propagate(coord)) {
+                last_result_ = WfcResult::Contradiction;
+                break;
+            }
+            
+            // Mark entropy cache as dirty
+            entropy_calc_.markDirty();
+        }
+        
+        // Check final state
+        if (iteration_count_ >= config_.max_iterations) {
+            last_result_ = WfcResult::MaxIterationsReached;
+        } else if (last_result_ != WfcResult::Contradiction && wave_.isComplete()) {
+            last_result_ = WfcResult::Success;
+            copyWaveToGrid();
+        }
+        
+        return last_result_;
+    }
+
+    // Get the result grid (only valid after successful run)
+    const Grid<N>& getGrid() const { return result_grid_; }
+
+    // Check if the algorithm completed successfully
+    bool isComplete() const { 
+        return last_result_ == WfcResult::Success && wave_.isComplete(); 
+    }
+
+    // Get iteration count from last run
+    Index getIterationCount() const { return iteration_count_; }
+
+    // Get last result
+    WfcResult getLastResult() const { return last_result_; }
+
+    // Get entropy statistics
+    typename EntropyCalculator<N>::EntropyStats getEntropyStats() const {
+        return entropy_calc_.getStats();
+    }
+
+    // Reset and prepare for another run
+    void reset() {
+        wave_.reset();
+        result_grid_.clear();
+        iteration_count_ = 0;
+        entropy_calc_.markDirty();
+        last_result_ = WfcResult::InvalidInput;
+    }
+
+    // Get the current wave state (for debugging)
+    const Wave<N>& getWave() const { return wave_; }
+
+    // Get configuration
+    const WfcConfig& getConfig() const { return config_; }
+
+    // Update configuration
+    void setConfig(const WfcConfig& config) { 
+        config_ = config; 
+        if (config_.seed != 0) {
+            rng_.seed(config_.seed);
+        }
+    }
+
+private:
+    // Copy the final wave state to the result grid
+    void copyWaveToGrid() {
+        for (Index i = 0; i < wave_.getTotalSize(); ++i) {
+            Coord<N> coord = wave_.indexToCoord(i);
+            if (wave_.isCollapsed(coord)) {
+                result_grid_[coord] = wave_.getCollapsedTile(coord);
+            } else {
+                result_grid_[coord] = static_cast<TileId>(-1); // Should not happen on success
+            }
+        }
+    }
+};
+
+// Type aliases for common dimensions
+using Wfc2D = Wfc<2>;
+using Wfc3D = Wfc<3>;
+using Wfc4D = Wfc<4>;
+
+} // namespace nd_wfc
diff --git a/src/main.cpp b/src/main.cpp
index 93bec98..eaab10c 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -1,6 +1,248 @@
 
 
+#include <iostream>
+#include <vector>
+#include <chrono>
+#include "nd-wfc/wfc.hpp"
+#include "nd-wfc/types.hpp"
+
+void printGrid2D(const nd_wfc::Grid2D& grid) {
+    const auto& size = grid.getSize();
+    std::cout << "\nGenerated 2D Grid (" << size[0] << "x" << size[1] << "):\n";
+    
+    for (nd_wfc::Index y = 0; y < size[1]; ++y) {
+        for (nd_wfc::Index x = 0; x < size[0]; ++x) {
+            nd_wfc::Coord<2> coord = {x, y};
+            nd_wfc::TileId tile = grid[coord];
+            
+            char symbol;
+            switch (tile) {
+                case 0: symbol = '.'; break; // Empty
+                case 1: symbol = 'G'; break; // Grass
+                case 2: symbol = 'W'; break; // Water
+                case 3: symbol = 'M'; break; // Mountain
+                default: symbol = '?'; break;
+            }
+            std::cout << symbol << ' ';
+        }
+        std::cout << '\n';
+    }
+    std::cout << std::endl;
+}
+
+void test2DSimple() {
+    std::cout << "=== Testing 2D WFC - Simple Terrain ===\n";
+    
+    // Create a 8x8 grid
+    nd_wfc::Size<2> size = {8, 8};
+    
+    // Define 4 tile types: Empty(0), Grass(1), Water(2), Mountain(3)
+    std::vector<nd_wfc::TileId> tiles = {0, 1, 2, 3};
+    
+    // Define adjacency constraints for terrain generation
+    std::vector<nd_wfc::Constraint> constraints = {
+        // Empty can connect to anything
+        {0, 0, nd_wfc::Direction::North}, {0, 0, nd_wfc::Direction::South},
+        {0, 0, nd_wfc::Direction::East}, {0, 0, nd_wfc::Direction::West},
+        {0, 1, nd_wfc::Direction::North}, {0, 1, nd_wfc::Direction::South},
+        {0, 1, nd_wfc::Direction::East}, {0, 1, nd_wfc::Direction::West},
+        {0, 2, nd_wfc::Direction::North}, {0, 2, nd_wfc::Direction::South},
+        {0, 2, nd_wfc::Direction::East}, {0, 2, nd_wfc::Direction::West},
+        {0, 3, nd_wfc::Direction::North}, {0, 3, nd_wfc::Direction::South},
+        {0, 3, nd_wfc::Direction::East}, {0, 3, nd_wfc::Direction::West},
+
+        // Grass can connect to Empty, Grass, and Mountain (not Water)
+        {1, 0, nd_wfc::Direction::North}, {1, 0, nd_wfc::Direction::South},
+        {1, 0, nd_wfc::Direction::East}, {1, 0, nd_wfc::Direction::West},
+        {1, 1, nd_wfc::Direction::North}, {1, 1, nd_wfc::Direction::South},
+        {1, 1, nd_wfc::Direction::East}, {1, 1, nd_wfc::Direction::West},
+        {1, 3, nd_wfc::Direction::North}, {1, 3, nd_wfc::Direction::South},
+        {1, 3, nd_wfc::Direction::East}, {1, 3, nd_wfc::Direction::West},
+
+        // Water can only connect to Empty and Water
+        {2, 0, nd_wfc::Direction::North}, {2, 0, nd_wfc::Direction::South},
+        {2, 0, nd_wfc::Direction::East}, {2, 0, nd_wfc::Direction::West},
+        {2, 2, nd_wfc::Direction::North}, {2, 2, nd_wfc::Direction::South},
+        {2, 2, nd_wfc::Direction::East}, {2, 2, nd_wfc::Direction::West},
+
+        // Mountain can connect to Empty, Grass, and Mountain (not Water)
+        {3, 0, nd_wfc::Direction::North}, {3, 0, nd_wfc::Direction::South},
+        {3, 0, nd_wfc::Direction::East}, {3, 0, nd_wfc::Direction::West},
+        {3, 1, nd_wfc::Direction::North}, {3, 1, nd_wfc::Direction::South},
+        {3, 1, nd_wfc::Direction::East}, {3, 1, nd_wfc::Direction::West},
+        {3, 3, nd_wfc::Direction::North}, {3, 3, nd_wfc::Direction::South},
+        {3, 3, nd_wfc::Direction::East}, {3, 3, nd_wfc::Direction::West},
+    };
+    
+    nd_wfc::WfcConfig config;
+    config.seed = 12345;
+    config.max_iterations = 10000;
+    
+    nd_wfc::Wfc2D wfc(size, tiles, constraints, config);
+    
+    auto start_time = std::chrono::high_resolution_clock::now();
+    nd_wfc::WfcResult result = wfc.run();
+    auto end_time = std::chrono::high_resolution_clock::now();
+    
+    auto duration = std::chrono::duration_cast<std::chrono::microseconds>(end_time - start_time);
+    
+    std::cout << "Result: ";
+    switch (result) {
+        case nd_wfc::WfcResult::Success:
+            std::cout << "Success!\n";
+            break;
+        case nd_wfc::WfcResult::Failure:
+            std::cout << "Failure\n";
+            break;
+        case nd_wfc::WfcResult::Contradiction:
+            std::cout << "Contradiction\n";
+            break;
+        case nd_wfc::WfcResult::MaxIterationsReached:
+            std::cout << "Max iterations reached\n";
+            break;
+        case nd_wfc::WfcResult::InvalidInput:
+            std::cout << "Invalid input\n";
+            break;
+    }
+    
+    std::cout << "Iterations: " << wfc.getIterationCount() << "\n";
+    std::cout << "Time: " << duration.count() << " microseconds\n";
+    
+    if (result == nd_wfc::WfcResult::Success) {
+        printGrid2D(wfc.getGrid());
+        
+        auto stats = wfc.getEntropyStats();
+        std::cout << "Final entropy stats - Min: " << static_cast<int>(stats.min_entropy)
+                  << ", Max: " << static_cast<int>(stats.max_entropy)
+                  << ", Uncollapsed: " << stats.uncollapsed_cells << "\n";
+    }
+}
+
+void test2DConstrainedPattern() {
+    std::cout << "\n=== Testing 2D WFC - Constrained Pattern ===\n";
+    
+    // Create a smaller grid for detailed analysis
+    nd_wfc::Size<2> size = {5, 5};
+    
+    // Two tile types: 0 (empty), 1 (filled)
+    std::vector<nd_wfc::TileId> tiles = {0, 1};
+    
+    // Constraints: create a checkerboard-like pattern
+    std::vector<nd_wfc::Constraint> constraints = {
+        // Tile 0 can only connect to tile 1
+        {0, 1, nd_wfc::Direction::North}, {0, 1, nd_wfc::Direction::South},
+        {0, 1, nd_wfc::Direction::East}, {0, 1, nd_wfc::Direction::West},
+        
+        // Tile 1 can only connect to tile 0
+        {1, 0, nd_wfc::Direction::North}, {1, 0, nd_wfc::Direction::South},
+        {1, 0, nd_wfc::Direction::East}, {1, 0, nd_wfc::Direction::West},
+    };
+    
+    nd_wfc::WfcConfig config;
+    config.seed = 54321;
+    config.max_iterations = 1000;
+    
+    nd_wfc::Wfc2D wfc(size, tiles, constraints, config);
+    
+    auto start_time = std::chrono::high_resolution_clock::now();
+    nd_wfc::WfcResult result = wfc.run();
+    auto end_time = std::chrono::high_resolution_clock::now();
+    
+    auto duration = std::chrono::duration_cast<std::chrono::microseconds>(end_time - start_time);
+    
+    std::cout << "Result: ";
+    switch (result) {
+        case nd_wfc::WfcResult::Success:
+            std::cout << "Success!\n";
+            break;
+        case nd_wfc::WfcResult::Contradiction:
+            std::cout << "Contradiction (expected for checkerboard on 5x5)\n";
+            break;
+        default:
+            std::cout << "Other: " << static_cast<int>(result) << "\n";
+            break;
+    }
+    
+    std::cout << "Iterations: " << wfc.getIterationCount() << "\n";
+    std::cout << "Time: " << duration.count() << " microseconds\n";
+    
+    if (result == nd_wfc::WfcResult::Success) {
+        printGrid2D(wfc.getGrid());
+    }
+}
 
 int main() {
+    std::cout << "N-Dimensional Wave Function Collapse - Performance Test\n";
+    std::cout << "======================================================\n";
+    
+    try {
+        test2DSimple();
+        test2DConstrainedPattern();
+        
+        // Performance test with larger grid
+        std::cout << "\n=== Performance Test - Large Grid ===\n";
+        nd_wfc::Size<2> large_size = {20, 20};
+        std::vector<nd_wfc::TileId> tiles = {0, 1, 2, 3};
+        
+        // Same terrain constraints as before
+        std::vector<nd_wfc::Constraint> constraints = {
+            {0, 0, nd_wfc::Direction::North}, {0, 0, nd_wfc::Direction::South},
+            {0, 0, nd_wfc::Direction::East}, {0, 0, nd_wfc::Direction::West},
+            {0, 1, nd_wfc::Direction::North}, {0, 1, nd_wfc::Direction::South},
+            {0, 1, nd_wfc::Direction::East}, {0, 1, nd_wfc::Direction::West},
+            {0, 2, nd_wfc::Direction::North}, {0, 2, nd_wfc::Direction::South},
+            {0, 2, nd_wfc::Direction::East}, {0, 2, nd_wfc::Direction::West},
+            {0, 3, nd_wfc::Direction::North}, {0, 3, nd_wfc::Direction::South},
+            {0, 3, nd_wfc::Direction::East}, {0, 3, nd_wfc::Direction::West},
+            {1, 0, nd_wfc::Direction::North}, {1, 0, nd_wfc::Direction::South},
+            {1, 0, nd_wfc::Direction::East}, {1, 0, nd_wfc::Direction::West},
+            {1, 1, nd_wfc::Direction::North}, {1, 1, nd_wfc::Direction::South},
+            {1, 1, nd_wfc::Direction::East}, {1, 1, nd_wfc::Direction::West},
+            {1, 3, nd_wfc::Direction::North}, {1, 3, nd_wfc::Direction::South},
+            {1, 3, nd_wfc::Direction::East}, {1, 3, nd_wfc::Direction::West},
+            {2, 0, nd_wfc::Direction::North}, {2, 0, nd_wfc::Direction::South},
+            {2, 0, nd_wfc::Direction::East}, {2, 0, nd_wfc::Direction::West},
+            {2, 2, nd_wfc::Direction::North}, {2, 2, nd_wfc::Direction::South},
+            {2, 2, nd_wfc::Direction::East}, {2, 2, nd_wfc::Direction::West},
+            {3, 0, nd_wfc::Direction::North}, {3, 0, nd_wfc::Direction::South},
+            {3, 0, nd_wfc::Direction::East}, {3, 0, nd_wfc::Direction::West},
+            {3, 1, nd_wfc::Direction::North}, {3, 1, nd_wfc::Direction::South},
+            {3, 1, nd_wfc::Direction::East}, {3, 1, nd_wfc::Direction::West},
+            {3, 3, nd_wfc::Direction::North}, {3, 3, nd_wfc::Direction::South},
+            {3, 3, nd_wfc::Direction::East}, {3, 3, nd_wfc::Direction::West},
+        };
+        
+        nd_wfc::WfcConfig config;
+        config.seed = 98765;
+        config.max_iterations = 50000;
+        
+        nd_wfc::Wfc2D large_wfc(large_size, tiles, constraints, config);
+        
+        auto start_time = std::chrono::high_resolution_clock::now();
+        nd_wfc::WfcResult result = large_wfc.run();
+        auto end_time = std::chrono::high_resolution_clock::now();
+        
+        auto duration = std::chrono::duration_cast<std::chrono::milliseconds>(end_time - start_time);
+        
+        std::cout << "Large grid (20x20) result: ";
+        switch (result) {
+            case nd_wfc::WfcResult::Success: std::cout << "Success!\n"; break;
+            case nd_wfc::WfcResult::Contradiction: std::cout << "Contradiction\n"; break;
+            case nd_wfc::WfcResult::MaxIterationsReached: std::cout << "Max iterations\n"; break;
+            default: std::cout << "Other: " << static_cast<int>(result) << "\n"; break;
+        }
+        
+        std::cout << "Iterations: " << large_wfc.getIterationCount() << "\n";
+        std::cout << "Time: " << duration.count() << " milliseconds\n";
+        std::cout << "Performance: " << (large_wfc.getIterationCount() * 1000.0 / duration.count()) 
+                  << " iterations/second\n";
+        
+        std::cout << "\nAll tests completed!\n";
+        
+    } catch (const std::exception& e) {
+        std::cerr << "Error: " << e.what() << std::endl;
+        return 1;
+    }
+    
     return 0;
 }
\ No newline at end of file
diff --git a/thirdparty/CMakeLists.txt b/thirdparty/CMakeLists.txt
index 3183db3..0517910 100644
--- a/thirdparty/CMakeLists.txt
+++ b/thirdparty/CMakeLists.txt
@@ -1,55 +1,63 @@
-# RapidJSON - header-only library
-if(ND_WFC_USE_SYSTEM_LIBS)
-    find_package(RapidJSON REQUIRED)
-else()
-    # Use the git submodule
-    # Disable examples and tests to avoid compilation issues with newer GCC
-    set(RAPIDJSON_BUILD_EXAMPLES OFF CACHE BOOL "Build RapidJSON examples" FORCE)
-    set(RAPIDJSON_BUILD_TESTS OFF CACHE BOOL "Build RapidJSON tests" FORCE)
-    set(RAPIDJSON_BUILD_DOC OFF CACHE BOOL "Build RapidJSON documentation" FORCE)
-    add_subdirectory(rapidjson)
+# Third-party dependencies - optional components for future features
+# For now, we'll make these optional since the core WFC doesn't need them
 
-    # Create a target for RapidJSON if not already created
-    if(NOT TARGET rapidjson)
-        add_library(rapidjson INTERFACE)
-        target_include_directories(rapidjson
-            INTERFACE
-                ${CMAKE_CURRENT_SOURCE_DIR}/rapidjson/include
-        )
+# RapidJSON - header-only library (optional)
+if(ND_WFC_USE_SYSTEM_LIBS)
+    find_package(RapidJSON QUIET)
+    if(RapidJSON_FOUND)
+        message(STATUS "Found system RapidJSON")
+    endif()
+else()
+    # Check if rapidjson submodule exists
+    if(EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/rapidjson/CMakeLists.txt")
+        set(RAPIDJSON_BUILD_EXAMPLES OFF CACHE BOOL "Build RapidJSON examples" FORCE)
+        set(RAPIDJSON_BUILD_TESTS OFF CACHE BOOL "Build RapidJSON tests" FORCE)
+        set(RAPIDJSON_BUILD_DOC OFF CACHE BOOL "Build RapidJSON documentation" FORCE)
+        add_subdirectory(rapidjson)
+        message(STATUS "Using bundled RapidJSON")
+    else()
+        message(STATUS "RapidJSON not found - skipping (optional dependency)")
     endif()
 endif()
 
-# Assimp - 3D model loading
+# Assimp - 3D model loading (optional)
 if(ND_WFC_USE_SYSTEM_LIBS)
-    find_package(assimp REQUIRED)
+    find_package(assimp QUIET)
+    if(assimp_FOUND)
+        message(STATUS "Found system Assimp")
+    endif()
 else()
-    # Use the git submodule
-    set(ASSIMP_BUILD_TESTS OFF CACHE BOOL "Build Assimp tests" FORCE)
-    set(ASSIMP_BUILD_ASSIMP_TOOLS OFF CACHE BOOL "Build Assimp tools" FORCE)
-    set(ASSIMP_BUILD_SAMPLES OFF CACHE BOOL "Build Assimp samples" FORCE)
-    set(ASSIMP_BUILD_DOCS OFF CACHE BOOL "Build Assimp documentation" FORCE)
-    set(ASSIMP_INSTALL OFF CACHE BOOL "Install Assimp" FORCE)
-    add_subdirectory(assimp)
-
-    # Assimp target should be created by the assimp CMakeLists.txt
+    # Check if assimp submodule exists
+    if(EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/assimp/CMakeLists.txt")
+        set(ASSIMP_BUILD_TESTS OFF CACHE BOOL "Build Assimp tests" FORCE)
+        set(ASSIMP_BUILD_ASSIMP_TOOLS OFF CACHE BOOL "Build Assimp tools" FORCE)
+        set(ASSIMP_BUILD_SAMPLES OFF CACHE BOOL "Build Assimp samples" FORCE)
+        set(ASSIMP_BUILD_DOCS OFF CACHE BOOL "Build Assimp documentation" FORCE)
+        set(ASSIMP_INSTALL OFF CACHE BOOL "Install Assimp" FORCE)
+        add_subdirectory(assimp)
+        message(STATUS "Using bundled Assimp")
+    else()
+        message(STATUS "Assimp not found - skipping (optional dependency)")
+    endif()
 endif()
 
-# SDL3 - for graphics and windowing
+# SDL3 - for graphics and windowing (optional)
 if(ND_WFC_USE_SYSTEM_LIBS)
-    find_package(SDL3 REQUIRED)
+    find_package(SDL3 QUIET)
+    if(SDL3_FOUND)
+        message(STATUS "Found system SDL3")
+    endif()
 else()
-    # Use the git submodule
-    set(SDL_SHARED OFF CACHE BOOL "Build SDL as a shared library" FORCE)
-    set(SDL_STATIC ON CACHE BOOL "Build SDL as a static library" FORCE)
-    set(SDL_TEST OFF CACHE BOOL "Build SDL tests" FORCE)
-    set(SDL_VIDEO OFF CACHE BOOL "Enable video subsystem" FORCE)  # Disable video for headless builds
-    set(SDL_GPU OFF CACHE BOOL "Enable GPU subsystem" FORCE)
-    set(SDL_RENDER OFF CACHE BOOL "Enable render subsystem" FORCE)
-    set(SDL_CAMERA OFF CACHE BOOL "Enable camera subsystem" FORCE)
-    set(SDL_UNIX_CONSOLE_BUILD ON CACHE BOOL "Build SDL for console (headless) applications" FORCE)
-    add_subdirectory(SDL)
-
-    # SDL3 target should be created by the SDL CMakeLists.txt
+    # Check if SDL submodule exists
+    if(EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/SDL/CMakeLists.txt")
+        set(SDL_SHARED OFF CACHE BOOL "Build SDL as a shared library" FORCE)
+        set(SDL_STATIC ON CACHE BOOL "Build SDL as a static library" FORCE)
+        set(SDL_TEST OFF CACHE BOOL "Build SDL tests" FORCE)
+        add_subdirectory(SDL)
+        message(STATUS "Using bundled SDL3")
+    else()
+        message(STATUS "SDL3 not found - skipping (optional dependency)")
+    endif()
 endif()
 
 # SDL_image - for image loading (temporarily disabled due to SDL3 dependency issue)