diff --git a/CHANGELOG.md b/CHANGELOG.md
index 507f40ec..5657af8c 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,11 +1,42 @@
 ## Unreleased
 
+### Breaking changes
+
+- The voxel query functions are now generic over the voxel storage instead of taking `&Voxels`:
+  `contact_manifolds_voxels_shape`, `contact_manifolds_voxels_ball`,
+  `contact_manifolds_voxels_composite_shape`, `contact_manifolds_voxels_voxels`,
+  `intersection_test_voxels_shape`, `intersection_test_shape_voxels`, `cast_shapes_voxels_shape`,
+  `cast_shapes_shape_voxels`, `cast_shapes_nonlinear_voxels_shape`,
+  `cast_shapes_nonlinear_shape_voxels`. Calls that pass a `&Voxels` keep working unchanged;
+  callers naming the functions with explicit turbofish generics gain one extra type parameter.
+- Reading a voxel's type, state, center, or grid coordinates from the item yielded by
+  `Voxels::voxels`, `Voxels::voxels_in_range`, and `Voxels::voxels_intersecting_local_aabb` should
+  go through the new `QueriedVoxel` trait methods (`voxel_type()`, `voxel_state()`, `center()`,
+  `grid_coords()`, `linear_id()`) for code that must also work with custom storages. The public
+  fields of `VoxelData` remain available.
+
 ### Added
 
 - `CompoundFlags::FIX_INTERNAL_EDGES` makes a `Compound` treat the edges (2D) or faces (3D) its parts share as
   interior to the union, so a body sliding across the cut between two parts of a convex
   decomposition no longer catches on it. `Compound::PartNormalConstraints` is now
   `CompoundPseudoNormals`, matching what `TriMesh` and `Polyline` already provide.
+- `VoxelQuery` trait, an abstraction over the storage of a shape made of axis-aligned, uniformly
+  sized voxels. Implementing it for a custom sparse structure (chunked grid, octree, VDB-like tree)
+  lets Parry's voxel collision algorithms (contact manifolds, intersection tests, linear and
+  nonlinear shape-casting) run directly on that structure without copying it into a `Voxels` shape. Implementors provide `voxel_size`, `domain`, and `voxels_in_range`; grid helpers
+  such as `voxel_at_point`, `voxel_center`, `voxel_aabb`, `voxel_range_intersecting_local_aabb`,
+  `align_aabb_to_grid`, and `local_aabb` have default implementations.
+- `QueriedVoxel` trait describing the per-voxel view handed out by a `VoxelQuery` storage. Views
+  may borrow from their storage so that `voxel_state()` can be computed lazily from local context,
+  while `voxel_type()` stays cheap for bulk iteration.
+- `Voxels` implements `VoxelQuery` with `VoxelData` as its voxel view, and `VoxelData` implements
+  `QueriedVoxel`.
+- `VoxelState::with_filled_neighbors(AxisMask)` builds the state of a non-empty voxel from the set
+  of its filled axis-aligned neighbors, for custom storages that only track occupancy.
+- `Default` implementations for `VoxelType` (`Empty`) and `VoxelState` (`EMPTY`).
+- `contact_manifolds_voxels_ball` is now re-exported from `parry::query`, alongside the other
+  voxel contact-manifold functions.
 
 ## 0.30.2
 
diff --git a/src/mass_properties/mass_properties_voxels.rs b/src/mass_properties/mass_properties_voxels.rs
index e08b907f..83b98bee 100644
--- a/src/mass_properties/mass_properties_voxels.rs
+++ b/src/mass_properties/mass_properties_voxels.rs
@@ -2,7 +2,7 @@ use crate::mass_properties::MassProperties;
 #[cfg(feature = "dim3")]
 use crate::math::Matrix;
 use crate::math::{Real, Vector};
-use crate::shape::Voxels;
+use crate::shape::{QueriedVoxel, VoxelType, Voxels};
 
 impl MassProperties {
     /// Computes the mass properties of a voxel grid.
@@ -163,7 +163,9 @@ impl MassProperties {
     ///
     /// - Only non-empty voxels contribute to mass
     /// - Empty voxels are ignored (zero mass, no inertia)
-    /// - The voxel state is checked using `vox.state.is_empty()`
+    /// - A voxel is considered empty if its
+    ///   [`QueriedVoxel::voxel_type`] is
+    ///   [`VoxelType::Empty`]
     ///
     /// # See Also
     ///
@@ -181,8 +183,8 @@ impl MassProperties {
         let block_ref_mprops = MassProperties::from_cuboid(density, voxels.voxel_size() / 2.0);
 
         for vox in voxels.voxels() {
-            if !vox.state.is_empty() {
-                com += vox.center;
+            if vox.voxel_type() != VoxelType::Empty {
+                com += vox.center();
                 num_not_empty += 1;
             }
         }
@@ -190,9 +192,9 @@ impl MassProperties {
         com /= num_not_empty as Real;
 
         for vox in voxels.voxels() {
-            if !vox.state.is_empty() {
+            if vox.voxel_type() != VoxelType::Empty {
                 angular_inertia +=
-                    block_ref_mprops.construct_shifted_inertia_matrix(vox.center - com);
+                    block_ref_mprops.construct_shifted_inertia_matrix(vox.center() - com);
             }
         }
 
diff --git a/src/query/contact_manifolds/contact_manifolds_voxels_ball.rs b/src/query/contact_manifolds/contact_manifolds_voxels_ball.rs
index fc16604c..291dfbb3 100644
--- a/src/query/contact_manifolds/contact_manifolds_voxels_ball.rs
+++ b/src/query/contact_manifolds/contact_manifolds_voxels_ball.rs
@@ -2,7 +2,8 @@ use crate::bounding_volume::BoundingVolume;
 use crate::math::{Pose, Real, Vector, VectorExt};
 use crate::query::{ContactManifold, PointQuery, TrackedContact};
 use crate::shape::{
-    Ball, Cuboid, OctantPattern, PackedFeatureId, Shape, VoxelState, VoxelType, Voxels,
+    Ball, Cuboid, OctantPattern, PackedFeatureId, QueriedVoxel, Shape, VoxelQuery, VoxelState,
+    VoxelType,
 };
 use alloc::vec::Vec;
 
@@ -31,10 +32,12 @@ pub fn contact_manifolds_voxels_ball_shapes<ManifoldData, ContactData>(
     }
 }
 
-/// Computes the contact manifold between a convex shape and a ball.
-pub fn contact_manifolds_voxels_ball<'a, ManifoldData, ContactData>(
+/// Computes the contact manifold between a voxels shape and a ball.
+///
+/// The voxels shape can be any voxel storage implementing [`VoxelQuery`].
+pub fn contact_manifolds_voxels_ball<'a, ManifoldData, ContactData, V>(
     pos12: &Pose,
-    voxels1: &'a Voxels,
+    voxels1: &'a V,
     ball2: &'a Ball,
     prediction: Real,
     manifolds: &mut Vec<ContactManifold<ManifoldData, ContactData>>,
@@ -42,6 +45,7 @@ pub fn contact_manifolds_voxels_ball<'a, ManifoldData, ContactData>(
 ) where
     ManifoldData: Default,
     ContactData: Default + Copy,
+    V: ?Sized + VoxelQuery,
 {
     // TODO: don’t generate one manifold per voxel.
     manifolds.clear();
@@ -55,7 +59,7 @@ pub fn contact_manifolds_voxels_ball<'a, ManifoldData, ContactData>(
     let aabb2 = ball2.aabb(pos12).loosened(prediction / 2.0);
     if let Some(aabb_intersection) = aabb1.intersection(&aabb2) {
         for vox1 in voxels1.voxels_intersecting_local_aabb(&aabb_intersection) {
-            match vox1.state.voxel_type() {
+            match vox1.voxel_type() {
                 #[cfg(feature = "dim2")]
                 VoxelType::Vertex | VoxelType::Face => { /* Ok */ }
                 #[cfg(feature = "dim3")]
@@ -65,9 +69,9 @@ pub fn contact_manifolds_voxels_ball<'a, ManifoldData, ContactData>(
 
             detect_hit_voxel_ball(
                 *pos12,
-                vox1.center,
+                vox1.center(),
                 radius1,
-                vox1.state,
+                vox1.voxel_state(),
                 center2,
                 radius2,
                 prediction,
diff --git a/src/query/contact_manifolds/contact_manifolds_voxels_composite_shape.rs b/src/query/contact_manifolds/contact_manifolds_voxels_composite_shape.rs
index af948a09..ca74a58b 100644
--- a/src/query/contact_manifolds/contact_manifolds_voxels_composite_shape.rs
+++ b/src/query/contact_manifolds/contact_manifolds_voxels_composite_shape.rs
@@ -7,7 +7,9 @@ use crate::query::{
     ContactManifold, ContactManifoldsWorkspace, PersistentQueryDispatcher, PointQuery,
     TypedWorkspaceData, WorkspaceData,
 };
-use crate::shape::{CompositeShape, Cuboid, Shape, SupportMap, VoxelType, Voxels};
+use crate::shape::{
+    CompositeShape, Cuboid, QueriedVoxel, Shape, SupportMap, VoxelQuery, VoxelType,
+};
 use crate::utils::hashmap::Entry;
 use crate::utils::PoseOpt;
 use alloc::{boxed::Box, vec::Vec};
@@ -55,10 +57,12 @@ pub fn contact_manifolds_voxels_composite_shape_shapes<ManifoldData, ContactData
 }
 
 /// Computes the contact manifold between voxels and a composite shape.
-pub fn contact_manifolds_voxels_composite_shape<ManifoldData, ContactData>(
+///
+/// The voxels shape can be any voxel storage implementing [`VoxelQuery`].
+pub fn contact_manifolds_voxels_composite_shape<ManifoldData, ContactData, V>(
     dispatcher: &dyn PersistentQueryDispatcher<ManifoldData, ContactData>,
     pos12: &Pose,
-    voxels1: &Voxels,
+    voxels1: &V,
     shape2: &dyn CompositeShape,
     prediction: Real,
     manifolds: &mut Vec<ContactManifold<ManifoldData, ContactData>>,
@@ -67,6 +71,7 @@ pub fn contact_manifolds_voxels_composite_shape<ManifoldData, ContactData>(
 ) where
     ManifoldData: Default + Clone,
     ContactData: Default + Copy,
+    V: ?Sized + VoxelQuery,
 {
     VoxelsShapeContactManifoldsWorkspace::<3>::ensure_exists(workspace);
     let workspace: &mut VoxelsShapeContactManifoldsWorkspace<3> =
@@ -89,7 +94,7 @@ pub fn contact_manifolds_voxels_composite_shape<ManifoldData, ContactData>(
 
     if let Some(intersection_aabb1) = aabb1.intersection(&aabb2_1) {
         for vox1 in voxels1.voxels_intersecting_local_aabb(&intersection_aabb1) {
-            let vox_type1 = vox1.state.voxel_type();
+            let vox_type1 = vox1.voxel_type();
 
             // TODO: would be nice to have a strategy to handle interior voxels for depenetration.
             if vox_type1 == VoxelType::Empty || vox_type1 == VoxelType::Interior {
@@ -135,7 +140,7 @@ pub fn contact_manifolds_voxels_composite_shape<ManifoldData, ContactData>(
                                 timestamp: new_timestamp,
                             };
 
-                            let vox_id = vox1.linear_id.flat_id() as u32;
+                            let vox_id = vox1.linear_id();
                             let (id1, id2) = if flipped {
                                 (leaf2, vox_id)
                             } else {
@@ -237,7 +242,8 @@ pub fn contact_manifolds_voxels_composite_shape<ManifoldData, ContactData>(
                             // interior of the infinitely expanded canonical shape by checking if
                             // the opposite normal had led to a better vector.
                             let cuboid1 = Cuboid::new(radius1);
-                            let sp1 = cuboid1.local_support_point(-penetration_dir1) + vox1.center;
+                            let sp1 =
+                                cuboid1.local_support_point(-penetration_dir1) + vox1.center();
                             let sm2 = part_shape2
                                 .as_support_map()
                                 .expect("Unsupported collision pair.");
@@ -254,9 +260,9 @@ pub fn contact_manifolds_voxels_composite_shape<ManifoldData, ContactData>(
                         }
 
                         let pt_in_voxel_space = if flipped {
-                            manifold.subshape_pos2().transform_point(pt.local_p2) - vox1.center
+                            manifold.subshape_pos2().transform_point(pt.local_p2) - vox1.center()
                         } else {
-                            manifold.subshape_pos1().transform_point(pt.local_p1) - vox1.center
+                            manifold.subshape_pos1().transform_point(pt.local_p1) - vox1.center()
                         };
                         sub_detector.selected_contacts |=
                             (test_voxel.contains_local_point(pt_in_voxel_space) as u32) << i;
diff --git a/src/query/contact_manifolds/contact_manifolds_voxels_shape.rs b/src/query/contact_manifolds/contact_manifolds_voxels_shape.rs
index f0b50e6d..71a6d933 100644
--- a/src/query/contact_manifolds/contact_manifolds_voxels_shape.rs
+++ b/src/query/contact_manifolds/contact_manifolds_voxels_shape.rs
@@ -4,7 +4,7 @@ use crate::query::{
     ContactManifold, ContactManifoldsWorkspace, PersistentQueryDispatcher, PointQuery,
     TypedWorkspaceData, WorkspaceData,
 };
-use crate::shape::{AxisMask, Cuboid, Shape, SupportMap, VoxelData, VoxelType, Voxels};
+use crate::shape::{AxisMask, Cuboid, QueriedVoxel, Shape, SupportMap, VoxelQuery, VoxelType};
 use crate::utils::hashmap::{Entry, HashMap};
 use crate::utils::PoseOpt;
 use alloc::{boxed::Box, vec::Vec};
@@ -113,10 +113,12 @@ pub fn contact_manifolds_voxels_shape_shapes<ManifoldData, ContactData>(
 }
 
 /// Computes the contact manifold between a convex shape and a voxels shape.
-pub fn contact_manifolds_voxels_shape<ManifoldData, ContactData>(
+///
+/// The voxels shape can be any voxel storage implementing [`VoxelQuery`].
+pub fn contact_manifolds_voxels_shape<ManifoldData, ContactData, V>(
     dispatcher: &dyn PersistentQueryDispatcher<ManifoldData, ContactData>,
     pos12: &Pose,
-    voxels1: &Voxels,
+    voxels1: &V,
     shape2: &dyn Shape,
     prediction: Real,
     manifolds: &mut Vec<ContactManifold<ManifoldData, ContactData>>,
@@ -125,6 +127,7 @@ pub fn contact_manifolds_voxels_shape<ManifoldData, ContactData>(
 ) where
     ManifoldData: Default + Clone,
     ContactData: Default + Copy,
+    V: ?Sized + VoxelQuery,
 {
     VoxelsShapeContactManifoldsWorkspace::<2>::ensure_exists(workspace);
     let workspace: &mut VoxelsShapeContactManifoldsWorkspace<2> =
@@ -148,7 +151,7 @@ pub fn contact_manifolds_voxels_shape<ManifoldData, ContactData>(
 
     if let Some(intersection_aabb1) = aabb1.intersection(&aabb2_1) {
         for vox1 in voxels1.voxels_intersecting_local_aabb(&intersection_aabb1) {
-            let vox_type1 = vox1.state.voxel_type();
+            let vox_type1 = vox1.voxel_type();
 
             // TODO: would be nice to have a strategy to handle interior voxels for depenetration.
             if vox_type1 == VoxelType::Empty || vox_type1 == VoxelType::Interior {
@@ -185,7 +188,7 @@ pub fn contact_manifolds_voxels_shape<ManifoldData, ContactData>(
                             timestamp: new_timestamp,
                         };
 
-                        let vid = vox1.linear_id.flat_id() as u32;
+                        let vid = vox1.linear_id();
                         let (id1, id2) = if flipped { (0, vid) } else { (vid, 0) };
                         manifolds.push(ContactManifold::with_data(
                             id1,
@@ -282,7 +285,7 @@ pub fn contact_manifolds_voxels_shape<ManifoldData, ContactData>(
                     // interior of the infinitely expanded canonical shape by checking if
                     // the opposite normal had led to a better vector.
                     let cuboid1 = Cuboid::new(radius1);
-                    let sp1 = cuboid1.local_support_point(-penetration_dir1) + vox1.center;
+                    let sp1 = cuboid1.local_support_point(-penetration_dir1) + vox1.center();
                     let sm2 = shape2
                         .as_support_map()
                         .expect("Unsupported collision pair.");
@@ -298,9 +301,9 @@ pub fn contact_manifolds_voxels_shape<ManifoldData, ContactData>(
                 }
 
                 let pt_in_voxel_space = if flipped {
-                    manifold.subshape_pos2().transform_point(pt.local_p2) - vox1.center
+                    manifold.subshape_pos2().transform_point(pt.local_p2) - vox1.center()
                 } else {
-                    manifold.subshape_pos1().transform_point(pt.local_p1) - vox1.center
+                    manifold.subshape_pos1().transform_point(pt.local_p1) - vox1.center()
                 };
                 sub_detector.selected_contacts |=
                     (test_voxel.contains_local_point(pt_in_voxel_space) as u32) << i;
@@ -334,8 +337,11 @@ pub(crate) struct CanonicalVoxelShape {
 }
 
 impl CanonicalVoxelShape {
-    pub fn from_voxel(voxels: &Voxels, vox: &VoxelData) -> Self {
-        let mut key_low = vox.grid_coords;
+    pub fn from_voxel<'a, 'b, V: ?Sized + VoxelQuery>(
+        voxels: &V,
+        vox: &'a impl QueriedVoxel<'b>,
+    ) -> Self {
+        let mut key_low = vox.grid_coords();
         let mut key_high = key_low;
 
         // NOTE: the mins/maxs here are offset by 1 so we can expand past the last voxel if it
@@ -344,7 +350,7 @@ impl CanonicalVoxelShape {
         let mins = voxels.domain()[0] - IVector::splat(1);
         let maxs = voxels.domain()[1];
         let counts = maxs - mins;
-        let mask1 = vox.state.free_faces();
+        let mask1 = vox.voxel_state().free_faces();
 
         let adjust_canon = |axis: AxisMask, i: usize, key: &mut IVector, val: Int| {
             if !mask1.contains(axis) {
@@ -381,17 +387,22 @@ impl CanonicalVoxelShape {
         }
     }
 
-    pub fn cuboid(&self, voxels: &Voxels, vox: &VoxelData, domain2_1: Aabb) -> (Vector, Cuboid) {
+    pub fn cuboid<'a, 'b, V: ?Sized + VoxelQuery>(
+        &self,
+        voxels: &V,
+        vox: &'a impl QueriedVoxel<'b>,
+        domain2_1: Aabb,
+    ) -> (Vector, Cuboid) {
         let radius = voxels.voxel_size() / 2.0;
         let mut canonical_mins = voxels.voxel_center(self.range[0]);
         let mut canonical_maxs = voxels.voxel_center(self.range[1]);
 
         for k in 0..DIM {
-            if self.range[0].ivget(k) != vox.grid_coords.ivget(k) {
+            if self.range[0].ivget(k) != vox.grid_coords().ivget(k) {
                 canonical_mins.vset(k, canonical_mins.vget(k).max(domain2_1.mins.vget(k)));
             }
 
-            if self.range[1].ivget(k) != vox.grid_coords.ivget(k) {
+            if self.range[1].ivget(k) != vox.grid_coords().ivget(k) {
                 canonical_maxs.vset(k, canonical_maxs.vget(k).min(domain2_1.maxs.vget(k)));
             }
         }
diff --git a/src/query/contact_manifolds/contact_manifolds_voxels_voxels.rs b/src/query/contact_manifolds/contact_manifolds_voxels_voxels.rs
index 630ab3a6..96a68c85 100644
--- a/src/query/contact_manifolds/contact_manifolds_voxels_voxels.rs
+++ b/src/query/contact_manifolds/contact_manifolds_voxels_voxels.rs
@@ -6,7 +6,7 @@ use crate::query::{
     ContactManifold, ContactManifoldsWorkspace, PersistentQueryDispatcher, PointQuery,
     TypedWorkspaceData, WorkspaceData,
 };
-use crate::shape::{Cuboid, Shape, SupportMap, VoxelData, VoxelType, Voxels};
+use crate::shape::{Cuboid, QueriedVoxel, Shape, SupportMap, VoxelQuery, VoxelType};
 use crate::utils::hashmap::Entry;
 use crate::utils::PoseOpt;
 use alloc::{boxed::Box, vec::Vec};
@@ -41,18 +41,22 @@ pub fn contact_manifolds_voxels_voxels_shapes<ManifoldData, ContactData>(
     }
 }
 
-/// Computes the contact manifold between a convex shape and a ball.
-pub fn contact_manifolds_voxels_voxels<'a, ManifoldData, ContactData>(
+/// Computes the contact manifold between two voxels shapes.
+///
+/// The voxels shapes can be any voxel storages implementing [`VoxelQuery`].
+pub fn contact_manifolds_voxels_voxels<'a, ManifoldData, ContactData, V1, V2>(
     dispatcher: &dyn PersistentQueryDispatcher<ManifoldData, ContactData>,
     pos12: &Pose,
-    voxels1: &'a Voxels,
-    voxels2: &'a Voxels,
+    voxels1: &'a V1,
+    voxels2: &'a V2,
     prediction: Real,
     manifolds: &mut Vec<ContactManifold<ManifoldData, ContactData>>,
     workspace: &mut Option<ContactManifoldsWorkspace>,
 ) where
     ManifoldData: Default + Clone,
     ContactData: Default + Copy,
+    V1: ?Sized + VoxelQuery,
+    V2: ?Sized + VoxelQuery,
 {
     VoxelsShapeContactManifoldsWorkspace::<4>::ensure_exists(workspace);
     let workspace: &mut VoxelsShapeContactManifoldsWorkspace<4> =
@@ -77,19 +81,25 @@ pub fn contact_manifolds_voxels_voxels<'a, ManifoldData, ContactData>(
         aabb1.aligned_intersections(pos12, &aabb2)
     {
         let domain_margin = (radius1 + radius2) * 10.0;
-        let full_domain2_1 = voxels2.compute_aabb(pos12).add_half_extents(domain_margin);
+        let full_domain2_1 = voxels2
+            .local_aabb()
+            .transform_by(pos12)
+            .add_half_extents(domain_margin);
         let domain2_1 = full_domain2_1
             .intersection(&aabb1.add_half_extents(domain_margin))
             .unwrap_or(full_domain2_1);
-        let full_domain1_2 = voxels1.compute_aabb(&pos21).add_half_extents(domain_margin);
+        let full_domain1_2 = voxels1
+            .local_aabb()
+            .transform_by(&pos21)
+            .add_half_extents(domain_margin);
         let domain1_2 = full_domain1_2
             .intersection(&aabb2.add_half_extents(domain_margin))
             .unwrap_or(full_domain1_2);
 
         let mut detect_hit = |canon1: CanonicalVoxelShape,
                               canon2: CanonicalVoxelShape,
-                              vox1: &VoxelData,
-                              vox2: &VoxelData| {
+                              vox1: &V1::Voxel<'_>,
+                              vox2: &V2::Voxel<'_>| {
             // Compute canonical shapes and dispatch.
             let workspace_key = [
                 canon1.workspace_key[0],
@@ -129,8 +139,8 @@ pub fn contact_manifolds_voxels_voxels<'a, ManifoldData, ContactData>(
                         };
 
                         manifolds.push(ContactManifold::with_data(
-                            vox1.linear_id.flat_id() as u32,
-                            vox2.linear_id.flat_id() as u32,
+                            vox1.linear_id(),
+                            vox2.linear_id(),
                             ManifoldData::default(),
                         ));
 
@@ -204,9 +214,9 @@ pub fn contact_manifolds_voxels_voxels<'a, ManifoldData, ContactData>(
                     // the opposite normal had led to a better vector.
                     let cuboid1 = Cuboid::new(radius1);
                     let cuboid2 = Cuboid::new(radius2);
-                    let sp1 = cuboid1.local_support_point(-penetration_dir1) + vox1.center;
+                    let sp1 = cuboid1.local_support_point(-penetration_dir1) + vox1.center();
                     let sp2 = cuboid2.support_point(
-                        &(pos12 * Pose::from_translation(vox2.center)),
+                        &(pos12 * Pose::from_translation(vox2.center())),
                         penetration_dir1,
                     );
                     let test_dist = (sp2 - sp1).dot(-penetration_dir1);
@@ -220,9 +230,9 @@ pub fn contact_manifolds_voxels_voxels<'a, ManifoldData, ContactData>(
                 }
 
                 let pt_in_voxel_space1 =
-                    manifold.subshape_pos1().transform_point(pt.local_p1) - vox1.center;
+                    manifold.subshape_pos1().transform_point(pt.local_p1) - vox1.center();
                 let pt_in_voxel_space2 =
-                    manifold.subshape_pos2().transform_point(pt.local_p2) - vox2.center;
+                    manifold.subshape_pos2().transform_point(pt.local_p2) - vox2.center();
                 sub_detector.selected_contacts |=
                     ((test_voxel1.contains_local_point(pt_in_voxel_space1) as u32) << i)
                         & ((test_voxel2.contains_local_point(pt_in_voxel_space2) as u32) << i);
@@ -230,7 +240,7 @@ pub fn contact_manifolds_voxels_voxels<'a, ManifoldData, ContactData>(
         };
 
         for vox1 in voxels1.voxels_intersecting_local_aabb(&intersection_aabb1) {
-            let type1 = vox1.state.voxel_type();
+            let type1 = vox1.voxel_type();
             match type1 {
                 #[cfg(feature = "dim2")]
                 VoxelType::Vertex => { /* Ok */ }
@@ -241,10 +251,10 @@ pub fn contact_manifolds_voxels_voxels<'a, ManifoldData, ContactData>(
 
             let canon1 = CanonicalVoxelShape::from_voxel(voxels1, &vox1);
             let centered_aabb1_2 = Cuboid::new(radius1 + Vector::splat(prediction))
-                .compute_aabb(&(pos21 * Pose::from_translation(vox1.center)));
+                .compute_aabb(&(pos21 * Pose::from_translation(vox1.center())));
 
             for vox2 in voxels2.voxels_intersecting_local_aabb(&centered_aabb1_2) {
-                let type2 = vox2.state.voxel_type();
+                let type2 = vox2.voxel_type();
 
                 #[cfg(feature = "dim2")]
                 match (type1, type2) {
@@ -271,7 +281,7 @@ pub fn contact_manifolds_voxels_voxels<'a, ManifoldData, ContactData>(
         }
 
         for vox2 in voxels2.voxels_intersecting_local_aabb(&intersection_aabb2) {
-            let type2 = vox2.state.voxel_type();
+            let type2 = vox2.voxel_type();
             match type2 {
                 #[cfg(feature = "dim2")]
                 VoxelType::Vertex => { /* Ok */ }
@@ -282,10 +292,10 @@ pub fn contact_manifolds_voxels_voxels<'a, ManifoldData, ContactData>(
 
             let canon2 = CanonicalVoxelShape::from_voxel(voxels2, &vox2);
             let centered_aabb2_1 = Cuboid::new(radius2 + Vector::splat(prediction))
-                .compute_aabb(&(pos12 * Pose::from_translation(vox2.center)));
+                .compute_aabb(&(pos12 * Pose::from_translation(vox2.center())));
 
             for vox1 in voxels1.voxels_intersecting_local_aabb(&centered_aabb2_1) {
-                let type1 = vox1.state.voxel_type();
+                let type1 = vox1.voxel_type();
 
                 #[cfg(feature = "dim2")]
                 match (type1, type2) {
diff --git a/src/query/contact_manifolds/mod.rs b/src/query/contact_manifolds/mod.rs
index 7fe44fb2..89885dc4 100644
--- a/src/query/contact_manifolds/mod.rs
+++ b/src/query/contact_manifolds/mod.rs
@@ -163,7 +163,9 @@ pub use self::contact_manifolds_pfm_pfm::{
 pub use self::contact_manifolds_trimesh_shape::{
     contact_manifolds_trimesh_shape, contact_manifolds_trimesh_shape_shapes,
 };
-pub use self::contact_manifolds_voxels_ball::contact_manifolds_voxels_ball_shapes;
+pub use self::contact_manifolds_voxels_ball::{
+    contact_manifolds_voxels_ball, contact_manifolds_voxels_ball_shapes,
+};
 pub use self::contact_manifolds_voxels_composite_shape::{
     contact_manifolds_voxels_composite_shape, contact_manifolds_voxels_composite_shape_shapes,
 };
diff --git a/src/query/intersection_test/intersection_test_voxels_shape.rs b/src/query/intersection_test/intersection_test_voxels_shape.rs
index 5e36affc..9ed63b39 100644
--- a/src/query/intersection_test/intersection_test_voxels_shape.rs
+++ b/src/query/intersection_test/intersection_test_voxels_shape.rs
@@ -1,6 +1,6 @@
 use crate::math::Pose;
 use crate::query::PersistentQueryDispatcher;
-use crate::shape::{Cuboid, Shape, VoxelType, Voxels};
+use crate::shape::{Cuboid, QueriedVoxel, Shape, VoxelQuery, VoxelType};
 
 /// Checks for any intersection between voxels and an arbitrary shape, both represented as a `Shape` trait-object.
 pub fn intersection_test_voxels_shape_shapes(
@@ -19,10 +19,12 @@ pub fn intersection_test_voxels_shape_shapes(
 }
 
 /// Checks for any intersection between voxels and an arbitrary shape.
-pub fn intersection_test_voxels_shape(
+///
+/// The voxels shape can be any voxel storage implementing [`VoxelQuery`].
+pub fn intersection_test_voxels_shape<V: ?Sized + VoxelQuery>(
     dispatcher: &dyn PersistentQueryDispatcher,
     pos12: &Pose,
-    voxels1: &Voxels,
+    voxels1: &V,
     shape2: &dyn Shape,
 ) -> bool {
     let radius1 = voxels1.voxel_size() / 2.0;
@@ -31,13 +33,13 @@ pub fn intersection_test_voxels_shape(
 
     if let Some(intersection_aabb1) = aabb1.intersection(&aabb2_1) {
         for vox1 in voxels1.voxels_intersecting_local_aabb(&intersection_aabb1) {
-            let vox_type1 = vox1.state.voxel_type();
+            let vox_type1 = vox1.voxel_type();
 
             if vox_type1 == VoxelType::Empty {
                 continue;
             }
 
-            let center1 = vox1.center;
+            let center1 = vox1.center();
             let cuboid1 = Cuboid::new(radius1);
             let cuboid_pose12 = Pose::from_translation(-center1) * pos12;
 
@@ -54,11 +56,13 @@ pub fn intersection_test_voxels_shape(
 }
 
 /// Checks for any intersection between voxels and an arbitrary shape.
-pub fn intersection_test_shape_voxels(
+///
+/// The voxels shape can be any voxel storage implementing [`VoxelQuery`].
+pub fn intersection_test_shape_voxels<V: ?Sized + VoxelQuery>(
     dispatcher: &dyn PersistentQueryDispatcher,
     pos12: &Pose,
     shape1: &dyn Shape,
-    voxels2: &Voxels,
+    voxels2: &V,
 ) -> bool {
     intersection_test_voxels_shape(dispatcher, &pos12.inverse(), voxels2, shape1)
 }
diff --git a/src/query/nonlinear_shape_cast/nonlinear_shape_cast_voxels_shape.rs b/src/query/nonlinear_shape_cast/nonlinear_shape_cast_voxels_shape.rs
index 123a8e32..b2961c64 100644
--- a/src/query/nonlinear_shape_cast/nonlinear_shape_cast_voxels_shape.rs
+++ b/src/query/nonlinear_shape_cast/nonlinear_shape_cast_voxels_shape.rs
@@ -1,13 +1,15 @@
 use crate::bounding_volume::BoundingVolume;
 use crate::math::{IVector, IVectorExt, Real, Vector, VectorExt};
 use crate::query::{NonlinearRigidMotion, QueryDispatcher, ShapeCastHit};
-use crate::shape::{Cuboid, Shape, Voxels};
+use crate::shape::{Cuboid, QueriedVoxel, Shape, VoxelQuery, VoxelType};
 
 /// Time Of Impact of a voxels shape with any other shape, under a rigid motion (translation + rotation).
-pub fn cast_shapes_nonlinear_voxels_shape<D>(
+///
+/// The voxels shape can be any voxel storage implementing [`VoxelQuery`].
+pub fn cast_shapes_nonlinear_voxels_shape<D, V>(
     dispatcher: &D,
     motion1: &NonlinearRigidMotion,
-    g1: &Voxels,
+    g1: &V,
     motion2: &NonlinearRigidMotion,
     g2: &dyn Shape,
     start_time: Real,
@@ -16,6 +18,7 @@ pub fn cast_shapes_nonlinear_voxels_shape<D>(
 ) -> Option<ShapeCastHit>
 where
     D: ?Sized + QueryDispatcher,
+    V: ?Sized + VoxelQuery,
 {
     use num_traits::Bounded;
 
@@ -58,9 +61,9 @@ where
 
     let mut check_voxels_in_range = |search_domain: [IVector; 2]| {
         for vox in g1.voxels_in_range(search_domain[0], search_domain[1]) {
-            if !vox.state.is_empty() {
+            if vox.voxel_type() != VoxelType::Empty {
                 // PERF: could we check the canonical shape instead, and deduplicate accordingly?
-                let center = g1.voxel_center(vox.grid_coords);
+                let center = g1.voxel_center(vox.grid_coords());
                 let cuboid = Cuboid::new(g1.voxel_size() / 2.0);
                 let vox_motion1 = motion1.prepend_translation(center);
                 if let Some(new_hit) = dispatcher
@@ -170,19 +173,22 @@ where
     hit
 }
 
-/// Time Of Impact of any shape with a composite shape, under a rigid motion (translation + rotation).
-pub fn cast_shapes_nonlinear_shape_voxels<D>(
+/// Time Of Impact of any shape with a voxels shape, under a rigid motion (translation + rotation).
+///
+/// The voxels shape can be any voxel storage implementing [`VoxelQuery`].
+pub fn cast_shapes_nonlinear_shape_voxels<D, V>(
     dispatcher: &D,
     motion1: &NonlinearRigidMotion,
     g1: &dyn Shape,
     motion2: &NonlinearRigidMotion,
-    g2: &Voxels,
+    g2: &V,
     start_time: Real,
     end_time: Real,
     stop_at_penetration: bool,
 ) -> Option<ShapeCastHit>
 where
     D: ?Sized + QueryDispatcher,
+    V: ?Sized + VoxelQuery,
 {
     cast_shapes_nonlinear_voxels_shape(
         dispatcher,
diff --git a/src/query/point/point_voxels.rs b/src/query/point/point_voxels.rs
index 87a5b147..6b42e025 100644
--- a/src/query/point/point_voxels.rs
+++ b/src/query/point/point_voxels.rs
@@ -1,6 +1,6 @@
 use crate::math::{Real, Vector};
 use crate::query::{PointProjection, PointQuery};
-use crate::shape::{Cuboid, FeatureId, Voxels, VoxelsChunkRef};
+use crate::shape::{Cuboid, FeatureId, QueriedVoxel, Voxels, VoxelsChunkRef};
 
 impl PointQuery for Voxels {
     #[inline]
@@ -49,17 +49,17 @@ impl<'a> VoxelsChunkRef<'a> {
         let mut result_vox_id = 0;
 
         for vox in self.voxels() {
-            let mut candidate = base_cuboid.project_local_point(pt - vox.center, solid);
-            candidate.point += vox.center;
+            let mut candidate = base_cuboid.project_local_point(pt - vox.center(), solid);
+            candidate.point += vox.center();
 
             let candidate_dist = (candidate.point - pt).length();
             if candidate_dist < smallest_dist {
                 result = candidate;
-                result_vox_id = vox.linear_id.flat_id();
+                result_vox_id = vox.linear_id;
                 smallest_dist = candidate_dist;
             }
         }
 
-        (smallest_dist < Real::MAX).then_some((result, result_vox_id as u32))
+        (smallest_dist < Real::MAX).then_some((result, result_vox_id))
     }
 }
diff --git a/src/query/shape_cast/shape_cast_voxels_shape.rs b/src/query/shape_cast/shape_cast_voxels_shape.rs
index baa8ada1..55ce67e6 100644
--- a/src/query/shape_cast/shape_cast_voxels_shape.rs
+++ b/src/query/shape_cast/shape_cast_voxels_shape.rs
@@ -1,18 +1,21 @@
 use crate::math::{IVector, IVectorExt, Pose, Real, Vector, VectorExt};
 use crate::query::{QueryDispatcher, ShapeCastHit, ShapeCastOptions};
-use crate::shape::{Cuboid, Shape, Voxels};
+use crate::shape::{Cuboid, QueriedVoxel, Shape, VoxelQuery, VoxelType};
 
 /// Time Of Impact of a voxels shape with any other shape, under a translational movement.
-pub fn cast_shapes_voxels_shape<D>(
+///
+/// The voxels shape can be any voxel storage implementing [`VoxelQuery`].
+pub fn cast_shapes_voxels_shape<D, V>(
     dispatcher: &D,
     pos12: &Pose,
     vel12: Vector,
-    g1: &Voxels,
+    g1: &V,
     g2: &dyn Shape,
     options: ShapeCastOptions,
 ) -> Option<ShapeCastHit>
 where
     D: ?Sized + QueryDispatcher,
+    V: ?Sized + VoxelQuery,
 {
     use num_traits::Bounded;
 
@@ -22,9 +25,9 @@ where
 
     let mut check_voxels_in_range = |search_domain: [IVector; 2]| {
         for vox in g1.voxels_in_range(search_domain[0], search_domain[1]) {
-            if !vox.state.is_empty() {
+            if vox.voxel_type() != VoxelType::Empty {
                 // PERF: could we check the canonical shape instead, and deduplicate accordingly?
-                let center = g1.voxel_center(vox.grid_coords);
+                let center = g1.voxel_center(vox.grid_coords());
                 let cuboid = Cuboid::new(g1.voxel_size() / 2.0);
                 let vox_pos12 = Pose::from_translation(center).inverse() * pos12;
                 if let Some(mut new_hit) = dispatcher
@@ -133,17 +136,20 @@ where
     hit
 }
 
-/// Time Of Impact of any shape with a composite shape, under a rigid motion (translation + rotation).
-pub fn cast_shapes_shape_voxels<D>(
+/// Time Of Impact of any shape with a voxels shape, under a translational movement.
+///
+/// The voxels shape can be any voxel storage implementing [`VoxelQuery`].
+pub fn cast_shapes_shape_voxels<D, V>(
     dispatcher: &D,
     pos12: &Pose,
     vel12: Vector,
     g1: &dyn Shape,
-    g2: &Voxels,
+    g2: &V,
     options: ShapeCastOptions,
 ) -> Option<ShapeCastHit>
 where
     D: ?Sized + QueryDispatcher,
+    V: ?Sized + VoxelQuery,
 {
     cast_shapes_voxels_shape(
         dispatcher,
diff --git a/src/shape/mod.rs b/src/shape/mod.rs
index e6b95f01..62c169b7 100644
--- a/src/shape/mod.rs
+++ b/src/shape/mod.rs
@@ -22,7 +22,10 @@ pub use self::{
     compound::Compound,
     polyline::Polyline,
     shared_shape::SharedShape,
-    voxels::{AxisMask, OctantPattern, VoxelData, VoxelState, VoxelType, Voxels, VoxelsChunkRef},
+    voxels::{
+        AxisMask, OctantPattern, QueriedVoxel, VoxelData, VoxelQuery, VoxelState, VoxelType,
+        Voxels, VoxelsChunkRef,
+    },
 };
 
 // `PolylineFlags` is a 2D-only feature.
diff --git a/src/shape/voxels/mod.rs b/src/shape/voxels/mod.rs
index ea8b7c0a..24acf7f1 100644
--- a/src/shape/voxels/mod.rs
+++ b/src/shape/voxels/mod.rs
@@ -1,8 +1,10 @@
+pub use voxel_query::*;
 pub use voxels::*;
 pub use voxels_chunk::*;
 
 use voxels_consts::*;
 
+mod voxel_query;
 mod voxels;
 mod voxels_chunk;
 mod voxels_consts;
diff --git a/src/shape/voxels/voxel_query.rs b/src/shape/voxels/voxel_query.rs
new file mode 100644
index 00000000..08f10b95
--- /dev/null
+++ b/src/shape/voxels/voxel_query.rs
@@ -0,0 +1,227 @@
+use crate::math::{ivect_to_vect, vect_to_ivect, IVector, Vector};
+
+use crate::bounding_volume::Aabb;
+use crate::shape::{VoxelData, VoxelState, VoxelType, Voxels};
+
+/// Abstraction over the storage of a shape made of axis-aligned, uniformly sized voxels.
+///
+/// Parry's voxel collision-detection algorithms (contact manifolds, intersection tests,
+/// linear and nonlinear shape-casting) are
+/// written against this trait rather than against the concrete [`Voxels`] shape. Implementing
+/// it for a custom sparse data-structure (chunked grid, octree, VDB-like tree, etc.) lets these
+/// algorithms run directly on that structure without copying it into a [`Voxels`] shape,
+/// typically by calling the generic query functions from a custom
+/// [`QueryDispatcher`](crate::query::QueryDispatcher).
+///
+/// # Grid conventions
+///
+/// Voxels are identified by their integer grid coordinates `key`. The voxel with coordinates
+/// `key` covers the world-space (well, shape-local-space) range
+/// `[key * voxel_size, (key + 1) * voxel_size]`, so its center is at
+/// `(key + 0.5) * voxel_size`. Grid ranges are always given as semi-open intervals
+/// `[mins, maxs)`: `mins` is included, `maxs` is excluded.
+///
+/// # Voxel views and neighborhood states
+///
+/// Lookups and iterators don't yield a fixed data struct: they yield storage-defined voxel
+/// *views* ([`Self::Voxel`], bounded by [`QueriedVoxel`]). A view exposes cheap per-voxel
+/// data — grid coordinates, center, and the coarse [`QueriedVoxel::voxel_type`], which a
+/// sparse storage can pack in two bits per stored voxel (with empty voxels simply absent).
+///
+/// Contact-manifold computation additionally needs to know *which* of a voxel's immediate
+/// axis-aligned neighbors are filled — a [`VoxelState`] — to avoid hitting the "internal
+/// edges" between adjacent voxels. It obtains this from [`QueriedVoxel::voxel_state`], and
+/// only for the few voxels that are actual contact candidates, never during bulk iteration.
+/// Since views can borrow from their storage, they can compute the state on demand from
+/// local context (e.g. leaf-local reads in a sparse tree);
+/// [`VoxelState::with_filled_neighbors`] builds the state from occupancy alone, while storages
+/// like [`Voxels`] that persist the state (one byte per voxel) just hand out the stored value.
+///
+/// # Note for implementors
+///
+/// This trait is not dyn-compatible (`voxels_in_range` returns `impl Iterator`). The generic
+/// query functions are monomorphized for each storage type. To plug a custom storage into a
+/// physics pipeline, wrap it in a type implementing [`Shape`](crate::shape::Shape) (typically
+/// with [`ShapeType::Custom`](crate::shape::ShapeType::Custom)) and dispatch to the generic
+/// voxel query functions from a custom `QueryDispatcher`.
+pub trait VoxelQuery {
+    /// The view type this storage hands out for a single voxel.
+    ///
+    /// Views can borrow from the storage (e.g. hold a cursor into a sparse tree), letting
+    /// [`QueriedVoxel::voxel_state`] read neighborhood information from local context
+    /// instead of independent whole-storage lookups.
+    type Voxel<'a>: QueriedVoxel<'a>
+    where
+        Self: 'a;
+
+    /// The size of each voxel along each local coordinate axis.
+    fn voxel_size(&self) -> Vector;
+
+    /// The semi-open range `[mins, maxs)` of grid coordinates covered by this shape.
+    ///
+    /// This must be a conservative bound: every non-empty voxel must lie within the returned
+    /// range, but the range may also cover empty voxels.
+    fn domain(&self) -> [IVector; 2];
+
+    /// Iterates through the voxels within the given semi-open grid coordinate range.
+    ///
+    /// Implementations must yield every non-empty voxel with grid coordinates in
+    /// `[mins, maxs)` exactly once. They may additionally yield empty voxels within that
+    /// range (callers filter on [`QueriedVoxel::voxel_type`]), but must never yield a
+    /// voxel outside of the range.
+    fn voxels_in_range(
+        &self,
+        mins: IVector,
+        maxs: IVector,
+    ) -> impl Iterator<Item = Self::Voxel<'_>>;
+
+    /// Iterates through every voxel of this shape.
+    ///
+    /// This is equivalent to [`Self::voxels_in_range`] applied to the whole [`Self::domain`].
+    fn voxels(&self) -> impl Iterator<Item = Self::Voxel<'_>> {
+        let [mins, maxs] = self.domain();
+        self.voxels_in_range(mins, maxs)
+    }
+
+    /// Iterates through every voxel intersecting the given local-space AABB.
+    fn voxels_intersecting_local_aabb(&self, aabb: &Aabb) -> impl Iterator<Item = Self::Voxel<'_>> {
+        let [mins, maxs] = self.voxel_range_intersecting_local_aabb(aabb);
+        self.voxels_in_range(mins, maxs)
+    }
+
+    /// The grid coordinates of the voxel containing the given local-space point.
+    ///
+    /// The returned coordinates are valid regardless of whether the corresponding voxel
+    /// is filled, empty, or outside of [`Self::domain`].
+    fn voxel_at_point(&self, point: Vector) -> IVector {
+        vect_to_ivect((point / self.voxel_size()).floor())
+    }
+
+    /// The local-space center of the voxel with the given grid coordinates.
+    fn voxel_center(&self, key: IVector) -> Vector {
+        (ivect_to_vect(key) + Vector::splat(0.5)) * self.voxel_size()
+    }
+
+    /// The local-space AABB of the voxel with the given grid coordinates.
+    fn voxel_aabb(&self, key: IVector) -> Aabb {
+        let center = self.voxel_center(key);
+        Aabb::from_half_extents(center, self.voxel_size() / 2.0)
+    }
+
+    /// The semi-open range of grid coordinates of the voxels intersecting the given AABB.
+    ///
+    /// The returned range covers both empty and non-empty voxels, and is not limited to the
+    /// bounds defined by [`Self::domain`].
+    fn voxel_range_intersecting_local_aabb(&self, aabb: &Aabb) -> [IVector; 2] {
+        let mins = vect_to_ivect((aabb.mins / self.voxel_size()).floor());
+        let maxs = vect_to_ivect((aabb.maxs / self.voxel_size()).ceil());
+        [mins, maxs]
+    }
+
+    /// The local-space AABB of the given semi-open range of voxel grid coordinates.
+    fn voxel_range_aabb(&self, mins: IVector, maxs: IVector) -> Aabb {
+        Aabb {
+            mins: ivect_to_vect(mins) * self.voxel_size(),
+            maxs: ivect_to_vect(maxs) * self.voxel_size(),
+        }
+    }
+
+    /// Aligns the given AABB with the voxelized grid.
+    ///
+    /// The returned AABB has corners lying at the grid intersections (i.e. matches voxel
+    /// corners) and fully contains the input `aabb`.
+    fn align_aabb_to_grid(&self, aabb: &Aabb) -> Aabb {
+        let mins = (aabb.mins / self.voxel_size()).floor() * self.voxel_size();
+        let maxs = (aabb.maxs / self.voxel_size()).ceil() * self.voxel_size();
+        Aabb { mins, maxs }
+    }
+
+    /// The local-space AABB of this voxels shape.
+    fn local_aabb(&self) -> Aabb {
+        let [mins, maxs] = self.domain();
+        self.voxel_range_aabb(mins, maxs)
+    }
+}
+
+/// A single voxel handed out by a [`VoxelQuery`] storage.
+///
+/// This is the item type of the [`VoxelQuery`] lookups and iterators. Storages define their
+/// own implementor (see [`VoxelQuery::Voxel`]), which may borrow from the storage so that
+/// [`Self::voxel_state`] can be computed lazily from local context.
+///
+/// The coarse [`Self::voxel_type`] must be cheap: it is read during bulk iteration. The full
+/// [`Self::voxel_state`] is only requested for contact-candidate voxels; implementations
+/// must keep the two consistent (`self.voxel_state().voxel_type() == self.voxel_type()`).
+pub trait QueriedVoxel<'a> {
+    /// The type of this voxel: empty, or how it is exposed on the shape's surface.
+    fn voxel_type(&self) -> VoxelType;
+
+    /// The neighborhood state of this voxel, indicating which of its immediate
+    /// axis-aligned neighbors are filled.
+    fn voxel_state(&self) -> VoxelState;
+
+    /// A stable, storage-defined identifier of this voxel.
+    ///
+    /// For the [`Voxels`] shape this is the flattened form of [`Voxels::linear_index`].
+    /// This identifier can be invalidated after the voxels shape is modified (e.g. by a call
+    /// to [`Voxels::set_voxel`], or [`Voxels::crop`]).
+    /// For stable references to voxels, always use [`Self::grid_coords`]. Only meaningful
+    /// for non-empty voxels.
+    fn linear_id(&self) -> u32;
+
+    /// The voxel's integer grid coordinates.
+    fn grid_coords(&self) -> IVector;
+
+    /// The voxel's center position in the local-space of the voxels shape it is part of.
+    fn center(&self) -> Vector;
+}
+
+impl QueriedVoxel<'_> for VoxelData {
+    fn voxel_type(&self) -> VoxelType {
+        self.state.voxel_type()
+    }
+
+    fn voxel_state(&self) -> VoxelState {
+        self.state
+    }
+
+    fn linear_id(&self) -> u32 {
+        self.linear_id
+    }
+
+    fn grid_coords(&self) -> IVector {
+        self.grid_coords
+    }
+
+    fn center(&self) -> Vector {
+        self.center
+    }
+}
+impl VoxelQuery for Voxels {
+    type Voxel<'a> = VoxelData;
+
+    #[inline]
+    fn voxel_size(&self) -> Vector {
+        self.voxel_size()
+    }
+
+    #[inline]
+    fn domain(&self) -> [IVector; 2] {
+        self.domain()
+    }
+
+    #[inline]
+    fn voxels_in_range(&self, mins: IVector, maxs: IVector) -> impl Iterator<Item = VoxelData> {
+        self.voxels_in_range(mins, maxs)
+    }
+
+    #[inline]
+    fn voxels(&self) -> impl Iterator<Item = VoxelData> {
+        self.voxels()
+    }
+
+    #[inline]
+    fn local_aabb(&self) -> Aabb {
+        self.local_aabb()
+    }
+}
diff --git a/src/shape/voxels/voxels.rs b/src/shape/voxels/voxels.rs
index 892b9f50..1a206733 100644
--- a/src/shape/voxels/voxels.rs
+++ b/src/shape/voxels/voxels.rs
@@ -55,9 +55,10 @@ use alloc::{vec, vec::Vec};
 /// println!("Voxel type: {:?}", voxel_type);
 /// # }
 /// ```
-#[derive(Copy, Clone, Debug, PartialEq, Eq)]
+#[derive(Copy, Clone, Debug, PartialEq, Eq, Default)]
 pub enum VoxelType {
     /// The voxel is empty.
+    #[default]
     Empty,
     /// The voxel is a vertex if all three coordinate axis directions have at
     /// least one empty neighbor.
@@ -218,6 +219,12 @@ impl OctantPattern {
 #[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]
 pub struct VoxelState(pub(super) u8);
 
+impl Default for VoxelState {
+    fn default() -> Self {
+        Self::EMPTY
+    }
+}
+
 impl VoxelState {
     /// The value of empty voxels.
     pub const EMPTY: VoxelState = VoxelState(EMPTY_FACE_MASK);
@@ -228,6 +235,14 @@ impl VoxelState {
         Self(state)
     }
 
+    /// The state of a **non-empty** voxel given the set of its non-empty axis-aligned neighbors.
+    /// ```
+    pub const fn with_filled_neighbors(filled_neighbors: AxisMask) -> Self {
+        // The `AxisMask` bits match the internal neighborhood bit layout:
+        // bit `2 * axis` is the positive direction, bit `2 * axis + 1` the negative one.
+        Self(filled_neighbors.bits())
+    }
+
     /// Is this voxel empty?
     pub const fn is_empty(self) -> bool {
         self.0 == EMPTY_FACE_MASK
@@ -261,7 +276,9 @@ impl VoxelState {
 /// Information associated to a voxel.
 ///
 /// This structure provides complete information about a single voxel including its position
-/// in both grid coordinates and world space, as well as its state (empty/filled and neighborhood).
+/// in both grid coordinates and world space, as well as its state (empty/filled and
+/// neighborhood). It is the view type ([`VoxelQuery::Voxel`](crate::shape::VoxelQuery::Voxel))
+/// handed out by the [`Voxels`] shape.
 ///
 /// # Note
 ///
@@ -283,21 +300,21 @@ impl VoxelState {
 ///
 /// // Iterate through all voxels
 /// for voxel in voxels.voxels() {
-///     if !voxel.state.is_empty() {
-///         println!("Voxel at grid position {:?}", voxel.grid_coords);
-///         println!("  World center: {:?}", voxel.center);
-///         println!("  Type: {:?}", voxel.state.voxel_type());
-///     }
+///     println!("Voxel at grid position {:?}", voxel.grid_coords);
+///     println!("  World center: {:?}", voxel.center);
+///     println!("  Type: {:?}", voxel.state.voxel_type());
 /// }
 /// # }
 /// ```
 #[derive(Copy, Clone, Debug, PartialEq)]
 pub struct VoxelData {
-    /// The temporary index in the internal voxels' storage.
+    /// A stable, storage-defined identifier of this voxel.
     ///
-    /// This index can be invalidated after a call to [`Voxels::set_voxel`], or
-    /// [`Voxels::crop`].
-    pub linear_id: VoxelIndex,
+    /// For the [`Voxels`] shape this is the flattened form of [`Voxels::linear_index`].
+    /// This identifier can be invalidated after the voxels shape is modified (e.g. by a call
+    /// to [`Voxels::set_voxel`], or [`Voxels::crop`]).
+    /// For stable references to voxels, always use `grid_coords`.
+    pub linear_id: u32,
     /// The voxel's integer grid coordinates.
     pub grid_coords: IVector,
     /// The voxel's center position in the local-space of the [`Voxels`] shape it is part of.
@@ -430,10 +447,8 @@ pub struct VoxelData {
 ///
 /// // Iterate through all non-empty voxels
 /// for voxel in voxels.voxels() {
-///     if !voxel.state.is_empty() {
-///         println!("Voxel at grid {:?}, world center {:?}",
-///                  voxel.grid_coords, voxel.center);
-///     }
+///     println!("Voxel at grid {:?}, world center {:?}",
+///              voxel.grid_coords, voxel.center);
 /// }
 /// # }
 /// ```
@@ -477,9 +492,7 @@ pub struct VoxelData {
 ///
 /// // Find voxels intersecting an AABB
 /// let query_aabb = Aabb::new(Vector::new(-0.5, -0.5, -0.5), Vector::new(1.5, 1.5, 1.5));
-/// let count = voxels.voxels_intersecting_local_aabb(&query_aabb)
-///     .filter(|v| !v.state.is_empty())
-///     .count();
+/// let count = voxels.voxels_intersecting_local_aabb(&query_aabb).count();
 /// println!("Found {} voxels in AABB", count);
 ///
 /// // Get the overall domain bounds
@@ -646,9 +659,7 @@ impl Voxels {
     /// let voxels = Voxels::from_points(Vector::new(1.0, 1.0, 1.0), &points);
     ///
     /// // Only 3 unique voxels created (first two points merged)
-    /// let filled_count = voxels.voxels()
-    ///     .filter(|v| !v.state.is_empty())
-    ///     .count();
+    /// let filled_count = voxels.voxels().count();
     /// assert_eq!(filled_count, 3);
     /// # }
     /// ```
@@ -714,9 +725,7 @@ impl Voxels {
     ///
     /// // Iterate through filled voxels (more efficient than iterating domain)
     /// for voxel in voxels.voxels() {
-    ///     if !voxel.state.is_empty() {
-    ///         println!("Filled voxel at {:?}", voxel.grid_coords);
-    ///     }
+    ///     println!("Filled voxel at {:?}", voxel.grid_coords);
     /// }
     /// # }
     /// ```
@@ -897,8 +906,6 @@ impl Voxels {
     }
 
     /// Iterates through every voxel intersecting the given aabb.
-    ///
-    /// Returns the voxel’s linearized id, center, and state.
     pub fn voxels_intersecting_local_aabb(
         &self,
         aabb: &Aabb,
@@ -919,7 +926,8 @@ impl Voxels {
         )
     }
 
-    /// Iterate through the data of all the voxels within the given (semi-open) voxel grid indices.
+    /// Iterate through the data of all the non-empty voxels within the given (semi-open)
+    /// voxel grid indices.
     ///
     /// Note that this only yields non-empty voxels within the range. This does not
     /// include any voxel that falls outside [`Self::domain`].
diff --git a/src/shape/voxels/voxels_chunk.rs b/src/shape/voxels/voxels_chunk.rs
index 47603f09..e00bb7ac 100644
--- a/src/shape/voxels/voxels_chunk.rs
+++ b/src/shape/voxels/voxels_chunk.rs
@@ -137,9 +137,7 @@ impl VoxelsChunk {
 ///
 /// // Query voxels within this chunk
 /// for voxel in chunk_ref.voxels() {
-///     if !voxel.state.is_empty() {
-///         println!("Voxel at {:?}", voxel.grid_coords);
-///     }
+///     println!("Voxel at {:?}", voxel.grid_coords);
 /// }
 ///
 /// // Get chunk's AABB
@@ -259,7 +257,8 @@ impl<'a> VoxelsChunkRef<'a> {
                     linear_id: VoxelIndex {
                         chunk_id: self.my_id,
                         id_in_chunk,
-                    },
+                    }
+                    .flat_id() as u32,
                     grid_coords,
                     center,
                     state,
@@ -303,7 +302,8 @@ impl<'a> VoxelsChunkRef<'a> {
                         linear_id: VoxelIndex {
                             chunk_id: self.my_id,
                             id_in_chunk,
-                        },
+                        }
+                        .flat_id() as u32,
                         grid_coords,
                         center,
                         state,
diff --git a/src/shape/voxels/voxels_edition.rs b/src/shape/voxels/voxels_edition.rs
index 3c962b02..a0d8138e 100644
--- a/src/shape/voxels/voxels_edition.rs
+++ b/src/shape/voxels/voxels_edition.rs
@@ -89,9 +89,7 @@ impl Voxels {
     /// }
     ///
     /// // Count filled voxels
-    /// let filled = voxels.voxels()
-    ///     .filter(|v| !v.state.is_empty())
-    ///     .count();
+    /// let filled = voxels.voxels().count();
     /// assert_eq!(filled, 9);
     /// # }
     /// ```
@@ -199,9 +197,7 @@ impl Voxels {
     /// voxels.crop(IVector::new(1, 0, 0), IVector::new(2, 0, 0));
     ///
     /// // Only two voxels remain
-    /// let count = voxels.voxels()
-    ///     .filter(|v| !v.state.is_empty())
-    ///     .count();
+    /// let count = voxels.voxels().count();
     /// assert_eq!(count, 2);
     /// # }
     /// ```
diff --git a/src/transformation/to_outline/voxels_to_outline.rs b/src/transformation/to_outline/voxels_to_outline.rs
index dd9fb0f6..b0635b95 100644
--- a/src/transformation/to_outline/voxels_to_outline.rs
+++ b/src/transformation/to_outline/voxels_to_outline.rs
@@ -1,6 +1,6 @@
 use crate::bounding_volume::Aabb;
 use crate::math::Vector;
-use crate::shape::{VoxelType, Voxels};
+use crate::shape::{QueriedVoxel, VoxelType, Voxels};
 use alloc::{vec, vec::Vec};
 
 impl Voxels {
@@ -29,9 +29,10 @@ impl Voxels {
         let vtx = aabb.vertices();
 
         for vox in self.voxels() {
-            match vox.state.voxel_type() {
+            let state = vox.voxel_state();
+            match vox.voxel_type() {
                 VoxelType::Vertex => {
-                    let mask = vox.state.feature_mask();
+                    let mask = state.feature_mask();
 
                     for edge in Aabb::EDGES_VERTEX_IDS {
                         if mask & (1 << edge.0) != 0 || mask & (1 << edge.1) != 0 {
@@ -41,7 +42,7 @@ impl Voxels {
                 }
                 VoxelType::Edge => {
                     let vtx = aabb.vertices();
-                    let mask = vox.state.feature_mask();
+                    let mask = state.feature_mask();
 
                     for (i, edge) in Aabb::EDGES_VERTEX_IDS.iter().enumerate() {
                         if mask & (1 << i) != 0 {
diff --git a/src/transformation/to_polyline/voxels_to_polyline.rs b/src/transformation/to_polyline/voxels_to_polyline.rs
index 727f3f40..e0cc7a10 100644
--- a/src/transformation/to_polyline/voxels_to_polyline.rs
+++ b/src/transformation/to_polyline/voxels_to_polyline.rs
@@ -1,6 +1,6 @@
 use crate::bounding_volume::Aabb;
 use crate::math::{Vector, Vector2};
-use crate::shape::{VoxelType, Voxels};
+use crate::shape::{QueriedVoxel, VoxelType, Voxels};
 use alloc::{vec, vec::Vec};
 
 impl Voxels {
@@ -25,9 +25,10 @@ impl Voxels {
         let vtx = aabb.vertices();
 
         for vox in self.voxels() {
-            match vox.state.voxel_type() {
+            let state = vox.voxel_state();
+            match state.voxel_type() {
                 VoxelType::Vertex => {
-                    let mask = vox.state.feature_mask();
+                    let mask = state.feature_mask();
 
                     for edge in Aabb::FACES_VERTEX_IDS {
                         if mask & (1 << edge.0) != 0 || mask & (1 << edge.1) != 0 {
@@ -37,7 +38,7 @@ impl Voxels {
                 }
                 VoxelType::Face => {
                     let vtx = aabb.vertices();
-                    let mask = vox.state.feature_mask();
+                    let mask = state.feature_mask();
 
                     for (i, edge) in Aabb::FACES_VERTEX_IDS.iter().enumerate() {
                         if mask & (1 << i) != 0 {
diff --git a/src/transformation/to_trimesh/voxels_to_trimesh.rs b/src/transformation/to_trimesh/voxels_to_trimesh.rs
index 676cc68b..10d97f54 100644
--- a/src/transformation/to_trimesh/voxels_to_trimesh.rs
+++ b/src/transformation/to_trimesh/voxels_to_trimesh.rs
@@ -1,6 +1,6 @@
-use crate::bounding_volume::Aabb;
 use crate::math::Vector;
 use crate::shape::Voxels;
+use crate::{bounding_volume::Aabb, shape::QueriedVoxel};
 use alloc::{vec, vec::Vec};
 
 impl Voxels {
@@ -15,7 +15,8 @@ impl Voxels {
         let mut vtx = vec![];
         let mut idx = vec![];
         for vox in self.voxels() {
-            let mask = vox.state.free_faces();
+            let state = vox.voxel_state();
+            let mask = state.free_faces();
             for i in 0..6 {
                 if mask.bits() & (1 << i) != 0 {
                     let fvid = Aabb::FACES_VERTEX_IDS[i];