diff --git a/splashsurf/src/reconstruction.rs b/splashsurf/src/reconstruction.rs
index d67a406..a93f459 100644
--- a/splashsurf/src/reconstruction.rs
+++ b/splashsurf/src/reconstruction.rs
@@ -363,9 +363,13 @@ mod arguments {
             let compact_support_radius = args.particle_radius * 2.0 * args.smoothing_length;
             let cube_size = args.particle_radius * args.cube_size;
 
-            let spatial_decomposition = if !args.octree_decomposition.into_bool() {
-                None
-            } else {
+            let spatial_decomposition = if args.subdomain_grid.into_bool() {
+                Some(splashsurf_lib::SpatialDecomposition::UniformGrid(
+                    splashsurf_lib::GridDecompositionParameters {
+                        subdomain_num_cubes_per_dim: args.subdomain_cubes,
+                    },
+                ))
+            } else if args.octree_decomposition.into_bool() {
                 let subdivision_criterion = if let Some(max_particles) = args.octree_max_particles {
                     splashsurf_lib::SubdivisionCriterion::MaxParticleCount(max_particles)
                 } else {
@@ -388,12 +392,16 @@ mod arguments {
                     }
                 };
 
-                Some(splashsurf_lib::SpatialDecompositionParameters {
-                    subdivision_criterion,
-                    ghost_particle_safety_factor,
-                    enable_stitching,
-                    particle_density_computation,
-                })
+                Some(splashsurf_lib::SpatialDecomposition::Octree(
+                    splashsurf_lib::OctreeDecompositionParameters {
+                        subdivision_criterion,
+                        ghost_particle_safety_factor,
+                        enable_stitching,
+                        particle_density_computation,
+                    },
+                ))
+            } else {
+                None
             };
 
             // Assemble all parameters for the surface reconstruction
@@ -404,10 +412,6 @@ mod arguments {
                 cube_size,
                 iso_surface_threshold: args.surface_threshold,
                 domain_aabb,
-                subdomain_num_cubes_per_dim: args
-                    .subdomain_grid
-                    .into_bool()
-                    .then_some(args.subdomain_cubes),
                 enable_multi_threading: args.parallelize_over_particles.into_bool(),
                 spatial_decomposition,
             };
diff --git a/splashsurf_lib/benches/benches/bench_full.rs b/splashsurf_lib/benches/benches/bench_full.rs
index 8aba88c..3fe6573 100644
--- a/splashsurf_lib/benches/benches/bench_full.rs
+++ b/splashsurf_lib/benches/benches/bench_full.rs
@@ -4,9 +4,9 @@ use splashsurf_lib::io::particles_from_file;
 #[allow(dead_code)]
 use splashsurf_lib::io::vtk_format::write_vtk;
 use splashsurf_lib::{
-    reconstruct_surface, reconstruct_surface_inplace, Parameters,
-    ParticleDensityComputationStrategy, SpatialDecompositionParameters, SubdivisionCriterion,
-    SurfaceReconstruction,
+    reconstruct_surface, reconstruct_surface_inplace, GridDecompositionParameters,
+    OctreeDecompositionParameters, Parameters, ParticleDensityComputationStrategy,
+    SpatialDecomposition, SubdivisionCriterion, SurfaceReconstruction,
 };
 use std::time::Duration;
 
@@ -103,7 +103,6 @@ pub fn surface_reconstruction_dam_break(c: &mut Criterion) {
         iso_surface_threshold: 0.6,
         domain_aabb: None,
         enable_multi_threading: true,
-        subdomain_num_cubes_per_dim: None,
         spatial_decomposition: None,
     };
 
@@ -124,13 +123,15 @@ pub fn surface_reconstruction_dam_break(c: &mut Criterion) {
     group.bench_function("surface_reconstruction_dam_break_par_octree", |b| {
         b.iter(|| {
             let mut parameters = parameters.clone();
-            parameters.spatial_decomposition = Some(SpatialDecompositionParameters {
-                subdivision_criterion: SubdivisionCriterion::MaxParticleCountAuto,
-                ghost_particle_safety_factor: Some(1.0),
-                enable_stitching: false,
-                particle_density_computation:
-                    ParticleDensityComputationStrategy::SynchronizeSubdomains,
-            });
+            parameters.spatial_decomposition = Some(SpatialDecomposition::Octree(
+                OctreeDecompositionParameters {
+                    subdivision_criterion: SubdivisionCriterion::MaxParticleCountAuto,
+                    ghost_particle_safety_factor: Some(1.0),
+                    enable_stitching: false,
+                    particle_density_computation:
+                        ParticleDensityComputationStrategy::SynchronizeSubdomains,
+                },
+            ));
 
             reconstruction =
                 reconstruct_surface::<i64, _>(particle_positions.as_slice(), &parameters).unwrap()
@@ -142,13 +143,15 @@ pub fn surface_reconstruction_dam_break(c: &mut Criterion) {
         |b| {
             b.iter(|| {
                 let mut parameters = parameters.clone();
-                parameters.spatial_decomposition = Some(SpatialDecompositionParameters {
-                    subdivision_criterion: SubdivisionCriterion::MaxParticleCountAuto,
-                    ghost_particle_safety_factor: Some(1.0),
-                    enable_stitching: true,
-                    particle_density_computation:
-                        ParticleDensityComputationStrategy::SynchronizeSubdomains,
-                });
+                parameters.spatial_decomposition = Some(SpatialDecomposition::Octree(
+                    OctreeDecompositionParameters {
+                        subdivision_criterion: SubdivisionCriterion::MaxParticleCountAuto,
+                        ghost_particle_safety_factor: Some(1.0),
+                        enable_stitching: true,
+                        particle_density_computation:
+                            ParticleDensityComputationStrategy::SynchronizeSubdomains,
+                    },
+                ));
 
                 reconstruction =
                     reconstruct_surface::<i64, _>(particle_positions.as_slice(), &parameters)
@@ -160,7 +163,11 @@ pub fn surface_reconstruction_dam_break(c: &mut Criterion) {
     group.bench_function("surface_reconstruction_dam_break_par_grid_64", |b| {
         b.iter(|| {
             let mut parameters = parameters.clone();
-            parameters.subdomain_num_cubes_per_dim = Some(64);
+            parameters.spatial_decomposition = Some(SpatialDecomposition::UniformGrid(
+                GridDecompositionParameters {
+                    subdomain_num_cubes_per_dim: 64,
+                },
+            ));
             reconstruction =
                 reconstruct_surface::<i64, _>(particle_positions.as_slice(), &parameters).unwrap()
         })
@@ -194,7 +201,6 @@ pub fn surface_reconstruction_double_dam_break(c: &mut Criterion) {
         iso_surface_threshold: 0.6,
         domain_aabb: None,
         enable_multi_threading: true,
-        subdomain_num_cubes_per_dim: None,
         spatial_decomposition: None,
     };
 
@@ -215,13 +221,15 @@ pub fn surface_reconstruction_double_dam_break(c: &mut Criterion) {
     group.bench_function("surface_reconstruction_double_dam_break_par_octree", |b| {
         b.iter(|| {
             let mut parameters = parameters.clone();
-            parameters.spatial_decomposition = Some(SpatialDecompositionParameters {
-                subdivision_criterion: SubdivisionCriterion::MaxParticleCountAuto,
-                ghost_particle_safety_factor: Some(1.0),
-                enable_stitching: false,
-                particle_density_computation:
-                    ParticleDensityComputationStrategy::SynchronizeSubdomains,
-            });
+            parameters.spatial_decomposition = Some(SpatialDecomposition::Octree(
+                OctreeDecompositionParameters {
+                    subdivision_criterion: SubdivisionCriterion::MaxParticleCountAuto,
+                    ghost_particle_safety_factor: Some(1.0),
+                    enable_stitching: false,
+                    particle_density_computation:
+                        ParticleDensityComputationStrategy::SynchronizeSubdomains,
+                },
+            ));
 
             reconstruction =
                 reconstruct_surface::<i64, _>(particle_positions.as_slice(), &parameters).unwrap()
@@ -233,13 +241,15 @@ pub fn surface_reconstruction_double_dam_break(c: &mut Criterion) {
         |b| {
             b.iter(|| {
                 let mut parameters = parameters.clone();
-                parameters.spatial_decomposition = Some(SpatialDecompositionParameters {
-                    subdivision_criterion: SubdivisionCriterion::MaxParticleCountAuto,
-                    ghost_particle_safety_factor: Some(1.0),
-                    enable_stitching: true,
-                    particle_density_computation:
-                        ParticleDensityComputationStrategy::SynchronizeSubdomains,
-                });
+                parameters.spatial_decomposition = Some(SpatialDecomposition::Octree(
+                    OctreeDecompositionParameters {
+                        subdivision_criterion: SubdivisionCriterion::MaxParticleCountAuto,
+                        ghost_particle_safety_factor: Some(1.0),
+                        enable_stitching: true,
+                        particle_density_computation:
+                            ParticleDensityComputationStrategy::SynchronizeSubdomains,
+                    },
+                ));
 
                 reconstruction =
                     reconstruct_surface::<i64, _>(particle_positions.as_slice(), &parameters)
@@ -251,7 +261,11 @@ pub fn surface_reconstruction_double_dam_break(c: &mut Criterion) {
     group.bench_function("surface_reconstruction_double_dam_break_par_grid_64", |b| {
         b.iter(|| {
             let mut parameters = parameters.clone();
-            parameters.subdomain_num_cubes_per_dim = Some(64);
+            parameters.spatial_decomposition = Some(SpatialDecomposition::UniformGrid(
+                GridDecompositionParameters {
+                    subdomain_num_cubes_per_dim: 64,
+                },
+            ));
             reconstruction =
                 reconstruct_surface::<i64, _>(particle_positions.as_slice(), &parameters).unwrap()
         })
@@ -285,7 +299,6 @@ pub fn surface_reconstruction_double_dam_break_inplace(c: &mut Criterion) {
         iso_surface_threshold: 0.6,
         domain_aabb: None,
         enable_multi_threading: true,
-        subdomain_num_cubes_per_dim: None,
         spatial_decomposition: None,
     };
 
@@ -315,13 +328,15 @@ pub fn surface_reconstruction_double_dam_break_inplace(c: &mut Criterion) {
         |b| {
             b.iter(|| {
                 let mut parameters = parameters.clone();
-                parameters.spatial_decomposition = Some(SpatialDecompositionParameters {
-                    subdivision_criterion: SubdivisionCriterion::MaxParticleCountAuto,
-                    ghost_particle_safety_factor: Some(1.0),
-                    enable_stitching: false,
-                    particle_density_computation:
-                        ParticleDensityComputationStrategy::SynchronizeSubdomains,
-                });
+                parameters.spatial_decomposition = Some(SpatialDecomposition::Octree(
+                    OctreeDecompositionParameters {
+                        subdivision_criterion: SubdivisionCriterion::MaxParticleCountAuto,
+                        ghost_particle_safety_factor: Some(1.0),
+                        enable_stitching: false,
+                        particle_density_computation:
+                            ParticleDensityComputationStrategy::SynchronizeSubdomains,
+                    },
+                ));
 
                 reconstruct_surface_inplace::<i64, _>(
                     particle_positions.as_slice(),
@@ -338,13 +353,15 @@ pub fn surface_reconstruction_double_dam_break_inplace(c: &mut Criterion) {
         |b| {
             b.iter(|| {
                 let mut parameters = parameters.clone();
-                parameters.spatial_decomposition = Some(SpatialDecompositionParameters {
-                    subdivision_criterion: SubdivisionCriterion::MaxParticleCountAuto,
-                    ghost_particle_safety_factor: Some(1.0),
-                    enable_stitching: true,
-                    particle_density_computation:
-                        ParticleDensityComputationStrategy::SynchronizeSubdomains,
-                });
+                parameters.spatial_decomposition = Some(SpatialDecomposition::Octree(
+                    OctreeDecompositionParameters {
+                        subdivision_criterion: SubdivisionCriterion::MaxParticleCountAuto,
+                        ghost_particle_safety_factor: Some(1.0),
+                        enable_stitching: true,
+                        particle_density_computation:
+                            ParticleDensityComputationStrategy::SynchronizeSubdomains,
+                    },
+                ));
 
                 reconstruct_surface_inplace::<i64, _>(
                     particle_positions.as_slice(),
@@ -361,7 +378,11 @@ pub fn surface_reconstruction_double_dam_break_inplace(c: &mut Criterion) {
         |b| {
             b.iter(|| {
                 let mut parameters = parameters.clone();
-                parameters.subdomain_num_cubes_per_dim = Some(64);
+                parameters.spatial_decomposition = Some(SpatialDecomposition::UniformGrid(
+                    GridDecompositionParameters {
+                        subdomain_num_cubes_per_dim: 64,
+                    },
+                ));
                 reconstruct_surface_inplace::<i64, _>(
                     particle_positions.as_slice(),
                     &parameters,
diff --git a/splashsurf_lib/benches/benches/bench_mesh.rs b/splashsurf_lib/benches/benches/bench_mesh.rs
index 24a864b..baff8d6 100644
--- a/splashsurf_lib/benches/benches/bench_mesh.rs
+++ b/splashsurf_lib/benches/benches/bench_mesh.rs
@@ -2,8 +2,9 @@ use criterion::{criterion_group, Criterion};
 use splashsurf_lib::io::particles_from_file;
 use splashsurf_lib::nalgebra::Vector3;
 use splashsurf_lib::{
-    reconstruct_surface, Parameters, ParticleDensityComputationStrategy,
-    SpatialDecompositionParameters, SubdivisionCriterion, SurfaceReconstruction,
+    reconstruct_surface, OctreeDecompositionParameters, Parameters,
+    ParticleDensityComputationStrategy, SpatialDecomposition, SubdivisionCriterion,
+    SurfaceReconstruction,
 };
 use std::path::Path;
 use std::time::Duration;
@@ -23,13 +24,15 @@ fn reconstruct_particles<P: AsRef<Path>>(particle_file: P) -> SurfaceReconstruct
         iso_surface_threshold: 0.6,
         domain_aabb: None,
         enable_multi_threading: true,
-        subdomain_num_cubes_per_dim: None,
-        spatial_decomposition: Some(SpatialDecompositionParameters {
-            subdivision_criterion: SubdivisionCriterion::MaxParticleCountAuto,
-            ghost_particle_safety_factor: None,
-            enable_stitching: true,
-            particle_density_computation: ParticleDensityComputationStrategy::SynchronizeSubdomains,
-        }),
+        spatial_decomposition: Some(SpatialDecomposition::Octree(
+            OctreeDecompositionParameters {
+                subdivision_criterion: SubdivisionCriterion::MaxParticleCountAuto,
+                ghost_particle_safety_factor: None,
+                enable_stitching: true,
+                particle_density_computation:
+                    ParticleDensityComputationStrategy::SynchronizeSubdomains,
+            },
+        )),
     };
 
     reconstruct_surface::<i64, _>(particle_positions.as_slice(), &parameters).unwrap()
diff --git a/splashsurf_lib/benches/benches/bench_subdomain_grid.rs b/splashsurf_lib/benches/benches/bench_subdomain_grid.rs
index 4adba28..e14ec3f 100644
--- a/splashsurf_lib/benches/benches/bench_subdomain_grid.rs
+++ b/splashsurf_lib/benches/benches/bench_subdomain_grid.rs
@@ -1,7 +1,10 @@
 use criterion::{criterion_group, Criterion, SamplingMode};
 use nalgebra::Vector3;
 use splashsurf_lib::io::particles_from_file;
-use splashsurf_lib::{reconstruct_surface, Parameters, SurfaceReconstruction};
+use splashsurf_lib::{
+    reconstruct_surface, GridDecompositionParameters, Parameters, SpatialDecomposition,
+    SurfaceReconstruction,
+};
 use std::time::Duration;
 
 fn parameters_canyon() -> Parameters<f32> {
@@ -17,8 +20,11 @@ fn parameters_canyon() -> Parameters<f32> {
         iso_surface_threshold: 0.6,
         domain_aabb: None,
         enable_multi_threading: true,
-        subdomain_num_cubes_per_dim: Some(32),
-        spatial_decomposition: None,
+        spatial_decomposition: Some(SpatialDecomposition::UniformGrid(
+            GridDecompositionParameters {
+                subdomain_num_cubes_per_dim: 32,
+            },
+        )),
     };
 
     parameters
@@ -40,7 +46,11 @@ pub fn grid_canyon(c: &mut Criterion) {
         b.iter(|| {
             let mut parameters = parameters.clone();
             parameters.cube_size = 1.5 * parameters.particle_radius;
-            parameters.subdomain_num_cubes_per_dim = Some(32);
+            parameters.spatial_decomposition = Some(SpatialDecomposition::UniformGrid(
+                GridDecompositionParameters {
+                    subdomain_num_cubes_per_dim: 32,
+                },
+            ));
             reconstruction =
                 reconstruct_surface::<i64, _>(particle_positions.as_slice(), &parameters).unwrap()
         })
@@ -50,7 +60,11 @@ pub fn grid_canyon(c: &mut Criterion) {
         b.iter(|| {
             let mut parameters = parameters.clone();
             parameters.cube_size = 1.0 * parameters.particle_radius;
-            parameters.subdomain_num_cubes_per_dim = Some(48);
+            parameters.spatial_decomposition = Some(SpatialDecomposition::UniformGrid(
+                GridDecompositionParameters {
+                    subdomain_num_cubes_per_dim: 48,
+                },
+            ));
             reconstruction =
                 reconstruct_surface::<i64, _>(particle_positions.as_slice(), &parameters).unwrap()
         })
@@ -60,7 +74,11 @@ pub fn grid_canyon(c: &mut Criterion) {
         b.iter(|| {
             let mut parameters = parameters.clone();
             parameters.cube_size = 0.75 * parameters.particle_radius;
-            parameters.subdomain_num_cubes_per_dim = Some(48);
+            parameters.spatial_decomposition = Some(SpatialDecomposition::UniformGrid(
+                GridDecompositionParameters {
+                    subdomain_num_cubes_per_dim: 64,
+                },
+            ));
             reconstruction =
                 reconstruct_surface::<i64, _>(particle_positions.as_slice(), &parameters).unwrap()
         })
@@ -89,7 +107,11 @@ pub fn grid_optimal_num_cubes_canyon(c: &mut Criterion) {
             group.bench_function(format!("subdomain_num_cubes_{}", num_cubes), |b| {
                 b.iter(|| {
                     let mut parameters = parameters.clone();
-                    parameters.subdomain_num_cubes_per_dim = Some(num_cubes);
+                    parameters.spatial_decomposition = Some(SpatialDecomposition::UniformGrid(
+                        GridDecompositionParameters {
+                            subdomain_num_cubes_per_dim: num_cubes,
+                        },
+                    ));
                     reconstruction =
                         reconstruct_surface::<i64, _>(particle_positions.as_slice(), &parameters)
                             .unwrap()
diff --git a/splashsurf_lib/src/dense_subdomains.rs b/splashsurf_lib/src/dense_subdomains.rs
index 5ce3c7d..635a648 100644
--- a/splashsurf_lib/src/dense_subdomains.rs
+++ b/splashsurf_lib/src/dense_subdomains.rs
@@ -1,4 +1,4 @@
-use anyhow::Context;
+use anyhow::{anyhow, Context};
 use arrayvec::ArrayVec;
 use itertools::Itertools;
 use log::{info, trace};
@@ -21,10 +21,13 @@ use crate::neighborhood_search::{
 };
 use crate::uniform_grid::{EdgeIndex, UniformCartesianCubeGrid3d};
 use crate::{
-    new_map, new_parallel_map, profile, Aabb3d, MapType, Parameters, SurfaceReconstruction,
+    new_map, new_parallel_map, profile, Aabb3d, MapType, Parameters, SpatialDecomposition,
+    SurfaceReconstruction,
 };
 use crate::{Index, Real};
 
+// TODO: Implement single-threaded processing
+
 type GlobalIndex = u64;
 
 /// Converts any literal or expression to the Index type I (panics if value does not fit)
@@ -86,8 +89,12 @@ pub(crate) fn initialize_parameters<'a, I: Index, R: Real>(
 ) -> Result<ParametersSubdomainGrid<I, R>, anyhow::Error> {
     let chunk_size = 500;
 
+    let Some(SpatialDecomposition::UniformGrid(grid_parameters)) = &parameters.spatial_decomposition else {
+        return Err(anyhow!("spatial decomposition parameters for uniform grid are missing"))
+    };
+
     // A subdomain will be a cube consisting of this number of MC cubes along each coordinate axis
-    let subdomain_cubes_in = parameters.subdomain_num_cubes_per_dim.unwrap_or(64);
+    let subdomain_cubes_in = grid_parameters.subdomain_num_cubes_per_dim;
     let subdomain_cubes = I::from_u32(subdomain_cubes_in)
         .expect("number of subdomain cubes has to fit in index type");
     let subdomain_cubes_global = GlobalIndex::from_u32(subdomain_cubes_in)
diff --git a/splashsurf_lib/src/lib.rs b/splashsurf_lib/src/lib.rs
index bc8149b..5d3fc21 100644
--- a/splashsurf_lib/src/lib.rs
+++ b/splashsurf_lib/src/lib.rs
@@ -114,9 +114,53 @@ fn new_parallel_map<K: Eq + Hash, V>() -> ParallelMapType<K, V> {
     ParallelMapType::with_hasher(HashState::default())
 }
 
-/// Parameters for the spatial decomposition
+/// Approach used for spatial decomposition of the surface reconstruction and its parameters
 #[derive(Clone, Debug)]
-pub struct SpatialDecompositionParameters<R: Real> {
+pub enum SpatialDecomposition<R: Real> {
+    /// Use a uniform grid of subdomains with contiguous (dense) marching cubes grids per subdomain
+    ///
+    /// Only subdomains containing at least one particle will be processed.
+    /// The small contiguous grid per subdomain make this approach very cache efficient.
+    UniformGrid(GridDecompositionParameters),
+    /// Use an octree for decomposition with hashing based (sparse) marching cubes grids per subdomain
+    Octree(OctreeDecompositionParameters<R>),
+}
+
+/// Default parameters for the spatial decomposition use the uniform grid based decomposition approach
+impl<R: Real> Default for SpatialDecomposition<R> {
+    fn default() -> Self {
+        Self::UniformGrid(GridDecompositionParameters::default())
+    }
+}
+
+impl<R: Real> SpatialDecomposition<R> {
+    /// Tries to convert the parameters from one [`Real`] type to another [`Real`] type, returns `None` if conversion fails
+    pub fn try_convert<T: Real>(&self) -> Option<SpatialDecomposition<T>> {
+        match &self {
+            Self::UniformGrid(g) => Some(SpatialDecomposition::UniformGrid(g.clone())),
+            Self::Octree(o) => o.try_convert().map(|o| SpatialDecomposition::Octree(o)),
+        }
+    }
+}
+
+/// Parameters for the uniform grid-based spatial decomposition
+#[derive(Clone, Debug)]
+pub struct GridDecompositionParameters {
+    /// Each uniform subdomain will be a cube consisting of this number of MC cube cells along each coordinate axis
+    pub subdomain_num_cubes_per_dim: u32,
+}
+
+impl Default for GridDecompositionParameters {
+    fn default() -> Self {
+        Self {
+            subdomain_num_cubes_per_dim: 64,
+        }
+    }
+}
+
+/// Parameters for the octree-based spatial decomposition
+#[derive(Clone, Debug)]
+pub struct OctreeDecompositionParameters<R: Real> {
     /// Criterion used for subdivision of the octree cells
     pub subdivision_criterion: SubdivisionCriterion,
     /// Safety factor applied to the kernel radius when it's used as a margin to collect ghost particles in the leaf nodes
@@ -169,10 +213,21 @@ pub enum ParticleDensityComputationStrategy {
     IndependentSubdomains,
 }
 
-impl<R: Real> SpatialDecompositionParameters<R> {
+impl<R: Real> Default for OctreeDecompositionParameters<R> {
+    fn default() -> Self {
+        Self {
+            subdivision_criterion: SubdivisionCriterion::MaxParticleCountAuto,
+            ghost_particle_safety_factor: None,
+            enable_stitching: true,
+            particle_density_computation: ParticleDensityComputationStrategy::SynchronizeSubdomains,
+        }
+    }
+}
+
+impl<R: Real> OctreeDecompositionParameters<R> {
     /// Tries to convert the parameters from one [`Real`] type to another [`Real`] type, returns `None` if conversion fails
-    pub fn try_convert<T: Real>(&self) -> Option<SpatialDecompositionParameters<T>> {
-        Some(SpatialDecompositionParameters {
+    pub fn try_convert<T: Real>(&self) -> Option<OctreeDecompositionParameters<T>> {
+        Some(OctreeDecompositionParameters {
             subdivision_criterion: self.subdivision_criterion.clone(),
             ghost_particle_safety_factor: map_option!(
                 &self.ghost_particle_safety_factor,
@@ -197,16 +252,18 @@ pub struct Parameters<R: Real> {
     pub cube_size: R,
     /// Density threshold value to distinguish between the inside (above threshold) and outside (below threshold) of the fluid
     pub iso_surface_threshold: R,
-    /// Manually restrict the domain to the surface reconstruction.
+    /// Bounding box of particles to reconstruct
+    ///
+    /// All particles outside of this domain will be filtered out before the reconstruction.
+    /// The surface reconstruction always results in a closed mesh around the particles.
+    /// The final mesh can extend beyond this AABB due to the smoothing of the kernel.
     /// If not provided, the smallest AABB enclosing all particles is computed instead.
     pub domain_aabb: Option<Aabb3d<R>>,
     /// Whether to allow multi threading within the surface reconstruction procedure
     pub enable_multi_threading: bool,
-    /// Each subdomain will be a cube consisting of this number of MC cube cells along each coordinate axis
-    pub subdomain_num_cubes_per_dim: Option<u32>,
-    /// Parameters for the spatial decomposition (octree subdivision) of the particles.
-    /// If not provided, no octree is generated and a global approach is used instead.
-    pub spatial_decomposition: Option<SpatialDecompositionParameters<R>>,
+    /// Parameters for the spatial decomposition of the surface reconstruction
+    /// If not provided, no spatial decomposition is performed and a global approach is used instead.
+    pub spatial_decomposition: Option<SpatialDecomposition<R>>,
 }
 
 impl<R: Real> Parameters<R> {
@@ -220,7 +277,6 @@ impl<R: Real> Parameters<R> {
             iso_surface_threshold: self.iso_surface_threshold.try_convert()?,
             domain_aabb: map_option!(&self.domain_aabb, aabb => aabb.try_convert()?),
             enable_multi_threading: self.enable_multi_threading,
-            subdomain_num_cubes_per_dim: self.subdomain_num_cubes_per_dim,
             spatial_decomposition: map_option!(&self.spatial_decomposition, sd => sd.try_convert()?),
         })
     }
@@ -365,29 +421,37 @@ pub fn reconstruct_surface_inplace<'a, I: Index, R: Real>(
 
     output_surface.grid.log_grid_info();
 
-    if parameters.subdomain_num_cubes_per_dim.is_some() {
-        reconstruction::reconstruct_surface_subdomain_grid::<I, R>(
-            particle_positions,
-            parameters,
-            output_surface,
-        )?;
-    } else if parameters.spatial_decomposition.is_some() {
-        reconstruction_octree::reconstruct_surface_domain_decomposition(
-            particle_positions,
-            parameters,
-            output_surface,
-        )?;
-    } else {
-        reconstruction_octree::reconstruct_surface_global(
+    match &parameters.spatial_decomposition {
+        Some(SpatialDecomposition::UniformGrid(_)) => {
+            reconstruction::reconstruct_surface_subdomain_grid::<I, R>(
+                particle_positions,
+                parameters,
+                output_surface,
+            )?
+        }
+        Some(SpatialDecomposition::Octree(_)) => {
+            reconstruction_octree::reconstruct_surface_domain_decomposition(
+                particle_positions,
+                parameters,
+                output_surface,
+            )?
+        }
+        None => reconstruction_octree::reconstruct_surface_global(
             particle_positions,
             parameters,
             output_surface,
-        )?;
+        )?,
     }
 
     Ok(())
 }
 
+/*
+fn filter_particles<I: Index, R: Real>(particle_positions: &[Vector3<R>], particle_aabb: &Aabb3d<R>, enable_multi_threading: bool) -> Vec<Vector3<R>> {
+    use rayon::prelude::*;
+    let is_inside = particle_positions.par_iter().map(|p| particle_aabb.contains_point(p)).collect::<Vec<_>>();
+}*/
+
 /// Constructs the background grid for marching cubes based on the parameters supplied to the surface reconstruction
 pub fn grid_for_reconstruction<I: Index, R: Real>(
     particle_positions: &[Vector3<R>],
diff --git a/splashsurf_lib/src/reconstruction_octree.rs b/splashsurf_lib/src/reconstruction_octree.rs
index 4bfe66d..f97a58d 100644
--- a/splashsurf_lib/src/reconstruction_octree.rs
+++ b/splashsurf_lib/src/reconstruction_octree.rs
@@ -7,9 +7,9 @@ use crate::octree::{NodeData, Octree, OctreeNode};
 use crate::uniform_grid::{OwningSubdomainGrid, Subdomain, UniformGrid};
 use crate::workspace::LocalReconstructionWorkspace;
 use crate::{
-    density_map, marching_cubes, neighborhood_search, new_map, profile, utils, Index, Parameters,
-    ParticleDensityComputationStrategy, Real, ReconstructionError, SpatialDecompositionParameters,
-    SurfaceReconstruction,
+    density_map, marching_cubes, neighborhood_search, new_map, profile, utils, Index,
+    OctreeDecompositionParameters, Parameters, ParticleDensityComputationStrategy, Real,
+    ReconstructionError, SpatialDecomposition, SurfaceReconstruction,
 };
 use log::{debug, info, trace};
 use nalgebra::Vector3;
@@ -82,7 +82,7 @@ struct OctreeBasedSurfaceReconstruction<I: Index, R: Real> {
     /// General parameters for the surface reconstruction
     parameters: Parameters<R>,
     /// Spatial decomposition specific parameters extracted from the general parameters
-    spatial_decomposition: SpatialDecompositionParameters<R>,
+    spatial_decomposition: OctreeDecompositionParameters<R>,
     /// The implicit global grid for the entire surface reconstruction
     grid: UniformGrid<I, R>,
     /// Octree containing the individual particle subdomains, built during the construction of this helper type
@@ -100,8 +100,15 @@ impl<I: Index, R: Real> OctreeBasedSurfaceReconstruction<I, R> {
         // The grid was already generated by the calling public function
         let grid = output_surface.grid.clone();
 
+        let Some(SpatialDecomposition::Octree(decomposition_parameters)) = &parameters.spatial_decomposition else {
+            // TODO: Use default values instead?
+
+            // If there are no decomposition parameters, we cannot construct an octree.
+            return None;
+        };
+
         // Construct the octree
-        let octree = if let Some(decomposition_parameters) = &parameters.spatial_decomposition {
+        let octree = {
             let margin_factor = decomposition_parameters
                 .ghost_particle_safety_factor
                 .unwrap_or(R::one());
@@ -114,11 +121,6 @@ impl<I: Index, R: Real> OctreeBasedSurfaceReconstruction<I, R> {
                 parameters.enable_multi_threading,
                 decomposition_parameters.enable_stitching,
             )
-        } else {
-            // TODO: Use default values instead?
-
-            // If there are no decomposition parameters, we cannot construct an octree.
-            return None;
         };
 
         // Disable all multi-threading in sub-tasks for now (instead, entire sub-tasks are processed in parallel)
@@ -130,7 +132,7 @@ impl<I: Index, R: Real> OctreeBasedSurfaceReconstruction<I, R> {
 
         Some(Self {
             octree,
-            spatial_decomposition: parameters.spatial_decomposition.as_ref().unwrap().clone(),
+            spatial_decomposition: decomposition_parameters.clone(),
             grid,
             parameters,
         })
diff --git a/splashsurf_lib/tests/integration_tests/test_full.rs b/splashsurf_lib/tests/integration_tests/test_full.rs
index 524a36e..d8c570b 100644
--- a/splashsurf_lib/tests/integration_tests/test_full.rs
+++ b/splashsurf_lib/tests/integration_tests/test_full.rs
@@ -3,8 +3,9 @@ use splashsurf_lib::io::particles_from_file;
 use splashsurf_lib::io::vtk_format::write_vtk;
 use splashsurf_lib::marching_cubes::check_mesh_consistency;
 use splashsurf_lib::{
-    reconstruct_surface, Aabb3d, Parameters, ParticleDensityComputationStrategy, Real,
-    SpatialDecompositionParameters, SubdivisionCriterion,
+    reconstruct_surface, Aabb3d, GridDecompositionParameters, OctreeDecompositionParameters,
+    Parameters, ParticleDensityComputationStrategy, Real, SpatialDecomposition,
+    SubdivisionCriterion,
 };
 use std::path::Path;
 
@@ -37,32 +38,39 @@ fn params_with_aabb<R: Real>(
         iso_surface_threshold,
         domain_aabb,
         enable_multi_threading: false,
-        subdomain_num_cubes_per_dim: None,
         spatial_decomposition: None,
     };
 
     match strategy {
         Strategy::Global => {}
         Strategy::Octree => {
-            parameters.spatial_decomposition = Some(SpatialDecompositionParameters {
-                subdivision_criterion: SubdivisionCriterion::MaxParticleCountAuto,
-                ghost_particle_safety_factor: Some(R::one()),
-                enable_stitching: false,
-                particle_density_computation:
-                    ParticleDensityComputationStrategy::SynchronizeSubdomains,
-            });
+            parameters.spatial_decomposition = Some(SpatialDecomposition::Octree(
+                OctreeDecompositionParameters {
+                    subdivision_criterion: SubdivisionCriterion::MaxParticleCountAuto,
+                    ghost_particle_safety_factor: Some(R::one()),
+                    enable_stitching: false,
+                    particle_density_computation:
+                        ParticleDensityComputationStrategy::SynchronizeSubdomains,
+                },
+            ));
         }
         Strategy::OctreeStitching => {
-            parameters.spatial_decomposition = Some(SpatialDecompositionParameters {
-                subdivision_criterion: SubdivisionCriterion::MaxParticleCountAuto,
-                ghost_particle_safety_factor: Some(R::one() + R::one()),
-                enable_stitching: true,
-                particle_density_computation:
-                    ParticleDensityComputationStrategy::SynchronizeSubdomains,
-            });
+            parameters.spatial_decomposition = Some(SpatialDecomposition::Octree(
+                OctreeDecompositionParameters {
+                    subdivision_criterion: SubdivisionCriterion::MaxParticleCountAuto,
+                    ghost_particle_safety_factor: Some(R::one() + R::one()),
+                    enable_stitching: true,
+                    particle_density_computation:
+                        ParticleDensityComputationStrategy::SynchronizeSubdomains,
+                },
+            ));
         }
         Strategy::SubdomainGrid => {
-            parameters.subdomain_num_cubes_per_dim = Some(64);
+            parameters.spatial_decomposition = Some(SpatialDecomposition::UniformGrid(
+                GridDecompositionParameters {
+                    subdomain_num_cubes_per_dim: 64,
+                },
+            ))
         }
     }
 
@@ -101,11 +109,11 @@ fn default_params<R: Real>() -> Parameters<R> {
 }
 
 fn test_for_boundary<R: Real>(params: &Parameters<R>) -> bool {
-    params
-        .spatial_decomposition
-        .as_ref()
-        .map(|s| s.enable_stitching)
-        .unwrap_or(true)
+    match &params.spatial_decomposition {
+        Some(SpatialDecomposition::Octree(p)) => p.enable_stitching,
+        Some(SpatialDecomposition::UniformGrid(_)) => true,
+        None => true,
+    }
 }
 
 macro_rules! generate_test {