diff --git a/bindings/python/src/stochastic/CMakeLists.txt b/bindings/python/src/stochastic/CMakeLists.txt
index 842ca0a..75881d9 100644
--- a/bindings/python/src/stochastic/CMakeLists.txt
+++ b/bindings/python/src/stochastic/CMakeLists.txt
@@ -28,6 +28,7 @@ add_geode_python_binding(
         "sampling/direct/double_sampler.hpp"
         "sampling/random_engine.hpp"
         "sampling/distributions.hpp"
+        "spatial/spatial_domain.hpp"
         "stochastic.cpp"
     DEPENDENCIES
         ${PROJECT_NAME}::stochastic
diff --git a/bindings/python/src/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp b/bindings/python/src/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp
index e6e4fb9..6de8a53 100644
--- a/bindings/python/src/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp
+++ b/bindings/python/src/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp
@@ -53,7 +53,7 @@ namespace geode
         pybind11::class_< FractureSimulationRunner,
             std::shared_ptr< FractureSimulationRunner > >(
             module, "FractureSimulationRunner" )
-            .def( pybind11::init< const BoundingBox2D& >(),
+            .def( pybind11::init< const SpatialDomain< 2 >& >(),
                 pybind11::arg( "box" ) )
             .def( "add_fracture_set_descriptor",
                 &FractureSimulationRunner::add_fracture_set_descriptor,
diff --git a/bindings/python/src/stochastic/spatial/spatial_domain.hpp b/bindings/python/src/stochastic/spatial/spatial_domain.hpp
new file mode 100644
index 0000000..db1ea49
--- /dev/null
+++ b/bindings/python/src/stochastic/spatial/spatial_domain.hpp
@@ -0,0 +1,60 @@
+/*
+ * Copyright (c) 2019 - 2025 Geode-solutions
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy
+ * of this software and associated documentation files (the "Software"), to deal
+ * in the Software without restriction, including without limitation the rights
+ * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ * copies of the Software, and to permit persons to whom the Software is
+ * furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included in
+ * all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+ * SOFTWARE.
+ *
+ */
+
+#include <geode/stochastic/spatial/spatial_domain.hpp>
+
+namespace
+{
+    template < geode::index_t dimension >
+    void declare_spatial_domain( pybind11::module& module )
+    {
+        using Domain = geode::SpatialDomain< dimension >;
+        using BBox = geode::BoundingBox< dimension >;
+
+        const auto pyclass_name =
+            "SpatialDomain" + std::to_string( dimension ) + "D";
+
+        pybind11::class_< Domain >( module, pyclass_name.c_str() )
+            .def( pybind11::init< const BBox&, double >(),
+                pybind11::arg( "domain" ), pybind11::arg( "buffer_size" ),
+                R"doc(
+                Create a spatial domain composed of:
+
+                - a core domain (the VOI)
+                - a buffer zone
+                - the extended domain (domain + buffer)
+
+                Arguments:
+                    domain (BoundingBox): main domain / VOI
+                    buffer_size (float): buffer thickness
+            )doc" );
+    }
+} // namespace
+namespace geode
+{
+    void define_spatial_domain( pybind11::module& module )
+    {
+        declare_spatial_domain< 2 >( module );
+        declare_spatial_domain< 3 >( module );
+    }
+} // namespace geode
\ No newline at end of file
diff --git a/bindings/python/src/stochastic/stochastic.cpp b/bindings/python/src/stochastic/stochastic.cpp
index 7e8930a..eb5e689 100644
--- a/bindings/python/src/stochastic/stochastic.cpp
+++ b/bindings/python/src/stochastic/stochastic.cpp
@@ -33,6 +33,7 @@
 
 #include "sampling/distributions.hpp"
 #include "sampling/random_engine.hpp"
+#include "spatial/spatial_domain.hpp"
 
 PYBIND11_MODULE( opengeode_stochastic_py_stochastic, module )
 {
@@ -40,6 +41,8 @@ PYBIND11_MODULE( opengeode_stochastic_py_stochastic, module )
     pybind11::class_< geode::StochasticLibrary >( module, "StochasticLibrary" )
         .def( "initialize", &geode::StochasticLibrary::initialize );
 
+    geode::define_spatial_domain( module );
+
     geode::define_distributions( module );
     geode::define_random_engine( module );
     geode::define_double_sampler( module );
diff --git a/bindings/python/tests/stochastic/test-py-mh-fractures.py b/bindings/python/tests/stochastic/test-py-mh-fractures.py
index f7d466a..544df04 100644
--- a/bindings/python/tests/stochastic/test-py-mh-fractures.py
+++ b/bindings/python/tests/stochastic/test-py-mh-fractures.py
@@ -41,6 +41,7 @@ def test_fracture_simulator():
     box = og.BoundingBox2D()
     box.add_point(og.Point2D([0.0, 0.0]))
     box.add_point(og.Point2D([100.0, 100.0]))
+    domain = stochastic.SpatialDomain2D(box,10.)
 
     # --- Object set
     setA = stochastic.FractureSetDescription()
@@ -60,7 +61,7 @@ def test_fracture_simulator():
     setA.p20 = 0.06
     setA.minimal_spacing = 1.0
 
-    runner = stochastic.FractureSimulationRunner(box)
+    runner = stochastic.FractureSimulationRunner(domain)
     runner.add_fracture_set_descriptor(setA)
     runner.initialize()
 
@@ -90,7 +91,7 @@ def test_two_fracture_sets_simulator():
     box = og.BoundingBox2D()
     box.add_point(og.Point2D([0.0, 0.0]))
     box.add_point(og.Point2D([100.0, 100.0]))
-
+    domain = stochastic.SpatialDomain2D(box,10.)
     # --- Object set A
     setA = stochastic.FractureSetDescription()
     setA.name = "A"
@@ -117,7 +118,7 @@ def test_two_fracture_sets_simulator():
     setB.p20 = 0.03
     setB.minimal_spacing = 2.0
 
-    runner = stochastic.FractureSimulationRunner(box)
+    runner = stochastic.FractureSimulationRunner(domain)
     runner.add_fracture_set_descriptor(setA)
     runner.add_fracture_set_descriptor(setB)
     runner.initialize()
diff --git a/include/geode/stochastic/sampling/direct/object_set_sampler/point_set_sampler.hpp b/include/geode/stochastic/sampling/direct/object_set_sampler/point_set_sampler.hpp
index 391e332..2500aec 100644
--- a/include/geode/stochastic/sampling/direct/object_set_sampler/point_set_sampler.hpp
+++ b/include/geode/stochastic/sampling/direct/object_set_sampler/point_set_sampler.hpp
@@ -28,6 +28,7 @@
 
 #include <geode/stochastic/sampling/direct/object_set_sampler/object_set_sampler.hpp>
 #include <geode/stochastic/sampling/direct/point_uniform_sampler.hpp>
+#include <geode/stochastic/spatial/spatial_domain.hpp>
 
 namespace geode
 {
@@ -35,53 +36,64 @@ namespace geode
     class UniformPointSetSampler : public ObjectSetSampler< Point< dimension > >
     {
     public:
-        UniformPointSetSampler( const BoundingBox< dimension >& box )
-            : ObjectSetSampler< Point< dimension > >{}, box_( box )
+        UniformPointSetSampler( const SpatialDomain< dimension >& domain )
+            : ObjectSetSampler< Point< dimension > >{}, domain_( domain )
         {
-            auto volume = box_.n_volume();
+            auto volume = domain_.extended_n_volume();
             OPENGEODE_EXCEPTION( volume != 0.,
-                "[PointSetSampler] - Undefined Bounding Box (volume ==0)." );
+                "[PointSetSampler] - Undefined Extended Bounding "
+                "Box (volume ==0)." );
             this->log_pdf_ = -std::log( volume );
+            step_move_ = define_step_for_move();
+            OPENGEODE_EXCEPTION( step_move_ > 0.,
+                "[PointSetSampler] - Undefined step length for move (value == ",
+                step_move_, ")." );
         }
 
         Point< dimension > sample( RandomEngine& engine ) const override
         {
-            return PointUniformSampler::sample< dimension >( engine, box_ );
+            return PointUniformSampler::sample< dimension >(
+                engine, domain_.extended_box() );
         }
 
         Point< dimension > change(
             const Point< dimension >& obj, RandomEngine& engine ) const override
         {
-            double ratio = 0.1;
-            geode::Sphere< dimension > ball{ obj,
-                ratio * std::get< 1 >( box_.smallest_length() ) };
+            geode::Sphere< dimension > ball{ obj, step_move_ };
 
             auto new_point =
                 PointUniformSampler::sample< dimension >( engine, ball );
             constexpr index_t max_try{ 100 };
             for( const auto try_id : geode::Range{ max_try } )
             {
-                if( box_.contains( new_point ) )
+                if( domain_.extended_contains( new_point ) )
                 {
                     return new_point;
                 }
                 new_point =
                     PointUniformSampler::sample< dimension >( engine, ball );
             }
-            throw OpenGeodeException( absl::StrCat(
-                "[PointSampler] - Cannot find a point in the box: ",
-                box_.string() ) );
+            throw OpenGeodeException(
+                absl::StrCat( "[PointSampler] - Cannot find a point in the "
+                              "extended domain" ) );
             return obj;
         }
 
     private:
+        double define_step_for_move()
+        {
+            double ratio = 0.1;
+            return ratio * domain_.smallest_length();
+        }
+
         bool is_valid_object( const Point< dimension >& obj ) const override
         {
-            return box_.contains( obj );
+            return domain_.extended_contains( obj );
         }
 
     private:
-        BoundingBox< dimension > box_;
+        const SpatialDomain< dimension >& domain_;
+        double step_move_{ 0. };
     };
 
 } // namespace geode
\ No newline at end of file
diff --git a/include/geode/stochastic/sampling/direct/object_set_sampler/segment_set_sampler.hpp b/include/geode/stochastic/sampling/direct/object_set_sampler/segment_set_sampler.hpp
index 45fedbf..7e8de1a 100644
--- a/include/geode/stochastic/sampling/direct/object_set_sampler/segment_set_sampler.hpp
+++ b/include/geode/stochastic/sampling/direct/object_set_sampler/segment_set_sampler.hpp
@@ -29,30 +29,32 @@
 #include <geode/stochastic/sampling/direct/object_set_sampler/object_set_sampler.hpp>
 #include <geode/stochastic/sampling/direct/point_uniform_sampler.hpp>
 #include <geode/stochastic/sampling/direct/segment_uniform_sampler.hpp>
+#include <geode/stochastic/spatial/spatial_domain.hpp>
 
 namespace geode
 {
     class UniformSegmentSetSampler : public ObjectSetSampler< OwnerSegment2D >
     {
     public:
-        UniformSegmentSetSampler( const BoundingBox< 2 >& box,
+        UniformSegmentSetSampler( const SpatialDomain< 2 >& domain,
             const DoubleSampler::Distribution& length,
             const DoubleSampler::Distribution& azimuth )
             : ObjectSetSampler< OwnerSegment2D >{},
-              box_{ box },
+              domain_{ domain },
               length_{ length },
               azimuth_{ azimuth }
         {
-            auto volume = box_.n_volume();
+            auto volume = domain_.extended_n_volume();
             OPENGEODE_EXCEPTION( volume != 0.,
-                "[SegmentSetSampler] - Undefined Bounding Box (volume == 0)." );
+                "[SegmentSetSampler] - Undefined Extended Bounding "
+                "Box (volume ==0)." );
             this->log_pdf_ = -std::log( volume );
         }
 
         OwnerSegment2D sample( RandomEngine& engine ) const override
         {
             auto seg = SegmentUniformSampler::sample(
-                engine, box_, length_, azimuth_ );
+                engine, domain_.extended_box(), length_, azimuth_ );
             return seg;
         }
 
@@ -72,8 +74,8 @@ namespace geode
             constexpr index_t max_try{ 100 };
             for( const auto try_id : geode::Range{ max_try } )
             {
-                if( box_.contains( new_point )
-                    || box_.contains( extremities[other] ) )
+                if( domain_.extended_contains( new_point )
+                    || domain_.extended_contains( extremities[other] ) )
                 {
                     OwnerSegment2D new_segment{ obj };
                     new_segment.set_point( current, new_point );
@@ -82,8 +84,7 @@ namespace geode
                 new_point = PointUniformSampler::sample< 2 >( engine, ball );
             }
             throw OpenGeodeException( absl::StrCat(
-                "[PointSampler] - Cannot find a point in the box: ",
-                box_.string() ) );
+                "[SegmentSetSampler] - Cannot find a point in the box" ) );
             return obj;
         }
 
@@ -91,13 +92,12 @@ namespace geode
         bool is_valid_object( const OwnerSegment2D& obj ) const override
         {
             const auto& extremities = obj.vertices();
-
-            return box_.contains( extremities[0] )
-                   || box_.contains( extremities[1] );
+            return domain_.extended_contains( extremities[0] )
+                   || domain_.extended_contains( extremities[1] );
         }
 
     private:
-        BoundingBox2D box_;
+        const SpatialDomain< 2 >& domain_;
         DoubleSampler::Distribution length_;
         DoubleSampler::Distribution azimuth_;
     };
diff --git a/include/geode/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp b/include/geode/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp
index 770b2bf..d6523b1 100644
--- a/include/geode/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp
+++ b/include/geode/stochastic/sampling/mcmc/helpers/fracture_simulation_runner.hpp
@@ -32,6 +32,7 @@
 #include <geode/stochastic/sampling/mcmc/models/components/pairwise_term.hpp>
 #include <geode/stochastic/sampling/mcmc/proposal/classical_proposals.hpp>
 #include <geode/stochastic/spatial/pairwise_interactions.hpp>
+#include <geode/stochastic/spatial/spatial_domain.hpp>
 
 namespace geode
 {
@@ -72,7 +73,10 @@ namespace geode
     class FractureSimulationRunner : public SimulationRunner< OwnerSegment2D >
     {
     public:
-        FractureSimulationRunner( const BoundingBox2D& box ) : box_( box ) {}
+        FractureSimulationRunner( const SpatialDomain< 2 >& domain )
+            : SimulationRunner< OwnerSegment2D >( domain )
+        {
+        }
 
         void add_fracture_set_descriptor(
             const FractureSetDescription& descriptor )
@@ -101,7 +105,7 @@ namespace geode
                         set_desc.azimuth );
                 this->set_samplers_.push_back(
                     std::make_unique< UniformSegmentSetSampler >(
-                        box_, length_distribution, azimuth_distribution ) );
+                        domain_, length_distribution, azimuth_distribution ) );
 
                 add_birth_death_change_moves( this->set_samplers_.back(),
                     *proposal_kernel, set_id, set_desc.birth_ratio,
@@ -117,7 +121,8 @@ namespace geode
                     this->energy_terms_collection_.add_energy_term(
                         std::make_unique< DensityTerm< OwnerSegment2D > >(
                             absl::StrCat( set_desc.name, "_density" ),
-                            set_desc.p20, std::vector< uuid >{ set_id } ) ) );
+                            set_desc.p20, std::vector< uuid >{ set_id },
+                            this->domain_ ) ) );
                 // spacing
                 if( set_desc.minimal_spacing < GLOBAL_EPSILON )
                 {
@@ -133,7 +138,7 @@ namespace geode
                         std::make_unique< PairwiseTerm< OwnerSegment2D > >(
                             absl::StrCat( set_desc.name, "_min_spacing" ), 0.,
                             std::vector< uuid >{ set_id },
-                            std::move( interaction ) ) ) );
+                            std::move( interaction ), this->domain_ ) ) );
             }
 
             this->mh_sampler_ =
@@ -159,7 +164,6 @@ namespace geode
         }
 
     private:
-        BoundingBox2D box_;
         std::vector< FractureSetDescription > set_descriptors_;
     };
 
diff --git a/include/geode/stochastic/sampling/mcmc/helpers/simulation_runner.hpp b/include/geode/stochastic/sampling/mcmc/helpers/simulation_runner.hpp
index 38b0a60..908553b 100644
--- a/include/geode/stochastic/sampling/mcmc/helpers/simulation_runner.hpp
+++ b/include/geode/stochastic/sampling/mcmc/helpers/simulation_runner.hpp
@@ -27,6 +27,7 @@
 #include <geode/stochastic/sampling/mcmc/helpers/simulation_printer.hpp>
 #include <geode/stochastic/sampling/mcmc/metropolis_hasting_sampler.hpp>
 #include <geode/stochastic/sampling/mcmc/models/energy_term_collection.hpp>
+#include <geode/stochastic/spatial/spatial_domain.hpp>
 
 #include <absl/strings/str_join.h>
 
@@ -62,7 +63,8 @@ namespace geode
     class SimulationRunner
     {
     public:
-        SimulationRunner() = default;
+        SimulationRunner( const SpatialDomain< ObjectType::dim >& domain )
+            : domain_( domain ) {};
         virtual ~SimulationRunner() = default;
 
         virtual void initialize() = 0;
@@ -153,6 +155,7 @@ namespace geode
         }
 
     protected:
+        SpatialDomain< ObjectType::dim > domain_;
         std::vector< std::unique_ptr< geode::ObjectSetSampler< ObjectType > > >
             set_samplers_;
 
diff --git a/include/geode/stochastic/sampling/mcmc/models/components/density_term.hpp b/include/geode/stochastic/sampling/mcmc/models/components/density_term.hpp
index 967d4a2..7b39ee7 100644
--- a/include/geode/stochastic/sampling/mcmc/models/components/density_term.hpp
+++ b/include/geode/stochastic/sampling/mcmc/models/components/density_term.hpp
@@ -35,9 +35,10 @@ namespace geode
     public:
         explicit DensityTerm( std::string_view name,
             double lambda,
-            std::vector< uuid > targeted_set_ids )
+            std::vector< uuid > targeted_set_ids,
+            const SpatialDomain< ObjectType::dim >& domain )
             : EnergyTerm< ObjectType >(
-                  name, lambda, std::move( targeted_set_ids ) )
+                  name, lambda, std::move( targeted_set_ids ), domain )
         {
         }
 
@@ -50,38 +51,52 @@ namespace geode
         double delta_log_add( const ObjectSets< ObjectType >& /*state*/,
             const ObjectRef< ObjectType >& new_object ) const override
         {
-            if( !this->is_targeted_set( new_object.set_id ) )
+            if( !this->is_targeted_set( new_object.set_id )
+                || !this->is_anchored_in_domain( new_object.object ) )
             {
                 return 0.0;
             }
             return this->contribution( 1.0 );
         }
 
-        double delta_log_remove( const ObjectSets< ObjectType >& /*state*/,
+        double delta_log_remove( const ObjectSets< ObjectType >& state,
             const ObjectId& object_id ) const override
         {
-            if( !this->is_targeted_set( object_id.set_id ) )
+            if( !this->is_targeted_set( object_id.set_id )
+                || !this->is_anchored_in_domain(
+                    state.get_object( object_id ) ) )
             {
                 return 0.0;
             }
             return this->contribution( -1.0 );
         }
 
-        double delta_log_change( const ObjectSets< ObjectType >& /*state*/,
-            const ObjectId& /*old_object_id*/,
-            const ObjectRef< ObjectType >& /*new_object*/ ) const override
+        double delta_log_change( const ObjectSets< ObjectType >& state,
+            const ObjectId& old_object_id,
+            const ObjectRef< ObjectType >& new_object ) const override
         {
-            return 0.0;
+            auto old_in = this->is_anchored_in_domain(
+                state.get_object( old_object_id ) );
+            auto new_in = this->is_anchored_in_domain( new_object.object );
+            if( old_in == new_in )
+            {
+                return 0.0;
+            }
+            return this->contribution( new_in ? 1.0 : -1.0 );
         }
 
         double statistic( const ObjectSets< ObjectType >& state ) const override
         {
-            index_t count{ 0 };
-            for( const auto& set_id : this->targeted_set_ids() )
-            {
-                count += state.nb_objects_in_set( set_id );
-            }
-            return static_cast< double >( count );
+            double sum = 0.0;
+            this->for_each_targeted_object(
+                state, [&]( const ObjectId& obj_id ) {
+                    if( this->is_anchored_in_domain(
+                            state.get_object( obj_id ) ) )
+                    {
+                        sum += 1.0;
+                    }
+                } );
+            return sum;
         }
     };
 } // namespace geode
\ No newline at end of file
diff --git a/include/geode/stochastic/sampling/mcmc/models/components/energy_term.hpp b/include/geode/stochastic/sampling/mcmc/models/components/energy_term.hpp
index 89e5519..dbd3b44 100644
--- a/include/geode/stochastic/sampling/mcmc/models/components/energy_term.hpp
+++ b/include/geode/stochastic/sampling/mcmc/models/components/energy_term.hpp
@@ -30,6 +30,8 @@
 #include <geode/stochastic/common.hpp>
 #include <geode/stochastic/spatial/object_sets.hpp>
 
+#include <geode/stochastic/spatial/spatial_domain.hpp>
+
 #include <optional>
 
 namespace geode
@@ -88,9 +90,11 @@ namespace geode
     public:
         explicit EnergyTerm( std::string_view name,
             double param,
-            std::vector< uuid >&& targeted_set_ids )
+            std::vector< uuid >&& targeted_set_ids,
+            const SpatialDomain< ObjectType::dim >& domain )
             : energy_scale_{ param },
-              targeted_set_ids_{ std::move( targeted_set_ids ) }
+              targeted_set_ids_{ std::move( targeted_set_ids ) },
+              domain_( domain )
         {
             std::sort( targeted_set_ids_.begin(), targeted_set_ids_.end() );
             IdentifierBuilder builder( *this );
@@ -152,6 +156,18 @@ namespace geode
                 targeted_set_ids_.begin(), targeted_set_ids_.end(), set_id );
         }
 
+        bool is_anchored_in_domain( const ObjectType& obj ) const
+        {
+            return SpatialDomainChecker< ObjectType >::is_anchored_in_domain(
+                domain_, obj );
+        }
+
+        bool intersects_domain( const ObjectType& obj ) const
+        {
+            return SpatialDomainChecker< ObjectType >::intersects_domain(
+                domain_, obj );
+        }
+
         template < typename Func >
         void for_each_targeted_object(
             const ObjectSets< ObjectType >& state, Func&& do_apply ) const
@@ -169,5 +185,6 @@ namespace geode
     private:
         detail::EnergyScale energy_scale_;
         std::vector< uuid > targeted_set_ids_;
+        const SpatialDomain< ObjectType::dim > domain_;
     };
 } // namespace geode
\ No newline at end of file
diff --git a/include/geode/stochastic/sampling/mcmc/models/components/pairwise_term.hpp b/include/geode/stochastic/sampling/mcmc/models/components/pairwise_term.hpp
index 6c0b44c..f84162d 100644
--- a/include/geode/stochastic/sampling/mcmc/models/components/pairwise_term.hpp
+++ b/include/geode/stochastic/sampling/mcmc/models/components/pairwise_term.hpp
@@ -38,9 +38,10 @@ namespace geode
         explicit PairwiseTerm( std::string_view name,
             double gamma,
             std::vector< uuid > targeted_set_ids,
-            std::unique_ptr< PairwiseInteraction< ObjectType > > interaction )
+            std::unique_ptr< PairwiseInteraction< ObjectType > > interaction,
+            const SpatialDomain< ObjectType::dim >& domain )
             : EnergyTerm< ObjectType >(
-                  name, gamma, std::move( targeted_set_ids ) ),
+                  name, gamma, std::move( targeted_set_ids ), domain ),
               interaction_( std::move( interaction ) )
         {
         }
@@ -67,8 +68,11 @@ namespace geode
                 geode::ObjectRef< ObjectType > neigh_object{
                     state.get_object( neigh_id ), neigh_id.set_id
                 };
-
-                delta += interaction_->evaluate( new_object, neigh_object );
+                if( this->is_anchored_in_domain( new_object.object )
+                    || this->is_anchored_in_domain( neigh_object.object ) )
+                {
+                    delta += interaction_->evaluate( new_object, neigh_object );
+                }
             }
             return this->contribution( delta );
         }
@@ -92,8 +96,12 @@ namespace geode
                 ObjectRef< ObjectType > neigh_object{
                     state.get_object( neigh_id ), neigh_id.set_id
                 };
-                delta +=
-                    interaction_->evaluate( object_to_remove, neigh_object );
+                if( this->is_anchored_in_domain( object_to_remove.object )
+                    || this->is_anchored_in_domain( neigh_object.object ) )
+                {
+                    delta += interaction_->evaluate(
+                        object_to_remove, neigh_object );
+                }
             }
             return this->contribution( -delta );
         }
@@ -121,8 +129,12 @@ namespace geode
                 ObjectRef< ObjectType > neigh_object{
                     state.get_object( neigh_id ), neigh_id.set_id
                 };
-                delta -=
-                    interaction_->evaluate( object_to_remove, neigh_object );
+                if( this->is_anchored_in_domain( object_to_remove.object )
+                    || this->is_anchored_in_domain( neigh_object.object ) )
+                {
+                    delta -= interaction_->evaluate(
+                        object_to_remove, neigh_object );
+                }
             }
 
             // Add new object's interactions
@@ -138,9 +150,12 @@ namespace geode
                 ObjectRef< ObjectType > neigh_object{
                     state.get_object( neigh_id ), neigh_id.set_id
                 };
-                delta += interaction_->evaluate( new_object, neigh_object );
+                if( this->is_anchored_in_domain( new_object.object )
+                    || this->is_anchored_in_domain( neigh_object.object ) )
+                {
+                    delta += interaction_->evaluate( new_object, neigh_object );
+                }
             }
-
             return this->contribution( delta );
         }
 
@@ -150,18 +165,17 @@ namespace geode
             this->for_each_targeted_object(
                 state, [&]( const ObjectId& obj_id ) {
                     const auto& cur_obj = state.get_object( obj_id );
+                    if( !this->is_anchored_in_domain( cur_obj ) )
+                    {
+                        return;
+                    }
                     ObjectRef< ObjectType > object{ cur_obj, obj_id.set_id };
                     const auto neighbors =
                         state.neighbors( obj_id, this->targeted_set_ids(),
                             interaction_->neighborhood_searching_distance() );
                     for( const auto& neigh_obj_id : neighbors )
                     {
-                        if( neigh_obj_id.set_id < obj_id.set_id )
-                        {
-                            continue;
-                        }
-                        if( neigh_obj_id.set_id == obj_id.set_id
-                            && neigh_obj_id.index <= obj_id.index )
+                        if( !is_new_pair( cur_obj, obj_id, neigh_obj_id ) )
                         {
                             continue;
                         }
@@ -176,6 +190,27 @@ namespace geode
             return sum;
         }
 
+    private:
+        bool is_new_pair( const ObjectType& obj,
+            const ObjectId& obj_id,
+            const ObjectId& neigh_obj_id ) const
+        {
+            if( !this->is_anchored_in_domain( obj ) )
+            {
+                return true;
+            }
+            if( neigh_obj_id.set_id < obj_id.set_id )
+            {
+                return false;
+            }
+            if( neigh_obj_id.set_id == obj_id.set_id
+                && neigh_obj_id.index <= obj_id.index )
+            {
+                return false;
+            }
+            return true;
+        }
+
     private:
         std::unique_ptr< PairwiseInteraction< ObjectType > > interaction_;
     };
diff --git a/include/geode/stochastic/spatial/spatial_domain.hpp b/include/geode/stochastic/spatial/spatial_domain.hpp
new file mode 100644
index 0000000..024b9d3
--- /dev/null
+++ b/include/geode/stochastic/spatial/spatial_domain.hpp
@@ -0,0 +1,133 @@
+/*
+ * Copyright (c) 2019 - 2025 Geode-solutions
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy
+ * of this software and associated documentation files (the "Software"), to deal
+ * in the Software without restriction, including without limitation the rights
+ * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ * copies of the Software, and to permit persons to whom the Software is
+ * furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included in
+ * all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+ * SOFTWARE.
+ *
+ */
+
+#pragma once
+
+#include <geode/geometry/basic_objects/segment.hpp>
+#include <geode/geometry/bounding_box.hpp>
+#include <geode/geometry/point.hpp>
+
+namespace geode
+{
+    template < index_t dimension >
+    class SpatialDomain
+    {
+    public:
+        SpatialDomain(
+            const BoundingBox< dimension >& domain, double buffer_size )
+            : domain_{ domain },
+              buffer_size_{ buffer_size },
+              extended_domain_{ domain }
+        {
+            auto volume = domain_.n_volume();
+            OPENGEODE_EXCEPTION( volume > 0.,
+                "[SpatialDomain] - Undefined Spatial Domain (volume == ",
+                volume, ")." );
+            OPENGEODE_EXCEPTION( buffer_size_ >= 0.0,
+                "[SpatialDomain] buffer size must be non-negative ( buffer "
+                "== ",
+                buffer_size_, ")" );
+            if( buffer_size_ != 0. )
+            {
+                extended_domain_.extends( buffer_size_ );
+            }
+        }
+
+        const BoundingBox< dimension > box() const
+        {
+            return domain_;
+        }
+
+        bool contains( const Point< dimension >& point ) const
+        {
+            return domain_.contains( point );
+        }
+
+        double n_volume() const
+        {
+            return domain_.n_volume();
+        }
+
+        double smallest_length() const
+        {
+            return std::get< 1 >( domain_.smallest_length() );
+        }
+
+        bool extended_contains( const Point< dimension >& point ) const
+        {
+            return extended_domain_.contains( point );
+        }
+
+        double extended_n_volume() const
+        {
+            return extended_domain_.n_volume();
+        }
+
+        const BoundingBox< dimension > extended_box() const
+        {
+            return extended_domain_;
+        }
+
+    private:
+        BoundingBox< dimension > domain_;
+
+        double buffer_size_{ 0. };
+        BoundingBox< dimension > extended_domain_;
+    };
+
+    template < typename ObjectType >
+    struct SpatialDomainChecker
+    {
+        static bool is_anchored_in_domain(
+            const SpatialDomain< ObjectType::dim >& domain,
+            const ObjectType& obj )
+        {
+            return domain.contains( obj );
+        }
+
+        static bool intersects_domain(
+            const SpatialDomain< ObjectType::dim >& domain,
+            const ObjectType& obj )
+        {
+            return domain.contains( obj );
+        }
+    };
+
+    // Specialization for segment
+    template <>
+    struct SpatialDomainChecker< OwnerSegment2D >
+    {
+        static bool is_anchored_in_domain(
+            const SpatialDomain< 2 >& domain, const OwnerSegment2D& seg )
+        {
+            const auto& v = seg.vertices();
+            return domain.contains( v[0] );
+        }
+
+        static bool intersects_domain(
+            const SpatialDomain< 2 >& domain, const OwnerSegment2D& seg )
+        {
+            return domain.box().intersects( seg );
+        }
+    };
+} // namespace geode
\ No newline at end of file
diff --git a/src/geode/stochastic/CMakeLists.txt b/src/geode/stochastic/CMakeLists.txt
index cebfadc..a1baa80 100644
--- a/src/geode/stochastic/CMakeLists.txt
+++ b/src/geode/stochastic/CMakeLists.txt
@@ -40,6 +40,7 @@ add_geode_library(
         "spatial/object_neighborhood.hpp"
         "spatial/object_helpers.hpp"
         "spatial/pairwise_interactions.hpp"
+        "spatial/spatial_domain.hpp"
         "spatial/details/RTree.hpp"
         "sampling/direct/object_set_sampler/object_set_sampler.hpp"
         "sampling/direct/object_set_sampler/point_set_sampler.hpp"
diff --git a/tests/stochastic/CMakeLists.txt b/tests/stochastic/CMakeLists.txt
index 0e8624b..d003af1 100644
--- a/tests/stochastic/CMakeLists.txt
+++ b/tests/stochastic/CMakeLists.txt
@@ -42,6 +42,14 @@ add_geode_test(
         ${PROJECT_NAME}::stochastic
 )
 
+add_geode_test(
+    SOURCE "spatial/test-spatial-domain.cpp"
+    DEPENDENCIES
+        OpenGeode::basic
+        OpenGeode::geometry
+        ${PROJECT_NAME}::stochastic
+)
+
 add_geode_test(
     SOURCE "sampling/direct/test-ball-sampler.cpp"
     DEPENDENCIES
@@ -98,7 +106,7 @@ add_geode_test(
 )
 
 add_geode_test(
-    SOURCE "sampling/mcmc/models/components/test-intensity-term.cpp"
+    SOURCE "sampling/mcmc/models/components/test-density-term.cpp"
     DEPENDENCIES
         OpenGeode::basic
         OpenGeode::geometry
diff --git a/tests/stochastic/sampling/mcmc/models/components/test-density-term.cpp b/tests/stochastic/sampling/mcmc/models/components/test-density-term.cpp
new file mode 100644
index 0000000..b71368d
--- /dev/null
+++ b/tests/stochastic/sampling/mcmc/models/components/test-density-term.cpp
@@ -0,0 +1,151 @@
+/*
+ * Copyright (c) 2019 - 2025 Geode-solutions
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy
+ * of this software and associated documentation files (the "Software"), to deal
+ * in the Software without restriction, including without limitation the rights
+ * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ * copies of the Software, and to permit persons to whom the Software is
+ * furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included in
+ * all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+ * SOFTWARE.
+ *
+ */
+#include <geode/stochastic/sampling/mcmc/models/components/density_term.hpp>
+
+#include <geode/geometry/point.hpp>
+#include <geode/stochastic/spatial/object_sets.hpp>
+#include <geode/stochastic/spatial/spatial_domain.hpp>
+
+geode::uuid init_object_set( geode::ObjectSets< geode::Point2D >& pattern )
+{
+    geode::Point2D p1{ { 0., 0. } };
+    geode::Point2D p2{ { 1., 1. } };
+    geode::Point2D p_buffer{ { 1.3, 0.1 } };
+
+    auto set_id = pattern.add_set( "default_name" );
+    pattern.add_object( std::move( p1 ), set_id );
+    pattern.add_object( std::move( p2 ), set_id );
+    pattern.add_object( std::move( p_buffer ), set_id ); // buffer last
+
+    return set_id;
+}
+
+geode::SpatialDomain< 2 > init_domain()
+{
+    geode::BoundingBox2D box;
+    box.add_point( geode::Point2D{ { -1.0, -1.0 } } );
+    box.add_point( geode::Point2D{ { 1.0, 1.0 } } );
+    return geode::SpatialDomain< 2 >{ box, 0.5 };
+}
+
+void run_density_test( double lambda,
+    const geode::ObjectSets< geode::Point2D >& pattern,
+    const geode::uuid& set_id,
+    const geode::SpatialDomain< 2 >& domain )
+{
+    geode::DensityTerm< geode::Point2D > term(
+        "density", lambda, { set_id }, domain );
+
+    auto neg_log_lambda = -std::log( lambda );
+    double expected_add =
+        ( lambda > 0. ? neg_log_lambda
+                      : std::numeric_limits< double >::infinity() );
+    double expected_remove = ( lambda > 0. ? -neg_log_lambda : 0. );
+
+    // --- Total log
+    double expected_total =
+        ( lambda > 0. ? neg_log_lambda * 2.
+                      : std::numeric_limits< double >::infinity() );
+    double total = term.total_log( pattern );
+    OPENGEODE_EXCEPTION(
+        total == expected_total, "[DensityTerm] total_log wrong" );
+
+    // --- Delta add inside VOI
+    geode::Point2D p_inside{ { 0.5, 0.5 } };
+    geode::ObjectRef< geode::Point2D > ref_inside{ p_inside, set_id };
+    double delta = term.delta_log_add( pattern, ref_inside );
+    OPENGEODE_EXCEPTION(
+        delta == expected_add, "[DensityTerm] delta_log_add inside VOI wrong" );
+
+    // --- Delta add in buffer (outside VOI)
+    geode::Point2D p_buffer{ { 1.3, 0.0 } };
+    geode::ObjectRef< geode::Point2D > ref_buffer{ p_buffer, set_id };
+    delta = term.delta_log_add( pattern, ref_buffer );
+    OPENGEODE_EXCEPTION(
+        delta == 0., "[DensityTerm] delta_log_add outside VOI wrong" );
+
+    // --- Delta remove anchored object
+    geode::ObjectId obj_id{ 0, set_id };
+    delta = term.delta_log_remove( pattern, obj_id );
+    OPENGEODE_EXCEPTION(
+        delta == expected_remove, "[DensityTerm] delta_log_remove wrong" );
+
+    // --- Delta change anchored → buffer
+    geode::ObjectRef< geode::Point2D > new_buffer{ p_buffer, set_id };
+    delta = term.delta_log_change( pattern, obj_id, new_buffer );
+    OPENGEODE_EXCEPTION( delta == expected_remove,
+        "[DensityTerm] delta_log_change anchored→buffer wrong" );
+
+    // --- Delta change anchored → anchored
+    geode::Point2D p_anchored{ { 0.1, 0.1 } };
+    geode::ObjectRef< geode::Point2D > new_anchored{ p_anchored, set_id };
+    delta = term.delta_log_change( pattern, obj_id, new_anchored );
+    OPENGEODE_EXCEPTION(
+        delta == 0., "[DensityTerm] delta_log_change anchored→anchored wrong" );
+
+    // --- Delta change buffer → anchored
+    geode::ObjectId buffer_id{ 2, set_id };
+    delta = term.delta_log_change( pattern, buffer_id, ref_inside );
+    OPENGEODE_EXCEPTION( delta == expected_add,
+        "[DensityTerm] delta_log_change buffer→anchored wrong" );
+}
+
+int main()
+{
+    try
+    {
+        geode::StochasticLibrary::initialize();
+        geode::Logger::set_level( geode::Logger::LEVEL::debug );
+
+        geode::ObjectSets< geode::Point2D > pattern;
+        auto set_id = init_object_set( pattern );
+        auto domain = init_domain();
+
+        // Test different lambda values including near-zero
+        run_density_test( 0.5, pattern, set_id, domain );
+        run_density_test( geode::GLOBAL_EPSILON, pattern, set_id, domain );
+        run_density_test( 100.0021165, pattern, set_id, domain );
+        run_density_test( 0., pattern, set_id, domain ); // zero lambda
+    }
+    catch( ... )
+    {
+        return geode::geode_lippincott();
+    }
+
+    try
+    {
+        geode::StochasticLibrary::initialize();
+        geode::uuid set_id;
+        auto domain = init_domain();
+
+        geode::DensityTerm< geode::Point2D > term(
+            "zero", -geode::GLOBAL_EPSILON, { set_id }, domain );
+        geode::Logger::info( "TEST FAILED" );
+        return 1;
+    }
+    catch( geode::OpenGeodeException& expt )
+    {
+        geode::Logger::info( "TEST SUCCESS" );
+        return 0;
+    }
+}
\ No newline at end of file
diff --git a/tests/stochastic/sampling/mcmc/models/components/test-intensity-term.cpp b/tests/stochastic/sampling/mcmc/models/components/test-intensity-term.cpp
deleted file mode 100644
index 34e3f71..0000000
--- a/tests/stochastic/sampling/mcmc/models/components/test-intensity-term.cpp
+++ /dev/null
@@ -1,131 +0,0 @@
-/*
- * Copyright (c) 2019 - 2025 Geode-solutions
- *
- * Permission is hereby granted, free of charge, to any person obtaining a copy
- * of this software and associated documentation files (the "Software"), to deal
- * in the Software without restriction, including without limitation the rights
- * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
- * copies of the Software, and to permit persons to whom the Software is
- * furnished to do so, subject to the following conditions:
- *
- * The above copyright notice and this permission notice shall be included in
- * all copies or substantial portions of the Software.
- *
- * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
- * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
- * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
- * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
- * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
- * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
- * SOFTWARE.
- *
- */
-#include <geode/stochastic/sampling/mcmc/models/components/density_term.hpp>
-
-#include <geode/geometry/point.hpp>
-#include <geode/stochastic/spatial/object_sets.hpp>
-
-geode::uuid init_object_set( geode::ObjectSets< geode::Point2D > &pattern )
-{
-    geode::Point2D p1{ { 0., 0. } };
-    geode::Point2D p2{ { 1., 1. } };
-
-    auto set_id = pattern.add_set( "default_name" );
-    pattern.add_object( std::move( p1 ), set_id );
-    pattern.add_object( std::move( p2 ), set_id );
-
-    return set_id;
-}
-
-void test_normal_positive_intensity( double lambda,
-    const geode::ObjectSets< geode::Point2D > &pattern,
-    const geode::uuid &set_id )
-{
-    geode::DensityTerm< geode::Point2D > term(
-        "intensity", lambda, { set_id } );
-    auto neg_log_lambda = -std::log( lambda );
-
-    double total = term.total_log( pattern );
-    OPENGEODE_EXCEPTION( total == neg_log_lambda * 2.,
-        "[test intensity]- total_log wrong value." );
-
-    geode::Point2D p3{ { 2., 2. } };
-    geode::ObjectRef< geode::Point2D > p_ref{ p3, set_id };
-    double delta_add = term.delta_log_add( pattern, p_ref );
-    OPENGEODE_EXCEPTION( delta_add == neg_log_lambda * 1.,
-        "[test intensity]- delta_log_add wrong value." );
-
-    geode::ObjectId obj_id{ 0, set_id };
-    double delta_remove = term.delta_log_remove( pattern, obj_id );
-
-    OPENGEODE_EXCEPTION( delta_remove == neg_log_lambda * -1.,
-        "[test intensity]- delta_log_remove wrong value." );
-    double delta_change = term.delta_log_change( pattern, obj_id, p_ref );
-    OPENGEODE_EXCEPTION(
-        delta_change == 0., "[test intensity]- delta_log_change wrong value." );
-}
-
-void test_normal_zero_intensity( double lambda,
-    const geode::ObjectSets< geode::Point2D > &pattern,
-    const geode::uuid &set_id )
-{
-    geode::DensityTerm< geode::Point2D > term(
-        "intensity", lambda, { set_id } );
-    double total = term.total_log( pattern );
-
-    OPENGEODE_EXCEPTION(
-        std::isinf( total ), "[test zero intensity]- total_log wrong value." );
-
-    geode::Point2D p3{ { 2., 2. } };
-    geode::ObjectRef< geode::Point2D > p_ref{ p3, set_id };
-    double delta_add = term.delta_log_add( pattern, p_ref );
-    OPENGEODE_EXCEPTION( std::isinf( delta_add ),
-        "[test zero intensity]- delta_log_add wrong value." );
-
-    geode::ObjectId obj_id{ 0, set_id };
-
-    double delta_remove = term.delta_log_remove( pattern, obj_id );
-    OPENGEODE_EXCEPTION( delta_remove == 0.,
-        "[test zero intensity]- delta_log_remove wrong value." );
-    double delta_change = term.delta_log_change( pattern, obj_id, p_ref );
-    OPENGEODE_EXCEPTION( delta_change == 0.,
-        "[test zero intensity]- delta_log_change wrong value." );
-}
-
-int main()
-{
-    try
-    {
-        geode::StochasticLibrary::initialize();
-        geode::ObjectSets< geode::Point2D > pattern;
-        auto set_id = init_object_set( pattern );
-
-        test_normal_positive_intensity( 0.5, pattern, set_id );
-        test_normal_positive_intensity(
-            geode::GLOBAL_EPSILON, pattern, set_id );
-        test_normal_positive_intensity( 100.0021165, pattern, set_id );
-
-        test_normal_zero_intensity( 0., pattern, set_id );
-        test_normal_zero_intensity(
-            0.9999 * geode::GLOBAL_EPSILON, pattern, set_id );
-    }
-    catch( ... )
-    {
-        return geode::geode_lippincott();
-    }
-
-    try
-    {
-        geode::StochasticLibrary::initialize();
-        geode::uuid set_id;
-        geode::DensityTerm< geode::Point2D > term(
-            "zero", -geode::GLOBAL_EPSILON, { set_id } );
-        geode::Logger::info( "TEST FAILED" );
-        return 1;
-    }
-    catch( geode::OpenGeodeException &expt )
-    {
-        geode::Logger::info( "TEST SUCCESS" );
-        return 0;
-    }
-}
\ No newline at end of file
diff --git a/tests/stochastic/sampling/mcmc/models/components/test-pairwise-term.cpp b/tests/stochastic/sampling/mcmc/models/components/test-pairwise-term.cpp
index 4dce460..51aadb7 100644
--- a/tests/stochastic/sampling/mcmc/models/components/test-pairwise-term.cpp
+++ b/tests/stochastic/sampling/mcmc/models/components/test-pairwise-term.cpp
@@ -27,86 +27,164 @@
 
 #include <geode/geometry/point.hpp>
 #include <geode/stochastic/spatial/object_sets.hpp>
+#include <geode/stochastic/spatial/spatial_domain.hpp>
 
 geode::uuid init_object_set( geode::ObjectSets< geode::Point2D >& pattern )
 {
-    geode::Point2D p1{ { 0., 0. } };
-    geode::Point2D p2{ { 1., 1. } };
+    geode::Point2D p1{ { 0.1, 0.1 } }; // VOI
+    geode::Point2D p2{ { 0.9, 0.9 } }; // VOI
+    geode::Point2D p3{ { 1.3, 1.3 } }; // buffer
 
     auto set_id = pattern.add_set( "default_name" );
     pattern.add_object( std::move( p1 ), set_id );
     pattern.add_object( std::move( p2 ), set_id );
+    pattern.add_object( std::move( p3 ), set_id );
 
     return set_id;
 }
 
-void test_normal_positive_pairwise( double gamma,
+geode::SpatialDomain< 2 > init_domain()
+{
+    geode::BoundingBox2D box;
+    box.add_point( geode::Point2D{ { 0., 0. } } );
+    box.add_point( geode::Point2D{ { 1.0, 1.0 } } );
+    return geode::SpatialDomain< 2 >{ box, 0.5 };
+}
+
+void test_pairwise_term( double gamma,
     const geode::ObjectSets< geode::Point2D >& pattern,
-    const geode::uuid& set_id )
+    const geode::uuid& set_id,
+    const geode::SpatialDomain< 2 >& domain )
 {
     auto interaction =
         std::make_unique< geode::EuclideanCutoffInteraction< geode::Point2D > >(
-            2.1 );
-
+            1 );
     geode::PairwiseTerm< geode::Point2D > term(
-        "strauss", gamma, { set_id }, std::move( interaction ) );
-    auto neg_log_gamma = -std::log( gamma );
+        "strauss", gamma, { set_id }, std::move( interaction ), domain );
 
-    // p1 and p2 interact → 1 pair
-    auto total = term.total_log( pattern );
+    double neg_log_gamma = -std::log( gamma );
 
-    OPENGEODE_EXCEPTION( total == neg_log_gamma * 1.,
-        "[test pairwise] - total_log wrong value." );
-
-    // Adding a third point close to p1 → expect new interactions
-    geode::Point2D p3{ { 0.5, 0.5 } };
-
-    geode::ObjectRef< geode::Point2D > p_ref{ p3, set_id };
-    auto delta_add = term.delta_log_add( pattern, p_ref );
-    // p3 interacts with p1 and p2 → 2 new pairs
-    OPENGEODE_EXCEPTION( delta_add == neg_log_gamma * 2.,
-        "[test pairwise] - delta_log_add wrong value." );
-    geode::ObjectId obj_id{ 0, set_id };
-    auto delta_remove = term.delta_log_remove( pattern, obj_id );
-    // Removing p1 removes its interaction with p2 → 1 removed pair
-    OPENGEODE_EXCEPTION( delta_remove == neg_log_gamma * -1.,
-        "[test pairwise] - delta_log_remove wrong value." );
-
-    auto delta_change = term.delta_log_change( pattern, obj_id, p_ref );
-    // Replacing p1 with p3 changes interactions: p3 interacts with p2 → 1 pair
-    // Old p1 interacted with p2 → 1 pair → no net change
+    // --- total_log
+    double total = term.total_log( pattern );
+    // Only VOI-anchored interactions counted in statistic: p1-p2
     OPENGEODE_EXCEPTION(
-        delta_change == 0., "[test pairwise] - delta_log_change wrong value." );
+        total == term.contribution( 1 ), "[PairwiseTerm] total_log wrong" );
+
+    // --- delta_add VOI → VOI
+    geode::Point2D p4{ { 0.5, 0.5 } };
+    geode::ObjectRef< geode::Point2D > p4_ref{ p4, set_id };
+    double delta = term.delta_log_add( pattern, p4_ref );
+    // interacts with p1 and p2 → 2 pairs
+    OPENGEODE_EXCEPTION( delta == term.contribution( 2 ),
+        "[PairwiseTerm] delta_log_add VOI wrong" );
+
+    // --- delta_add buffer → buffer (outside VOI, inside buffer)
+    geode::Point2D p_buffer{ { 1.2, 1.2 } };
+    geode::ObjectRef< geode::Point2D > buffer_ref{ p_buffer, set_id };
+    delta = term.delta_log_add( pattern, buffer_ref );
+    // energy counts interactions: buffer interacts with p2,p3 (p3 is in
+    // buffer) → 2 pairs
+    OPENGEODE_EXCEPTION( delta == term.contribution( 1 ),
+        "[PairwiseTerm] delta_log_add buffer wrong" );
+
+    // --- delta_change VOI → VOI
+    geode::ObjectId obj_id{ 0, set_id }; // p1
+    delta = term.delta_log_change( pattern, obj_id, p4_ref );
+    // p1 replaced by p4: interacts with p2 → add 1 interaction
+    OPENGEODE_EXCEPTION( delta == term.contribution( 1 ),
+        "[PairwiseTerm] delta_log_change VOI->VOI wrong" );
+
+    // --- delta_change buffer → VOI
+    geode::ObjectId old_buffer{ 2, set_id }; // p3
+    delta = term.delta_log_change( pattern, old_buffer, p4_ref );
+    // old buffer interactions removed (-p3 pairs (-1)), new VOI interactions
+    // added
+    // (+p4 pairs (+2)) net change = 1
+    double expected_delta = term.contribution( 1 ); // adjust based on count
+    OPENGEODE_EXCEPTION( delta == expected_delta,
+        "[PairwiseTerm] delta_log_change buffer->VOI wrong" );
+
+    // --- delta_change VOI → buffer
+    delta = term.delta_log_change( pattern, obj_id, buffer_ref );
+    // p1 → buffer: p1 old interaction = 0
+    // p4  → p_buffer interaction with p2 +1
+    expected_delta = term.contribution( 1 );
+    OPENGEODE_EXCEPTION( delta == expected_delta,
+        "[PairwiseTerm] delta_log_change VOI->buffer wrong" );
+
+    // --- delta_remove VOI
+    delta = term.delta_log_remove( pattern, obj_id );
+    // p1 removed: no interaction with p2 removed
+    OPENGEODE_EXCEPTION( delta == term.contribution( 0 ),
+        "[PairwiseTerm] delta_log_remove VOI wrong" );
+
+    // --- statistic (only anchored objects in VOI)
+    double stat = term.statistic( pattern );
+    // p1,p2 anchored → 1 pair
+    OPENGEODE_EXCEPTION( stat == 1., "[PairwiseTerm] statistic wrong" );
 }
 
-void test_zero_pairwise( double gamma,
+void test_pairwise_term_zero_gamma( double gamma,
     const geode::ObjectSets< geode::Point2D >& pattern,
-    const geode::uuid& set_id )
+    const geode::uuid& set_id,
+    const geode::SpatialDomain< 2 >& domain )
 {
     auto interaction =
         std::make_unique< geode::EuclideanCutoffInteraction< geode::Point2D > >(
-            2.1 );
-
+            1 );
     geode::PairwiseTerm< geode::Point2D > term(
-        "interaction", gamma, { set_id }, std::move( interaction ) );
-
-    auto total = term.total_log( pattern );
-    OPENGEODE_EXCEPTION(
-        std::isinf( total ), "[test zero pairwise] - total_log wrong value." );
-
-    geode::Point2D p3{ { 0.5, 0.5 } };
-    geode::ObjectRef< geode::Point2D > p_ref{ p3, set_id };
-    auto delta_add = term.delta_log_add( pattern, p_ref );
-    OPENGEODE_EXCEPTION( std::isinf( delta_add ),
-        "[test zero pairwise] - delta_log_add wrong value." );
+        "strauss", gamma, { set_id }, std::move( interaction ), domain );
+
+    // --- total_log
+    double total = term.total_log( pattern );
+    OPENGEODE_EXCEPTION( std::isinf( total ),
+        "[PairwiseTerm] total_log with gamma<epsilon should be infinite" );
+
+    // --- delta_add VOI → VOI
+    geode::Point2D p4{ { 0.5, 0.5 } };
+    geode::ObjectRef< geode::Point2D > p4_ref{ p4, set_id };
+    double delta = term.delta_log_add( pattern, p4_ref );
+    OPENGEODE_EXCEPTION( std::isinf( delta ),
+        "[PairwiseTerm] delta_log_add with gamma<epsilon should be infinite" );
+
+    // --- delta_add buffer → buffer
+    geode::Point2D p_buffer{ { 1.2, 1.2 } };
+    geode::ObjectRef< geode::Point2D > buffer_ref{ p_buffer, set_id };
+    delta = term.delta_log_add( pattern, buffer_ref );
+    OPENGEODE_EXCEPTION( std::isinf( delta ),
+        "[PairwiseTerm] delta_log_add buffer with "
+        "gamma<epsilon should be infinite" );
+
+    // --- delta_change VOI → VOI
     geode::ObjectId obj_id{ 0, set_id };
-    auto delta_remove = term.delta_log_remove( pattern, obj_id );
-    OPENGEODE_EXCEPTION( delta_remove == 0.,
-        "[test zero pairwise] - delta_log_remove wrong value." );
-
-    auto delta_change = term.delta_log_change( pattern, obj_id, p_ref );
-    OPENGEODE_EXCEPTION( delta_change == 0.,
-        "[test zero pairwise] - delta_log_change wrong value." );
+    delta = term.delta_log_change( pattern, obj_id, p4_ref );
+    OPENGEODE_EXCEPTION( std::isinf( delta ),
+        "[PairwiseTerm] delta_log_change VOI->VOI with "
+        "gamma<epsilon should be infinite" );
+
+    // --- delta_change buffer → VOI
+    geode::ObjectId old_buffer{ 2, set_id };
+    delta = term.delta_log_change( pattern, old_buffer, p4_ref );
+    OPENGEODE_EXCEPTION( std::isinf( delta ),
+        "[PairwiseTerm] delta_log_change buffer->VOI with "
+        "gamma<epsilon should be infinite" );
+
+    // --- delta_change VOI → buffer
+    delta = term.delta_log_change( pattern, obj_id, buffer_ref );
+    OPENGEODE_EXCEPTION( std::isinf( delta ),
+        "[PairwiseTerm] delta_log_change VOI->buffer with "
+        "gamma<epsilon should be infinite" );
+
+    // --- delta_remove VOI
+    delta = term.delta_log_remove( pattern, obj_id );
+    DEBUG( delta );
+    OPENGEODE_EXCEPTION( delta == 0, "[PairwiseTerm] delta_log_remove VOI with "
+                                     "gamma<epsilon should be infinite" );
+
+    // --- statistic (only anchored objects in VOI)
+    double stat = term.statistic( pattern );
+    // p1,p2 anchored → 1 pair
+    OPENGEODE_EXCEPTION( stat == 1., "[PairwiseTerm] statistic wrong" );
 }
 
 int main()
@@ -118,13 +196,23 @@ int main()
 
         geode::ObjectSets< geode::Point2D > pattern;
         auto set_id = init_object_set( pattern );
+        auto domain = init_domain();
+
+        test_pairwise_term( 0.5, pattern, set_id, domain );
+        test_pairwise_term( geode::GLOBAL_EPSILON, pattern, set_id, domain );
+        test_pairwise_term( 2.0, pattern, set_id, domain );
+
+        test_pairwise_term_zero_gamma(
+            0.9999 * geode::GLOBAL_EPSILON, pattern, set_id, domain );
 
-        test_normal_positive_pairwise( 0.5, pattern, set_id );
-        test_normal_positive_pairwise( geode::GLOBAL_EPSILON, pattern, set_id );
-        test_normal_positive_pairwise( 2.0, pattern, set_id );
+        //            test_normal_positive_pairwise( 0.5, pattern, set_id );
+        //            test_normal_positive_pairwise(
+        //                geode::GLOBAL_EPSILON, pattern, set_id );
+        //            test_normal_positive_pairwise( 2.0, pattern, set_id );
 
-        test_zero_pairwise( 0., pattern, set_id );
-        test_zero_pairwise( 0.9999 * geode::GLOBAL_EPSILON, pattern, set_id );
+        //            test_zero_pairwise( 0., pattern, set_id );
+        //            test_zero_pairwise(
+        //                0.9999 * geode::GLOBAL_EPSILON, pattern, set_id );
         geode::Logger::info( "TEST SUCCESS" );
         return 0;
     }
diff --git a/tests/stochastic/sampling/mcmc/models/test-gibbs-energy.cpp b/tests/stochastic/sampling/mcmc/models/test-gibbs-energy.cpp
index 72f5cb1..6eab15c 100644
--- a/tests/stochastic/sampling/mcmc/models/test-gibbs-energy.cpp
+++ b/tests/stochastic/sampling/mcmc/models/test-gibbs-energy.cpp
@@ -33,6 +33,11 @@
 
 void test_gibbs_energy()
 {
+    geode::BoundingBox2D box;
+    box.add_point( geode::Point2D{ { -1.0, -1.0 } } );
+    box.add_point( geode::Point2D{ { 1.0, 1.0 } } );
+    geode::SpatialDomain domain( box, 0.5 );
+
     geode::ObjectSets< geode::Point2D > pattern;
     const auto set_id = pattern.add_set( "default_name" );
     geode::Point2D p1{ { 0., 0. } };
@@ -45,7 +50,7 @@ void test_gibbs_energy()
     // Add intensity term
     energy_terms.add_energy_term(
         std::make_unique< geode::DensityTerm< geode::Point2D > >(
-            "intensity", 0.5, std::vector< geode::uuid >{ set_id } ) );
+            "intensity", 0.5, std::vector< geode::uuid >{ set_id }, domain ) );
 
     // Add pairwise term with trivial interaction: always counts 1 for each pair
     auto interaction =
@@ -55,7 +60,7 @@ void test_gibbs_energy()
     energy_terms.add_energy_term(
         std::make_unique< geode::PairwiseTerm< geode::Point2D > >(
             "interaction", 0.8, std::vector< geode::uuid >{ set_id },
-            std::move( interaction ) ) );
+            std::move( interaction ), domain ) );
 
     OPENGEODE_EXCEPTION( energy_terms.size() == 2,
         "[test gibbs] Wrong number of components after adding terms." );
diff --git a/tests/stochastic/sampling/mcmc/proposal/test-proposal-kernel.cpp b/tests/stochastic/sampling/mcmc/proposal/test-proposal-kernel.cpp
index 15cd1ad..34178d5 100644
--- a/tests/stochastic/sampling/mcmc/proposal/test-proposal-kernel.cpp
+++ b/tests/stochastic/sampling/mcmc/proposal/test-proposal-kernel.cpp
@@ -56,8 +56,9 @@ namespace
         geode::BoundingBox2D box;
         box.add_point( min_point );
         box.add_point( max_point );
+        geode::SpatialDomain domain( box, 0. );
 
-        geode::UniformPointSetSampler< 2 > sampler( box );
+        geode::UniformPointSetSampler< 2 > sampler( domain );
 
         // Create classical birth-death-change kernel
         auto kernel = geode::create_birth_death_change_kernel< geode::Point2D >(
diff --git a/tests/stochastic/sampling/mcmc/test-metropolis-hasting-sampler.cpp b/tests/stochastic/sampling/mcmc/test-metropolis-hasting-sampler.cpp
index 7509217..2213302 100644
--- a/tests/stochastic/sampling/mcmc/test-metropolis-hasting-sampler.cpp
+++ b/tests/stochastic/sampling/mcmc/test-metropolis-hasting-sampler.cpp
@@ -159,10 +159,11 @@ int main()
         geode::BoundingBox2D box;
         box.add_point( min_point );
         box.add_point( max_point );
+        geode::SpatialDomain domain( box, 0. );
 
         geode::ObjectSets< geode::Point2D > state;
         const auto set_id = state.add_set( "default_name" );
-        geode::UniformPointSetSampler< 2 > sampler( box );
+        geode::UniformPointSetSampler< 2 > sampler( domain );
 
         double birth_prob = 0.3;
         double death_prob = 0.1;
@@ -174,7 +175,8 @@ int main()
         // Add intensity term
         energy_terms.add_energy_term(
             std::make_unique< geode::DensityTerm< geode::Point2D > >(
-                "intensity", 0.5, std::vector< geode::uuid >{ set_id } ) );
+                "intensity", 0.5, std::vector< geode::uuid >{ set_id },
+                domain ) );
 
         geode::MetropolisHastings< geode::Point2D > mh(
             energy_terms, std::move( kernel ) );
diff --git a/tests/stochastic/sampling/mcmc/test-mh-fractures.cpp b/tests/stochastic/sampling/mcmc/test-mh-fractures.cpp
index a8f9d14..1a72134 100644
--- a/tests/stochastic/sampling/mcmc/test-mh-fractures.cpp
+++ b/tests/stochastic/sampling/mcmc/test-mh-fractures.cpp
@@ -36,7 +36,8 @@ namespace
         geode::BoundingBox2D box;
         box.add_point( geode::Point2D{ { 0.0, 0.0 } } );
         box.add_point( geode::Point2D{ { 100.0, 100.0 } } );
-
+        // todo change
+        geode::SpatialDomain domain( box, 0. );
         // --- Object set
         geode::FractureSetDescription setA;
         setA.name = "A";
@@ -57,7 +58,7 @@ namespace
         setA.p20 = 0.05;
         setA.minimal_spacing = 1.;
 
-        geode::FractureSimulationRunner runner( box );
+        geode::FractureSimulationRunner runner( domain );
         runner.add_fracture_set_descriptor( setA );
 
         runner.initialize();
@@ -90,6 +91,8 @@ namespace
         geode::BoundingBox2D box;
         box.add_point( geode::Point2D{ { 0.0, 0.0 } } );
         box.add_point( geode::Point2D{ { 100.0, 100.0 } } );
+        // todo change
+        geode::SpatialDomain domain( box, 0. );
 
         // --- Object set
         geode::FractureSetDescription setA;
@@ -133,7 +136,7 @@ namespace
         setB.p20 = 0.05;
         setB.minimal_spacing = 2.;
 
-        geode::FractureSimulationRunner runner( box );
+        geode::FractureSimulationRunner runner( domain );
         runner.add_fracture_set_descriptor( setA );
         runner.add_fracture_set_descriptor( setB );
 
diff --git a/tests/stochastic/sampling/mcmc/test-mh-poisson.cpp b/tests/stochastic/sampling/mcmc/test-mh-poisson.cpp
index 9f760fe..c1e2b5a 100644
--- a/tests/stochastic/sampling/mcmc/test-mh-poisson.cpp
+++ b/tests/stochastic/sampling/mcmc/test-mh-poisson.cpp
@@ -49,8 +49,8 @@ namespace
         : public geode::SimulationRunner< geode::Point2D >
     {
     public:
-        PoissonSimulationRunner( const geode::BoundingBox2D& box )
-            : box_( box ) {};
+        PoissonSimulationRunner( const geode::SpatialDomain< 2 >& domain )
+            : geode::SimulationRunner< geode::Point2D >( domain ) {};
 
         void add_set_descriptor( const SetDescription& descriptor )
         {
@@ -79,7 +79,7 @@ namespace
 
                 this->set_samplers_.push_back(
                     std::make_unique< geode::UniformPointSetSampler< 2 > >(
-                        box_ ) );
+                        domain_ ) );
 
                 geode::add_birth_death_change_moves( this->set_samplers_.back(),
                     *proposal_kernel, set_id, set_desc.birth_ratio,
@@ -97,7 +97,8 @@ namespace
                             geode::DensityTerm< geode::Point2D > >(
                             absl::StrCat( energy_desc.name, "_density" ),
                             energy_desc.density,
-                            std::vector< geode::uuid >{ set_id } ) ) );
+                            std::vector< geode::uuid >{ set_id },
+                            this->domain_ ) ) );
 
                 this->ordered_target_statistics_.push_back(
                     energy_desc.target_count );
@@ -163,6 +164,7 @@ namespace
         geode::BoundingBox2D box;
         box.add_point( geode::Point2D{ { 0.0, 0.0 } } );
         box.add_point( geode::Point2D{ { 10.0, 10.0 } } );
+        geode::SpatialDomain domain( box, 0. );
 
         std::array< double, 4 > birth_ratio{ 0.1, 0.5, 2., 4. };
         std::array< double, 4 > change_ratio{ 0., 1., 1., 0. };
@@ -182,7 +184,7 @@ namespace
             densityA.density = 0.3;
             densityA.target_count = 30.0;
 
-            PoissonSimulationRunner runner( box );
+            PoissonSimulationRunner runner( domain );
             runner.add_set_descriptor( setA );
             runner.add_density_descriptor( densityA );
             runner.initialize();
@@ -216,6 +218,7 @@ namespace
         geode::BoundingBox2D box;
         box.add_point( geode::Point2D{ { 0.0, 0.0 } } );
         box.add_point( geode::Point2D{ { 10.0, 10.0 } } );
+        geode::SpatialDomain domain( box, 0. );
 
         // --- Set descriptions
         SetDescription set01{ "set01", 2.0, 3.0, 1.0 };
@@ -227,7 +230,7 @@ namespace
         PoissonDensityDescription density02{ "set02", 0.4, 40.0 };
         PoissonDensityDescription density03{ "set03", 0.3, 30.0 };
 
-        PoissonSimulationRunner runner( box );
+        PoissonSimulationRunner runner( domain );
         runner.add_set_descriptor( set01 );
         runner.add_set_descriptor( set02 );
         runner.add_set_descriptor( set03 );
diff --git a/tests/stochastic/sampling/mcmc/test-mh-strauss.cpp b/tests/stochastic/sampling/mcmc/test-mh-strauss.cpp
index 3ce35a4..8a0950c 100644
--- a/tests/stochastic/sampling/mcmc/test-mh-strauss.cpp
+++ b/tests/stochastic/sampling/mcmc/test-mh-strauss.cpp
@@ -61,7 +61,8 @@ namespace
         : public geode::SimulationRunner< geode::Point2D >
     {
     public:
-        StraussSimulationRunner( const geode::BoundingBox2D& box ) : box_( box )
+        StraussSimulationRunner( const geode::SpatialDomain< 2 >& domain )
+            : geode::SimulationRunner< geode::Point2D >( domain )
         {
         }
 
@@ -98,7 +99,7 @@ namespace
 
                 this->set_samplers_.push_back(
                     std::make_unique< geode::UniformPointSetSampler< 2 > >(
-                        box_ ) );
+                        this->domain_ ) );
 
                 geode::add_birth_death_change_moves( this->set_samplers_.back(),
                     *proposal_kernel, set_id, set_desc.birth_ratio,
@@ -115,7 +116,8 @@ namespace
                             geode::DensityTerm< geode::Point2D > >(
                             absl::StrCat( density_desc.name, "_density" ),
                             density_desc.density,
-                            std::vector< geode::uuid >{ set_id } ) ) );
+                            std::vector< geode::uuid >{ set_id },
+                            this->domain_ ) ) );
 
                 this->ordered_target_statistics_.push_back(
                     density_desc.target_count );
@@ -143,7 +145,7 @@ namespace
                                 absl::StrJoin( interaction_desc.names, "_" ),
                                 "_interaction" ),
                             interaction_desc.strength, set_ids,
-                            std::move( interaction ) ) ) );
+                            std::move( interaction ), this->domain_ ) ) );
 
                 this->ordered_target_statistics_.push_back(
                     interaction_desc.target_interaction_count );
@@ -184,7 +186,6 @@ namespace
         }
 
     private:
-        geode::BoundingBox2D box_;
         std::vector< SetDescription > set_descriptors_;
         std::vector< PoissonDensityDescription > density_descriptors_;
         std::vector< PairwiseInteractionDescription > interaction_descriptors_;
@@ -201,10 +202,12 @@ namespace
         geode::BoundingBox2D box;
         box.add_point( geode::Point2D{ { 0.0, 0.0 } } );
         box.add_point( geode::Point2D{ { 10.0, 10.0 } } );
+        // todo change!!
+        geode::SpatialDomain domain( box, 1. );
 
         std::array< double, 5 > gamma_values{ 0, 0.3, 0.5, 0.7, 1.0 };
-        std::array< double, 5 > nb_points{ 22.6, 27.4, 31.3, 36.1, 50. };
-        std::array< double, 5 > nb_interactions{ 0, 4, 8, 14, 36 };
+        std::array< double, 5 > nb_points{ 19.5, 24.4, 31.3, 36.1, 50. };
+        std::array< double, 5 > nb_interactions{ 0, 4, 8, 15, 43 };
 
         for( const auto config : geode::Range{ gamma_values.size() } )
         {
@@ -230,7 +233,7 @@ namespace
             // geode::PairwiseInteraction::SCOPE::INTRA;
             intraA.target_interaction_count = nb_interactions[config];
 
-            StraussSimulationRunner runner( box );
+            StraussSimulationRunner runner( domain );
             runner.add_set_descriptor( setA );
             runner.add_density_descriptor( densityA );
             runner.add_interaction_descriptor( intraA );
@@ -267,6 +270,8 @@ namespace
         geode::BoundingBox2D box;
         box.add_point( geode::Point2D{ { 0.0, 0.0 } } );
         box.add_point( geode::Point2D{ { 10.0, 10.0 } } );
+        // todo change!!
+        geode::SpatialDomain domain( box, 0. );
 
         std::array< double, 3 > gamma_values{ 0, 0.5, 1.0 };
         std::array< double, 3 > nb_points01{ 3.5, 5, 10.0 };
@@ -295,7 +300,7 @@ namespace
             PairwiseInteractionDescription intra02{ { "set02" }, 1., 2.,
                 nb_interactions02[config] };
 
-            StraussSimulationRunner runner( box );
+            StraussSimulationRunner runner( domain );
             runner.add_set_descriptor( set01 );
             runner.add_set_descriptor( set02 );
             runner.add_set_descriptor( set03 );
diff --git a/tests/stochastic/spatial/test-spatial-domain.cpp b/tests/stochastic/spatial/test-spatial-domain.cpp
new file mode 100644
index 0000000..973a185
--- /dev/null
+++ b/tests/stochastic/spatial/test-spatial-domain.cpp
@@ -0,0 +1,142 @@
+/*
+ * Copyright (c) 2019 - 2025 Geode-solutions
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy
+ * of this software and associated documentation files (the "Software"), to deal
+ * in the Software without restriction, including without limitation the rights
+ * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ * copies of the Software, and to permit persons to whom the Software is
+ * furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included in
+ * all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+ * SOFTWARE.
+ *
+ */
+#include <geode/stochastic/common.hpp>
+
+#include <geode/basic/assert.hpp>
+#include <geode/geometry/bounding_box.hpp>
+#include <geode/geometry/point.hpp>
+#include <geode/stochastic/spatial/spatial_domain.hpp>
+
+geode::SpatialDomain< 2 > init_domain()
+{
+    geode::BoundingBox2D box;
+    box.add_point( geode::Point2D{ { 0., 0. } } );
+    box.add_point( geode::Point2D{ { 1.0, 1.0 } } );
+    return geode::SpatialDomain< 2 >{ box, 0.5 };
+}
+
+void test_spatial_domain_2D()
+{
+    using namespace geode;
+
+    auto domain = init_domain();
+
+    // Test volumes
+    OPENGEODE_EXCEPTION(
+        domain.n_volume() == 1.0, "[Test] Domain volume wrong." );
+    OPENGEODE_EXCEPTION( domain.extended_n_volume() == 4,
+        "[Test] Extended domain volume wrong." ); // (1+0.5*2)^2
+
+    // Test points inside domain
+    geode::Point2D inside{ { 0.5, 0.5 } };
+    geode::Point2D boundary{ { 1., 1. } };
+    geode::Point2D outside{ { 1.4, 1.4 } };
+
+    OPENGEODE_EXCEPTION(
+        domain.contains( inside ), "[Test] Point inside should be contained." );
+    OPENGEODE_EXCEPTION( domain.contains( boundary ),
+        "[Test] Boundary point should be contained." );
+    OPENGEODE_EXCEPTION( !domain.contains( outside ),
+        "[Test] Point outside should not be contained." );
+
+    OPENGEODE_EXCEPTION( domain.extended_contains( outside ),
+        "[Test] Point outside original but in buffer should be contained." );
+
+    // Test SpatialDomainChecker for points
+    OPENGEODE_EXCEPTION( SpatialDomainChecker< Point2D >::is_anchored_in_domain(
+                             domain, inside ),
+        "[Test] anchored_in_domain failed for inside point" );
+    OPENGEODE_EXCEPTION(
+        !SpatialDomainChecker< Point2D >::is_anchored_in_domain(
+            domain, outside ),
+        "[Test] anchored_in_domain failed for outside point" );
+    OPENGEODE_EXCEPTION(
+        SpatialDomainChecker< Point2D >::intersects_domain( domain, boundary ),
+        "[Test] intersects_domain failed for boundary point" );
+
+    // Test segments
+    geode::OwnerSegment2D seg_inside{ geode::Point2D{ { 0.1, 0.1 } },
+        geode::Point2D{ { 0.9, 0.9 } } };
+    geode::OwnerSegment2D seg_partial_anchored{ geode::Point2D{ { 0.9, 0.9 } },
+        geode::Point2D{ { 1.6, 1.6 } } };
+    geode::OwnerSegment2D seg_partial{ geode::Point2D{ { -0.5, 0.9 } },
+        geode::Point2D{ { 0.6, 0.6 } } };
+    geode::OwnerSegment2D seg_outside{ geode::Point2D{ { 1.7, 1.7 } },
+        geode::Point2D{ { 2.0, 2.0 } } };
+    geode::OwnerSegment2D seg_cross{ geode::Point2D{ { -0.5, 0.5 } },
+        geode::Point2D{ { 1.5, 0.5 } } };
+
+    OPENGEODE_EXCEPTION(
+        SpatialDomainChecker< OwnerSegment2D >::is_anchored_in_domain(
+            domain, seg_inside ),
+        "[Test] Segment inside VOI should be anchored." );
+    OPENGEODE_EXCEPTION(
+        SpatialDomainChecker< OwnerSegment2D >::is_anchored_in_domain(
+            domain, seg_partial_anchored ),
+        "[Test] Segment with lower left extremity inside VOI should be "
+        "anchored." );
+    OPENGEODE_EXCEPTION(
+        !SpatialDomainChecker< OwnerSegment2D >::is_anchored_in_domain(
+            domain, seg_partial ),
+        "[Test] Segment with lower left extremity outside VOI is not anchored "
+        "by definition." );
+    OPENGEODE_EXCEPTION(
+        !SpatialDomainChecker< OwnerSegment2D >::is_anchored_in_domain(
+            domain, seg_outside ),
+        "[Test] Segment completely outside should not be anchored." );
+
+    OPENGEODE_EXCEPTION(
+        SpatialDomainChecker< OwnerSegment2D >::intersects_domain(
+            domain, seg_inside ),
+        "[Test] Segment inside VOI should intersect." );
+    OPENGEODE_EXCEPTION(
+        SpatialDomainChecker< OwnerSegment2D >::intersects_domain(
+            domain, seg_partial ),
+        "[Test] Segment partially in VOI should intersect." );
+    OPENGEODE_EXCEPTION(
+        !SpatialDomainChecker< OwnerSegment2D >::intersects_domain(
+            domain, seg_outside ),
+        "[Test] Segment completely outside should not intersect." );
+    OPENGEODE_EXCEPTION(
+        SpatialDomainChecker< OwnerSegment2D >::intersects_domain(
+            domain, seg_cross ),
+        "[Test] Segment crossing the voi should intersect." );
+}
+
+int main()
+{
+    try
+    {
+        geode::StochasticLibrary::initialize();
+        geode::Logger::set_level( geode::Logger::LEVEL::debug );
+
+        test_spatial_domain_2D();
+
+        geode::Logger::info( "TEST SpatialDomain SUCCESS" );
+        return 0;
+    }
+    catch( ... )
+    {
+        return geode::geode_lippincott();
+    }
+}
\ No newline at end of file