diff --git a/examples/bevy_wgpu_render.rs b/examples/bevy_wgpu_render.rs
index dcfb64f..c1df140 100644
--- a/examples/bevy_wgpu_render.rs
+++ b/examples/bevy_wgpu_render.rs
@@ -11,7 +11,7 @@ use shocovox_rs::octree::{
 const DISPLAY_RESOLUTION: [u32; 2] = [1024, 768];
 
 #[cfg(feature = "bevy_wgpu")]
-const ARRAY_DIMENSION: u32 = 64;
+const ARRAY_DIMENSION: u32 = 128;
 
 #[cfg(feature = "bevy_wgpu")]
 fn main() {
diff --git a/src/octree/raytracing/raytracing_on_cpu.rs b/src/octree/raytracing/raytracing_on_cpu.rs
index 651a639..06e2844 100644
--- a/src/octree/raytracing/raytracing_on_cpu.rs
+++ b/src/octree/raytracing/raytracing_on_cpu.rs
@@ -1,51 +1,38 @@
-use crate::octree::{
-    types::{NodeChildrenArray, NodeContent},
-    Cube, Octree, V3c, VoxelData,
+use crate::{
+    octree::{
+        types::{NodeChildrenArray, NodeContent},
+        Cube, Octree, V3c, VoxelData,
+    },
+    spatial::math::step_octant,
 };
 
-use crate::spatial::math::cube_impact_normal;
 use crate::spatial::{
     math::{
-        flat_projection, hash_direction, hash_region, octant_bitmask, offset_region,
+        cube_impact_normal, flat_projection, hash_direction, hash_region, octant_bitmask,
         position_in_bitmap_64bits,
     },
     raytracing::{
-        lut::RAY_TO_LEAF_OCCUPANCY_BITMASK_LUT, lut::RAY_TO_NODE_OCCUPANCY_BITMASK_LUT, Ray,
+        lut::{OOB_OCTANT, RAY_TO_LEAF_OCCUPANCY_BITMASK_LUT, RAY_TO_NODE_OCCUPANCY_BITMASK_LUT},
+        Ray,
     },
     FLOAT_ERROR_TOLERANCE,
 };
 
+#[derive(Debug)]
 struct NodeStackItem {
     pub(crate) bounds: Cube,
-    pub(crate) node: u32,
-    pub(crate) target_octant: u8,
-    pub(crate) child_center: V3c<f32>,
+    node: u32,
+    target_octant: u8,
 }
 
 impl NodeStackItem {
     pub(crate) fn new(bounds: Cube, node: u32, target_octant: u8) -> Self {
-        let child_center = Into::<V3c<f32>>::into(bounds.min_position)
-            + V3c::unit(bounds.size as f32 / 4.)
-            + Into::<V3c<f32>>::into(offset_region(target_octant)) * (bounds.size as f32 / 2.);
         Self {
             bounds,
             node,
             target_octant,
-            child_center,
         }
     }
-
-    pub(crate) fn add_point(&mut self, p: V3c<f32>) {
-        self.child_center = self.child_center + p;
-        self.target_octant = hash_region(
-            &(self.child_center - self.bounds.min_position.into()),
-            self.bounds.size as f32,
-        );
-    }
-
-    pub(crate) fn contains_target_center(&self) -> bool {
-        self.bounds.contains_point(&self.child_center)
-    }
 }
 
 impl<T: Default + PartialEq + Clone + std::fmt::Debug + VoxelData, const DIM: usize>
@@ -128,7 +115,7 @@ impl<T: Default + PartialEq + Clone + std::fmt::Debug + VoxelData, const DIM: us
         direction_lut_index: usize,
     ) -> Option<V3c<usize>> {
         let mut current_index = {
-            let pos = ray.point_at(*ray_current_distance) - V3c::<f32>::from(bounds.min_position);
+            let pos = ray.point_at(*ray_current_distance) - bounds.min_position;
             V3c::new(
                 (pos.x as i32).clamp(0, (DIM - 1) as i32),
                 (pos.y as i32).clamp(0, (DIM - 1) as i32),
@@ -251,6 +238,7 @@ impl<T: Default + PartialEq + Clone + std::fmt::Debug + VoxelData, const DIM: us
         let direction_lut_index = hash_direction(&ray.direction) as usize;
 
         let root_bounds = Cube::root_bounds(self.octree_size as f32);
+        let mut node_current_size = root_bounds.size;
         let mut ray_current_distance = 0.0; // No need to initialize, but it will shut the compiler
         let mut node_stack = Vec::new();
         if let Some(root_hit) = root_bounds.intersect_ray(&ray) {
@@ -315,7 +303,7 @@ impl<T: Default + PartialEq + Clone + std::fmt::Debug + VoxelData, const DIM: us
             }
 
             if leaf_miss
-                || !node_stack_top.contains_target_center() // If current target is OOB
+                || node_stack_top.target_octant == OOB_OCTANT
                 // In case the current Node is empty
                 || 0 == current_node_occupied_bits
                 // In case there is no overlap between the node occupancy and the potential slots the ray would hit
@@ -330,8 +318,9 @@ impl<T: Default + PartialEq + Clone + std::fmt::Debug + VoxelData, const DIM: us
                         &popped_target.bounds,
                         &ray_scale_factors,
                     );
-                    parent.add_point(step_vec);
+                    parent.target_octant = step_octant(parent.target_octant, step_vec);
                 }
+                node_current_size *= 2.;
                 continue; // Re-calculate current_bounds
             }
 
@@ -351,6 +340,7 @@ impl<T: Default + PartialEq + Clone + std::fmt::Debug + VoxelData, const DIM: us
                     target_child_key,
                     child_target_octant,
                 ));
+                node_current_size /= 2.;
             } else {
                 // ADVANCE
                 // target child is invalid, or it does not intersect with the ray
@@ -366,12 +356,14 @@ impl<T: Default + PartialEq + Clone + std::fmt::Debug + VoxelData, const DIM: us
                     if let Some(hit) = target_hit {
                         ray_current_distance = hit.exit_distance;
                     }
-                    node_stack.last_mut().unwrap().add_point(step_vec);
-                    target_octant = node_stack.last().unwrap().target_octant;
-                    target_bounds = current_bounds.child_bounds_for(target_octant);
-                    target_child_key = self.node_children[current_node_key][target_octant as u32];
+                    target_octant = step_octant(target_octant, step_vec);
+                    if OOB_OCTANT != target_octant {
+                        target_bounds = current_bounds.child_bounds_for(target_octant);
+                        target_child_key =
+                            self.node_children[current_node_key][target_octant as u32];
+                    }
 
-                    if !node_stack.last().unwrap().contains_target_center()
+                    if target_octant == OOB_OCTANT
                         || (self.nodes.key_is_valid(target_child_key as usize)
                             // current node is occupied at target octant
                             && 0 != current_node_occupied_bits & octant_bitmask(target_octant)
@@ -381,7 +373,7 @@ impl<T: Default + PartialEq + Clone + std::fmt::Debug + VoxelData, const DIM: us
                                 NodeContent::Internal(_) | NodeContent::Leaf(_)=> self.occupied_8bit(target_child_key)
                                     & RAY_TO_NODE_OCCUPANCY_BITMASK_LUT[hash_region(
                                         &(ray.point_at(ray_current_distance) - target_bounds.min_position),
-                                        target_bounds.size as f32,
+                                        target_bounds.size,
                                     ) as usize]
                                     [direction_lut_index as usize],
                             })
@@ -389,6 +381,7 @@ impl<T: Default + PartialEq + Clone + std::fmt::Debug + VoxelData, const DIM: us
                         // stop advancing because current target is either
                         // - OOB
                         // - or (not empty while inside bounds AND collides with the ray based on its occupancy bitmap)
+                        node_stack.last_mut().unwrap().target_octant = target_octant;
                         break;
                     }
                     target_hit = target_bounds.intersect_ray(&ray);
diff --git a/src/octree/raytracing/tests.rs b/src/octree/raytracing/tests.rs
index 116b13d..0f0e603 100644
--- a/src/octree/raytracing/tests.rs
+++ b/src/octree/raytracing/tests.rs
@@ -1,3 +1,63 @@
+use crate::octree::{Cube, V3c};
+use crate::spatial::raytracing::Ray;
+use crate::spatial::{math::plane_line_intersection, FLOAT_ERROR_TOLERANCE};
+
+/// Reference implementation to decide step to sibling boundary
+#[allow(dead_code)]
+pub(crate) fn get_step_to_next_sibling(current: &Cube, ray: &Ray) -> V3c<f32> {
+    //Find the point furthest from the ray
+    let midpoint = V3c::unit((current.size / 2.0) as f32) + current.min_position.into();
+    let ref_point = midpoint
+        + V3c::new(
+            (current.size as f32 / 2.).copysign(ray.direction.x),
+            (current.size as f32 / 2.).copysign(ray.direction.y),
+            (current.size as f32 / 2.).copysign(ray.direction.z),
+        );
+
+    // Find the min of the 3 plane intersections
+    let x_plane_distance = plane_line_intersection(
+        &ref_point,
+        &V3c::new(1., 0., 0.),
+        &ray.origin,
+        &ray.direction,
+    )
+    .unwrap_or(f32::MAX);
+    let y_plane_distance = plane_line_intersection(
+        &ref_point,
+        &V3c::new(0., 1., 0.),
+        &ray.origin,
+        &ray.direction,
+    )
+    .unwrap_or(f32::MAX);
+    let z_plane_distance = plane_line_intersection(
+        &ref_point,
+        &V3c::new(0., 0., 1.),
+        &ray.origin,
+        &ray.direction,
+    )
+    .unwrap_or(f32::MAX);
+    let min_d = x_plane_distance.min(y_plane_distance).min(z_plane_distance);
+
+    // Step along the axes with the minimum distances
+    V3c::new(
+        if (min_d - x_plane_distance).abs() < FLOAT_ERROR_TOLERANCE {
+            (current.size as f32).copysign(ray.direction.x)
+        } else {
+            0.
+        },
+        if (min_d - y_plane_distance).abs() < FLOAT_ERROR_TOLERANCE {
+            (current.size as f32).copysign(ray.direction.y)
+        } else {
+            0.
+        },
+        if (min_d - z_plane_distance).abs() < FLOAT_ERROR_TOLERANCE {
+            (current.size as f32).copysign(ray.direction.z)
+        } else {
+            0.
+        },
+    )
+}
+
 #[cfg(test)]
 mod wgpu_tests {
     #[test]
@@ -9,67 +69,12 @@ mod wgpu_tests {
 
 #[cfg(test)]
 mod octree_raytracing_tests {
-    use crate::octree::{Albedo, Cube, Octree, V3c};
+    use crate::octree::{raytracing::tests::get_step_to_next_sibling, Albedo, Cube, Octree, V3c};
     use crate::spatial::raytracing::Ray;
-    use crate::spatial::{math::plane_line_intersection, FLOAT_ERROR_TOLERANCE};
+    use crate::spatial::FLOAT_ERROR_TOLERANCE;
 
     use rand::{rngs::ThreadRng, Rng};
 
-    /// Reference implementation to decide step to sibling boundary
-    fn get_step_to_next_sibling(current: &Cube, ray: &Ray) -> V3c<f32> {
-        //Find the point furthest from the ray
-        let midpoint = V3c::unit((current.size / 2.0) as f32) + current.min_position.into();
-        let ref_point = midpoint
-            + V3c::new(
-                (current.size as f32 / 2.).copysign(ray.direction.x),
-                (current.size as f32 / 2.).copysign(ray.direction.y),
-                (current.size as f32 / 2.).copysign(ray.direction.z),
-            );
-
-        // Find the min of the 3 plane intersections
-        let x_plane_distance = plane_line_intersection(
-            &ref_point,
-            &V3c::new(1., 0., 0.),
-            &ray.origin,
-            &ray.direction,
-        )
-        .unwrap_or(f32::MAX);
-        let y_plane_distance = plane_line_intersection(
-            &ref_point,
-            &V3c::new(0., 1., 0.),
-            &ray.origin,
-            &ray.direction,
-        )
-        .unwrap_or(f32::MAX);
-        let z_plane_distance = plane_line_intersection(
-            &ref_point,
-            &V3c::new(0., 0., 1.),
-            &ray.origin,
-            &ray.direction,
-        )
-        .unwrap_or(f32::MAX);
-        let min_d = x_plane_distance.min(y_plane_distance).min(z_plane_distance);
-
-        // Step along the axes with the minimum distances
-        V3c::new(
-            if (min_d - x_plane_distance).abs() < FLOAT_ERROR_TOLERANCE {
-                (current.size as f32).copysign(ray.direction.x)
-            } else {
-                0.
-            },
-            if (min_d - y_plane_distance).abs() < FLOAT_ERROR_TOLERANCE {
-                (current.size as f32).copysign(ray.direction.y)
-            } else {
-                0.
-            },
-            if (min_d - z_plane_distance).abs() < FLOAT_ERROR_TOLERANCE {
-                (current.size as f32).copysign(ray.direction.z)
-            } else {
-                0.
-            },
-        )
-    }
-
     #[test]
     #[ignore = "May fail in edge cases"]
     fn compare_sibling_step_functions() {
diff --git a/src/spatial/math/mod.rs b/src/spatial/math/mod.rs
index 1034d29..ff7b0fb 100644
--- a/src/spatial/math/mod.rs
+++ b/src/spatial/math/mod.rs
@@ -1,7 +1,7 @@
 mod tests;
 pub mod vector;
 
-use crate::spatial::{math::vector::V3c, Cube};
+use crate::spatial::{math::vector::V3c, raytracing::lut::OCTANT_STEP_RESULT_LUT, Cube};
 
 ///####################################################################################
 /// Octant
@@ -42,6 +42,16 @@ pub(crate) fn hash_direction(direction: &V3c<f32>) -> u8 {
     hash_region(&offset, 2.)
 }
 
+pub(crate) fn step_octant(octant: u8, step: V3c<f32>) -> u8 {
+    let step_signum_index = V3c::new(
+        ((step.x as i32).signum() + 1) as usize,
+        ((step.y as i32).signum() + 1) as usize,
+        ((step.z as i32).signum() + 1) as usize,
+    );
+    OCTANT_STEP_RESULT_LUT[octant as usize][step_signum_index.x][step_signum_index.y]
+        [step_signum_index.z]
+}
+
 /// Maps 3 dimensional space limited by `size` to 1 dimension
 /// This mapping function supposes that the coordinates are bound inside
 /// a cube, each dimension `size` long.
diff --git a/src/spatial/math/vector.rs b/src/spatial/math/vector.rs
index 03f7942..448572e 100644
--- a/src/spatial/math/vector.rs
+++ b/src/spatial/math/vector.rs
@@ -29,6 +29,13 @@ impl V3c<f32> {
     pub fn normalized(self) -> V3c<f32> {
         self / self.length()
     }
+    pub fn signum(&self) -> V3c<f32> {
+        V3c {
+            x: self.x.signum(),
+            y: self.y.signum(),
+            z: self.z.signum(),
+        }
+    }
 }
 
 impl V3c<u32> {
diff --git a/src/spatial/raytracing/lut.rs b/src/spatial/raytracing/lut.rs
index 0eefab4..a3eefe9 100644
--- a/src/spatial/raytracing/lut.rs
+++ b/src/spatial/raytracing/lut.rs
@@ -136,6 +136,146 @@ fn generate_lut_8_bits() -> [[u8; 8]; 8] {
     bitmap_lut
 }
 
+pub const OOB_OCTANT: u8 = 8;
+pub const OCTANT_STEP_RESULT_LUT: [[[[u8; 3]; 3]; 3]; 8] = [
+    [
+        [
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+        ],
+        [
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+            [OOB_OCTANT, 0, 2],
+            [OOB_OCTANT, 4, 6],
+        ],
+        [
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+            [OOB_OCTANT, 1, 3],
+            [OOB_OCTANT, 5, 7],
+        ],
+    ],
+    [
+        [
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+            [OOB_OCTANT, 0, 2],
+            [OOB_OCTANT, 4, 6],
+        ],
+        [
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+            [OOB_OCTANT, 1, 3],
+            [OOB_OCTANT, 5, 7],
+        ],
+        [
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+        ],
+    ],
+    [
+        [
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+        ],
+        [
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+            [0, 2, OOB_OCTANT],
+            [4, 6, OOB_OCTANT],
+        ],
+        [
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+            [1, 3, OOB_OCTANT],
+            [5, 7, OOB_OCTANT],
+        ],
+    ],
+    [
+        [
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+            [0, 2, OOB_OCTANT],
+            [4, 6, OOB_OCTANT],
+        ],
+        [
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+            [1, 3, OOB_OCTANT],
+            [5, 7, OOB_OCTANT],
+        ],
+        [
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+        ],
+    ],
+    [
+        [
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+        ],
+        [
+            [OOB_OCTANT, 0, 2],
+            [OOB_OCTANT, 4, 6],
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+        ],
+        [
+            [OOB_OCTANT, 1, 3],
+            [OOB_OCTANT, 5, 7],
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+        ],
+    ],
+    [
+        [
+            [OOB_OCTANT, 0, 2],
+            [OOB_OCTANT, 4, 6],
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+        ],
+        [
+            [OOB_OCTANT, 1, 3],
+            [OOB_OCTANT, 5, 7],
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+        ],
+        [
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+        ],
+    ],
+    [
+        [
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+        ],
+        [
+            [0, 2, OOB_OCTANT],
+            [4, 6, OOB_OCTANT],
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+        ],
+        [
+            [1, 3, OOB_OCTANT],
+            [5, 7, OOB_OCTANT],
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+        ],
+    ],
+    [
+        [
+            [0, 2, OOB_OCTANT],
+            [4, 6, OOB_OCTANT],
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+        ],
+        [
+            [1, 3, OOB_OCTANT],
+            [5, 7, OOB_OCTANT],
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+        ],
+        [
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+            [OOB_OCTANT, OOB_OCTANT, OOB_OCTANT],
+        ],
+    ],
+];
+
 pub const RAY_TO_NODE_OCCUPANCY_BITMASK_LUT: [[u8; 8]; 8] = [
     [1, 3, 5, 15, 17, 51, 85, 255],
     [3, 2, 15, 10, 51, 34, 255, 170],