From 063c487ad24c85cc1922cca58294285a8c5f566e Mon Sep 17 00:00:00 2001
From: Kyle Grover <kyle@kylegrover.com>
Date: Sun, 16 Nov 2025 22:16:24 -0800
Subject: [PATCH 01/10] modify voxels in place to avoid mem copy

---
 app/geometry/sphere.py | 15 +++++++--------
 1 file changed, 7 insertions(+), 8 deletions(-)

diff --git a/app/geometry/sphere.py b/app/geometry/sphere.py
index efcf00a..e945971 100644
--- a/app/geometry/sphere.py
+++ b/app/geometry/sphere.py
@@ -116,25 +116,24 @@ def _maybe_expand_voxel_space(self, voxel_space, max_index, axis_number):
         return voxel_space
 
     def _deposit_sphere(self, radius, voxel_space, nozzle_height, target_volume):
-        copy_voxel_space = voxel_space.copy()
-
-        copy_voxel_space = self.deform_voxel_space_for_big_spheres(
-            copy_voxel_space, radius
+        # Modify voxel_space in-place to avoid redundant copying
+        voxel_space = self.deform_voxel_space_for_big_spheres(
+            voxel_space, radius
         )
 
         lower_indexes, upper_indexes = self.find_sphere_limits(radius, nozzle_height)
 
-        copy_voxel_space = self.fill_voxels(
-            copy_voxel_space, radius, lower_indexes, upper_indexes
+        voxel_space = self.fill_voxels(
+            voxel_space, radius, lower_indexes, upper_indexes
         )
 
         current_volume = GeometryMath.calculate_filled_volume(
-            copy_voxel_space, self.voxel_size
+            voxel_space, self.voxel_size
         )
 
         volume_overshoot = current_volume / target_volume - 1.0
 
-        return volume_overshoot, copy_voxel_space
+        return volume_overshoot, voxel_space
 
     def _increase_solver_tolerance(self, radius_a, radius_b):
         return radius_b - radius_a < self.voxel_size * 0.5

From 5707e139c89f2e26bbdcfcb7229e3459e8ec3fb7 Mon Sep 17 00:00:00 2001
From: Kyle Grover <kyle@kylegrover.com>
Date: Sun, 16 Nov 2025 22:18:48 -0800
Subject: [PATCH 02/10] vectorize voxel fill operations using numpy

---
 app/geometry/sphere.py | 21 ++++++++++++---------
 1 file changed, 12 insertions(+), 9 deletions(-)

diff --git a/app/geometry/sphere.py b/app/geometry/sphere.py
index e945971..3be99f0 100644
--- a/app/geometry/sphere.py
+++ b/app/geometry/sphere.py
@@ -41,15 +41,18 @@ def fill_voxels(self, voxel_space, radius, lower_indexes, upper_indexes):
             voxel_space, lower_indexes, upper_indexes
         )
 
-        for voxel in empty_voxels:
-            voxel_coordinate = GeometryMath.find_coordinates(voxel, self.voxel_size)
-
-            distance_to_centre = GeometryMath.distance(
-                voxel_coordinate, self.centre_coordinates
-            )
-
-            if distance_to_centre <= radius + self.voxel_size * 1e-8:
-                voxel_space[tuple(voxel)] = 1
+        # Convert to numpy array for vectorized operations
+        empty_voxels_np = np.array(empty_voxels)
+        # Calculate coordinates for all voxels at once
+        voxel_coords = self.voxel_size * (2 * (empty_voxels_np + 1) - 1) * 0.5
+        # Calculate squared distances to centre for all voxels
+        centre = np.array(self.centre_coordinates)
+        dists = np.linalg.norm(voxel_coords - centre, axis=1)
+        # Mask for voxels within radius
+        mask = dists <= radius + self.voxel_size * 1e-8
+        # Assign value 1 to all selected voxels
+        for idx in empty_voxels_np[mask]:
+            voxel_space[tuple(idx)] = 1
 
         return voxel_space
 

From 072267cff9dce0597e5bd67b70638c732414d47a Mon Sep 17 00:00:00 2001
From: Kyle Grover <kyle@kylegrover.com>
Date: Sun, 16 Nov 2025 22:43:42 -0800
Subject: [PATCH 03/10] check for no empty voxels

---
 app/geometry/sphere.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/app/geometry/sphere.py b/app/geometry/sphere.py
index 3be99f0..ba451d2 100644
--- a/app/geometry/sphere.py
+++ b/app/geometry/sphere.py
@@ -41,6 +41,9 @@ def fill_voxels(self, voxel_space, radius, lower_indexes, upper_indexes):
             voxel_space, lower_indexes, upper_indexes
         )
 
+        if not empty_voxels:
+            return voxel_space
+
         # Convert to numpy array for vectorized operations
         empty_voxels_np = np.array(empty_voxels)
         # Calculate coordinates for all voxels at once

From 4bf5918e6acca381c961c649e922776489bc9793 Mon Sep 17 00:00:00 2001
From: Kyle Grover <kyle@kylegrover.com>
Date: Sun, 16 Nov 2025 23:12:27 -0800
Subject: [PATCH 04/10] feat: enhance voxel space expansion logging and update
 bisection method to handle updated arguments

---
 app/geometry/sphere.py          |  9 +++++----
 app/solvers/bisection_method.py | 35 ++++++++++++++++++++-------------
 2 files changed, 26 insertions(+), 18 deletions(-)

diff --git a/app/geometry/sphere.py b/app/geometry/sphere.py
index ba451d2..6091038 100644
--- a/app/geometry/sphere.py
+++ b/app/geometry/sphere.py
@@ -111,11 +111,12 @@ def _maybe_expand_voxel_space(self, voxel_space, max_index, axis_number):
             return voxel_space
 
         number_to_be_added = max_index - index_size + 1
+        import logging
+        logging.info(f"Expanding voxel space on axis {axis_number}: current size {index_size}, required {max_index}, adding {number_to_be_added}")
 
-        new_size = list(size)
-        new_size[axis_number] = number_to_be_added
-
-        mat_add = np.zeros(new_size, dtype=np.int8)
+        mat_add_size = list(size)
+        mat_add_size[axis_number] = number_to_be_added
+        mat_add = np.zeros(mat_add_size, dtype=np.int8)
 
         voxel_space = np.concatenate((voxel_space, mat_add), axis=axis_number)
 
diff --git a/app/solvers/bisection_method.py b/app/solvers/bisection_method.py
index 953ce58..5bfbaa9 100644
--- a/app/solvers/bisection_method.py
+++ b/app/solvers/bisection_method.py
@@ -23,35 +23,42 @@ def execute(
         else:
             point_b = initial_point
 
-        _, _, point_a, point_b = self._loop_fb(fun, point_b, increment, args)
+        fb, out, point_a, point_b, updated_args = self._loop_fb(fun, point_b, increment, args)
 
-        return self._loop_fc(
-            fun, point_a, point_b, tolerance, fun_increase_tolerance, args
+        # _loop_fc returns (fc, out) where out may be the raw output from `fun`.
+        fc, out = self._loop_fc(
+            fun, point_a, point_b, tolerance, fun_increase_tolerance, updated_args
         )
 
+        # If `out` is a tuple like (volume_overshoot, voxel_space), unwrap it
+        if isinstance(out, tuple) and len(out) == 2:
+            return fc, out[1]
+
+        return fc, out
+
     def _loop_fb(self, fun, point_b, inc, args):
         fb = -1.0
-
         point_a = 0.0
-
+        updated_args = list(args)
         while fb < 0.0:
-
-            fb, out = fun(point_b, *args)
-
+            fb, out = fun(point_b, *updated_args)
+            # If the function returns an updated voxel_space as the second value,
+            # `out` will be that ndarray (it will have a `shape` attribute).
+            if hasattr(out, "shape"):
+                updated_args[0] = out
             if fb < 0:
                 point_a = point_b
                 point_b += inc
-
-        return fb, out, point_a, point_b
+        return fb, out, point_a, point_b, updated_args
 
     def _loop_fc(self, fun, point_a, point_b, tolerance, fun_increase_tolerance, args):
         fc = 2.0 * tolerance
-
+        updated_args = list(args)
         while abs(fc) > tolerance:
             point_c = (point_a + point_b) * 0.5
-
-            fc, out = fun(point_c, *args)
-
+            fc, out = fun(point_c, *updated_args)
+            if hasattr(out, "shape"):
+                updated_args[0] = out
             if fun_increase_tolerance(point_a, point_b):
                 logger.debug(
                     "[BisectionMethod]: increasing tolerance because point_b and point_a are too close"

From ffd29cb7842f750871e25aa3f363af4930f6dcdf Mon Sep 17 00:00:00 2001
From: Kyle Grover <kyle@kylegrover.com>
Date: Sun, 16 Nov 2025 23:24:24 -0800
Subject: [PATCH 05/10] preallocate voxel space

---
 app/geometry/sphere.py      |  9 +++++++--
 app/geometry/voxel_space.py | 40 ++++++++++++++++++++++++++++++++++++-
 2 files changed, 46 insertions(+), 3 deletions(-)

diff --git a/app/geometry/sphere.py b/app/geometry/sphere.py
index 6091038..94554e7 100644
--- a/app/geometry/sphere.py
+++ b/app/geometry/sphere.py
@@ -112,10 +112,15 @@ def _maybe_expand_voxel_space(self, voxel_space, max_index, axis_number):
 
         number_to_be_added = max_index - index_size + 1
         import logging
-        logging.info(f"Expanding voxel space on axis {axis_number}: current size {index_size}, required {max_index}, adding {number_to_be_added}")
 
+        # Add a buffer to reduce the number of reallocations. Buffer is 20%
+        # of current size (rounded), but at least the minimum required.
+        buffer_layers = max(int(index_size * 0.2), 1)
+        layers_to_add = max(number_to_be_added, buffer_layers)
+
+        # Create only the additional block to concatenate
         mat_add_size = list(size)
-        mat_add_size[axis_number] = number_to_be_added
+        mat_add_size[axis_number] = layers_to_add
         mat_add = np.zeros(mat_add_size, dtype=np.int8)
 
         voxel_space = np.concatenate((voxel_space, mat_add), axis=axis_number)
diff --git a/app/geometry/voxel_space.py b/app/geometry/voxel_space.py
index 819d576..2bb0f0a 100644
--- a/app/geometry/voxel_space.py
+++ b/app/geometry/voxel_space.py
@@ -47,8 +47,46 @@ def __init__(
         self._consider_acceleration = self._simulation.consider_acceleration
 
     def initialize_space(self):
+        # Try to intelligently preallocate Z dimension to avoid repeated expansions.
+        # Estimate the maximum per-step sphere radius from filaments and add a safety margin.
+        z_dim = self.dimensions["z"]
+
+        try:
+            max_radius = 0.0
+            step_size = self._simulation.step_size
+            for filament in self._instruction.filaments_coordinates:
+                coord_old = filament[0]
+                coord_new = filament[1]
+                volume = filament[2]
+
+                # compute filament length and estimated number of steps
+                length = GeometryMath.distance(coord_old, coord_new)
+                n_steps = int(round(length / step_size))
+                if n_steps < 1:
+                    n_steps = 1
+
+                per_step_volume = volume / n_steps
+                # estimate radius of sphere that would contain this volume
+                radius = (3.0 * per_step_volume / (4.0 * math.pi)) ** (1.0 / 3.0)
+                if radius > max_radius:
+                    max_radius = radius
+
+            # include configured sphere z offset
+            max_extra_z = int(math.ceil((max_radius + self._simulation.sphere_z_offset) / self._simulation.voxel_size))
+
+            # safety margin (20% of current z) to reduce chance of further expansion
+            safety_margin = int(max(1, z_dim * 0.2))
+
+            if max_extra_z > 0:
+                z_dim = z_dim + max_extra_z + safety_margin
+
+            logging.getLogger(__name__).info(f"Preallocating voxel space z-dimension: base={self.dimensions['z']}, extra={max_extra_z}, safety={safety_margin}, final={z_dim}")
+        except Exception:
+            # Fallback to original allocation if any issue occurs during estimation
+            z_dim = self.dimensions["z"]
+
         self.space = np.zeros(
-            (self.dimensions["x"], self.dimensions["y"], self.dimensions["z"]),
+            (self.dimensions["x"], self.dimensions["y"], z_dim),
             dtype=np.int8,
         )
 

From 426b2754dee74ec576afa122f8f43f62912c87a2 Mon Sep 17 00:00:00 2001
From: Kyle Grover <kyle@kylegrover.com>
Date: Sun, 16 Nov 2025 23:27:43 -0800
Subject: [PATCH 06/10] keep track of filled voxels

---
 app/geometry/geometry_math.py   | 15 +++++++
 app/geometry/sphere.py          | 72 +++++++++++++++++++++------------
 app/geometry/voxel_space.py     | 22 ++++++++--
 app/solvers/bisection_method.py |  9 +++--
 4 files changed, 84 insertions(+), 34 deletions(-)

diff --git a/app/geometry/geometry_math.py b/app/geometry/geometry_math.py
index b0dc014..4855fcf 100644
--- a/app/geometry/geometry_math.py
+++ b/app/geometry/geometry_math.py
@@ -49,6 +49,21 @@ def find_coordinates(indexes, voxel_size):
 
     @staticmethod
     def calculate_filled_volume(voxel_space, voxel_size):
+        # Accept either a raw ndarray or a VoxelSpace-like object with a
+        # `_filled_voxels_count` running counter. Prefer the running counter
+        # when available to avoid expensive full-array scans.
+        if hasattr(voxel_space, "_filled_voxels_count"):
+            return int(voxel_space._filled_voxels_count) * voxel_size ** 3
+
+        # If an object with `.space` is provided, try to use its counter first.
+        if hasattr(voxel_space, "space") and hasattr(voxel_space.space, "shape"):
+            # If the container itself tracks a counter, use it.
+            if hasattr(voxel_space, "_filled_voxels_count"):
+                return int(voxel_space._filled_voxels_count) * voxel_size ** 3
+            # Fallback to scanning the ndarray
+            return np.count_nonzero(voxel_space.space) * voxel_size ** 3
+
+        # Otherwise assume voxel_space is a numpy array
         return np.count_nonzero(voxel_space) * voxel_size**3
 
     def find_index(coordinate, voxel_size):
diff --git a/app/geometry/sphere.py b/app/geometry/sphere.py
index 94554e7..37b4828 100644
--- a/app/geometry/sphere.py
+++ b/app/geometry/sphere.py
@@ -21,7 +21,10 @@ def deposit_sphere(
     ):
         initial_radius = self.estimate_initial_radius(sphere_volume)
 
-        _, voxel_space = BisectionMethod().execute(
+        # `voxel_space` is expected to be a VoxelSpace object here. The
+        # bisection solver and inner functions will mutate and ultimately
+        # return that object so callers can read `.space` and counters.
+        _, voxel_space_out = BisectionMethod().execute(
             self._deposit_sphere,
             initial_point=initial_radius,
             tolerance=solver_tolerance,
@@ -34,30 +37,51 @@ def deposit_sphere(
             ),
         )
 
-        return voxel_space
+        return voxel_space_out
 
-    def fill_voxels(self, voxel_space, radius, lower_indexes, upper_indexes):
+    def fill_voxels(self, voxel_space_obj, radius, lower_indexes, upper_indexes):
+        # `voxel_space_obj` is a VoxelSpace instance. Operate on its `.space` ndarray
         empty_voxels = GeometryMath.find_empty_voxels_in_space(
-            voxel_space, lower_indexes, upper_indexes
+            voxel_space_obj.space, lower_indexes, upper_indexes
         )
 
         if not empty_voxels:
-            return voxel_space
+            return voxel_space_obj
 
         # Convert to numpy array for vectorized operations
         empty_voxels_np = np.array(empty_voxels)
         # Calculate coordinates for all voxels at once
         voxel_coords = self.voxel_size * (2 * (empty_voxels_np + 1) - 1) * 0.5
-        # Calculate squared distances to centre for all voxels
+        # Calculate distances to centre for all voxels
         centre = np.array(self.centre_coordinates)
         dists = np.linalg.norm(voxel_coords - centre, axis=1)
         # Mask for voxels within radius
         mask = dists <= radius + self.voxel_size * 1e-8
-        # Assign value 1 to all selected voxels
-        for idx in empty_voxels_np[mask]:
-            voxel_space[tuple(idx)] = 1
 
-        return voxel_space
+        if not np.any(mask):
+            return voxel_space_obj
+
+        # Filter indices that should be filled
+        fill_indices = empty_voxels_np[mask]
+
+        # Advanced index assignment (vectorized)
+        xi = fill_indices[:, 0].astype(int)
+        yj = fill_indices[:, 1].astype(int)
+        zk = fill_indices[:, 2].astype(int)
+
+        # Before setting, count how many of these are actually zero (defensive)
+        # They should be zero because `find_empty_voxels_in_space` returned empties,
+        # but reconfirm to be robust in case of race or prior modifications.
+        current_vals = voxel_space_obj.space[xi, yj, zk]
+        new_mask = current_vals == 0
+        n_new = int(np.count_nonzero(new_mask))
+        if n_new > 0:
+            voxel_space_obj.space[xi[new_mask], yj[new_mask], zk[new_mask]] = 1
+            # Update running counter on the VoxelSpace object
+            if hasattr(voxel_space_obj, "_filled_voxels_count"):
+                voxel_space_obj._filled_voxels_count += n_new
+
+        return voxel_space_obj
 
     def estimate_initial_radius(self, volume):
         return (3.0 * volume / (4.0 * math.pi)) ** (1.0 / 3.0)
@@ -91,6 +115,7 @@ def find_sphere_limits(self, radius, nozzle_height):
     """
 
     def deform_voxel_space_for_big_spheres(self, voxel_space, radius):
+        # voxel_space is expected to be a VoxelSpace instance; check and operate on its `.space`
         max_indexes = [
             self._find_index(coord + radius) for coord in self.centre_coordinates
         ]
@@ -102,16 +127,16 @@ def deform_voxel_space_for_big_spheres(self, voxel_space, radius):
 
         return voxel_space
 
-    def _maybe_expand_voxel_space(self, voxel_space, max_index, axis_number):
-        size = voxel_space.shape
+    def _maybe_expand_voxel_space(self, voxel_space_obj, max_index, axis_number):
+        # Operate on `voxel_space_obj.space` and update in-place by reassigning.
+        size = voxel_space_obj.space.shape
 
         index_size = size[axis_number]
 
         if max_index < index_size:
-            return voxel_space
+            return voxel_space_obj
 
         number_to_be_added = max_index - index_size + 1
-        import logging
 
         # Add a buffer to reduce the number of reallocations. Buffer is 20%
         # of current size (rounded), but at least the minimum required.
@@ -123,28 +148,23 @@ def _maybe_expand_voxel_space(self, voxel_space, max_index, axis_number):
         mat_add_size[axis_number] = layers_to_add
         mat_add = np.zeros(mat_add_size, dtype=np.int8)
 
-        voxel_space = np.concatenate((voxel_space, mat_add), axis=axis_number)
+        voxel_space_obj.space = np.concatenate((voxel_space_obj.space, mat_add), axis=axis_number)
 
-        return voxel_space
+        return voxel_space_obj
 
     def _deposit_sphere(self, radius, voxel_space, nozzle_height, target_volume):
-        # Modify voxel_space in-place to avoid redundant copying
-        voxel_space = self.deform_voxel_space_for_big_spheres(
-            voxel_space, radius
-        )
+        # `voxel_space` is a VoxelSpace instance; deform and fill it in-place.
+        voxel_space = self.deform_voxel_space_for_big_spheres(voxel_space, radius)
 
         lower_indexes, upper_indexes = self.find_sphere_limits(radius, nozzle_height)
 
-        voxel_space = self.fill_voxels(
-            voxel_space, radius, lower_indexes, upper_indexes
-        )
+        voxel_space = self.fill_voxels(voxel_space, radius, lower_indexes, upper_indexes)
 
-        current_volume = GeometryMath.calculate_filled_volume(
-            voxel_space, self.voxel_size
-        )
+        current_volume = GeometryMath.calculate_filled_volume(voxel_space, self.voxel_size)
 
         volume_overshoot = current_volume / target_volume - 1.0
 
+        # Return the computed overshoot and the (possibly mutated) VoxelSpace object
         return volume_overshoot, voxel_space
 
     def _increase_solver_tolerance(self, radius_a, radius_b):
diff --git a/app/geometry/voxel_space.py b/app/geometry/voxel_space.py
index 2bb0f0a..d4c75cf 100644
--- a/app/geometry/voxel_space.py
+++ b/app/geometry/voxel_space.py
@@ -45,6 +45,8 @@ def __init__(
         self._simulation = simulation_config
         self._printer = printer
         self._consider_acceleration = self._simulation.consider_acceleration
+        # Running count of filled voxels (to avoid repeated `np.count_nonzero` calls)
+        self._filled_voxels_count = 0
 
     def initialize_space(self):
         # Try to intelligently preallocate Z dimension to avoid repeated expansions.
@@ -89,6 +91,8 @@ def initialize_space(self):
             (self.dimensions["x"], self.dimensions["y"], z_dim),
             dtype=np.int8,
         )
+        # reset the running counter when allocating space
+        self._filled_voxels_count = 0
 
     def print(self):
         number_printed_filaments = 0
@@ -173,8 +177,9 @@ def _deposit_filament(
         filament_initial_coordinates,
         volumes,
     ):
-        total_deposited_volume = GeometryMath.calculate_filled_volume(
-            self.space, self._simulation.voxel_size
+        # Use the running counter to compute total deposited volume quickly
+        total_deposited_volume = (
+            self._filled_voxels_count * self._simulation.voxel_size ** 3
         )
 
         for step_n in range(0, number_simulation_steps):
@@ -221,11 +226,20 @@ def _deposit_filament(
                 f"Depositing filament: step = {step_n + 1}/{number_simulation_steps}"
             )
 
-            self.space = sphere.deposit_sphere(
-                voxel_space=self.space,
+            # Pass the VoxelSpace instance so sphere code can update the
+            # running counter and expand the `.space` ndarray in-place.
+            voxel_space_out = sphere.deposit_sphere(
+                voxel_space=self,
                 nozzle_height=nozzle_height,
                 sphere_volume=sphere_volume,
                 voxel_space_target_volume=volume_target,
                 solver_tolerance=self._simulation.solver_tolerance,
                 radius_increment=self._simulation.radius_increment,
             )
+
+            # Ensure our internal `.space` and counter reflect any mutations
+            # returned by the sphere logic (voxel_space_out is a VoxelSpace)
+            if hasattr(voxel_space_out, "space"):
+                self.space = voxel_space_out.space
+            if hasattr(voxel_space_out, "_filled_voxels_count"):
+                self._filled_voxels_count = voxel_space_out._filled_voxels_count
diff --git a/app/solvers/bisection_method.py b/app/solvers/bisection_method.py
index 5bfbaa9..80b5e09 100644
--- a/app/solvers/bisection_method.py
+++ b/app/solvers/bisection_method.py
@@ -42,9 +42,10 @@ def _loop_fb(self, fun, point_b, inc, args):
         updated_args = list(args)
         while fb < 0.0:
             fb, out = fun(point_b, *updated_args)
-            # If the function returns an updated voxel_space as the second value,
-            # `out` will be that ndarray (it will have a `shape` attribute).
-            if hasattr(out, "shape"):
+            # If the function returns an updated voxel_space, propagate it.
+            # The returned `out` may be an ndarray (has `shape`) or a
+            # VoxelSpace-like object with a `space` attribute. Accept both.
+            if hasattr(out, "shape") or hasattr(out, "space"):
                 updated_args[0] = out
             if fb < 0:
                 point_a = point_b
@@ -57,7 +58,7 @@ def _loop_fc(self, fun, point_a, point_b, tolerance, fun_increase_tolerance, arg
         while abs(fc) > tolerance:
             point_c = (point_a + point_b) * 0.5
             fc, out = fun(point_c, *updated_args)
-            if hasattr(out, "shape"):
+            if hasattr(out, "shape") or hasattr(out, "space"):
                 updated_args[0] = out
             if fun_increase_tolerance(point_a, point_b):
                 logger.debug(

From 47d7438a28ab2edfe6357701b9cd2b5873393e7b Mon Sep 17 00:00:00 2001
From: Kyle Grover <kyle@kylegrover.com>
Date: Mon, 17 Nov 2025 07:40:13 -0800
Subject: [PATCH 07/10] fix py tests

---
 app/geometry/sphere.py             | 58 +++++++++++++++++-------------
 app/physics/volume.py              |  2 +-
 tests/geometry/voxel_space_test.py | 11 ++++--
 3 files changed, 44 insertions(+), 27 deletions(-)

diff --git a/app/geometry/sphere.py b/app/geometry/sphere.py
index 37b4828..a3e181f 100644
--- a/app/geometry/sphere.py
+++ b/app/geometry/sphere.py
@@ -40,9 +40,12 @@ def deposit_sphere(
         return voxel_space_out
 
     def fill_voxels(self, voxel_space_obj, radius, lower_indexes, upper_indexes):
-        # `voxel_space_obj` is a VoxelSpace instance. Operate on its `.space` ndarray
+        # Accept either a VoxelSpace-like object (has `.space`) or a raw ndarray.
+        is_voxel_space = hasattr(voxel_space_obj, "space")
+        space = voxel_space_obj.space if is_voxel_space else voxel_space_obj
+
         empty_voxels = GeometryMath.find_empty_voxels_in_space(
-            voxel_space_obj.space, lower_indexes, upper_indexes
+            space, lower_indexes, upper_indexes
         )
 
         if not empty_voxels:
@@ -70,18 +73,20 @@ def fill_voxels(self, voxel_space_obj, radius, lower_indexes, upper_indexes):
         zk = fill_indices[:, 2].astype(int)
 
         # Before setting, count how many of these are actually zero (defensive)
-        # They should be zero because `find_empty_voxels_in_space` returned empties,
-        # but reconfirm to be robust in case of race or prior modifications.
-        current_vals = voxel_space_obj.space[xi, yj, zk]
+        current_vals = space[xi, yj, zk]
         new_mask = current_vals == 0
         n_new = int(np.count_nonzero(new_mask))
         if n_new > 0:
-            voxel_space_obj.space[xi[new_mask], yj[new_mask], zk[new_mask]] = 1
-            # Update running counter on the VoxelSpace object
-            if hasattr(voxel_space_obj, "_filled_voxels_count"):
+            # Assign into the ndarray `space` (in-place)
+            space[xi[new_mask], yj[new_mask], zk[new_mask]] = 1
+            # Update running counter only when we were given a VoxelSpace object
+            if is_voxel_space and hasattr(voxel_space_obj, "_filled_voxels_count"):
                 voxel_space_obj._filled_voxels_count += n_new
 
-        return voxel_space_obj
+        # Return the same type we were given
+        if is_voxel_space:
+            return voxel_space_obj
+        return space
 
     def estimate_initial_radius(self, volume):
         return (3.0 * volume / (4.0 * math.pi)) ** (1.0 / 3.0)
@@ -115,22 +120,23 @@ def find_sphere_limits(self, radius, nozzle_height):
     """
 
     def deform_voxel_space_for_big_spheres(self, voxel_space, radius):
-        # voxel_space is expected to be a VoxelSpace instance; check and operate on its `.space`
+        # Accept either a VoxelSpace object or a raw ndarray and return the same type.
         max_indexes = [
             self._find_index(coord + radius) for coord in self.centre_coordinates
         ]
 
+        vs = voxel_space
         for axis_number in range(0, 3):
-            voxel_space = self._maybe_expand_voxel_space(
-                voxel_space, max_indexes[axis_number], axis_number
-            )
+            vs = self._maybe_expand_voxel_space(vs, max_indexes[axis_number], axis_number)
 
-        return voxel_space
+        return vs
 
     def _maybe_expand_voxel_space(self, voxel_space_obj, max_index, axis_number):
-        # Operate on `voxel_space_obj.space` and update in-place by reassigning.
-        size = voxel_space_obj.space.shape
+        # Accept either a VoxelSpace-like object (has `.space`) or a raw ndarray.
+        is_voxel_space = hasattr(voxel_space_obj, "space")
+        space = voxel_space_obj.space if is_voxel_space else voxel_space_obj
 
+        size = space.shape
         index_size = size[axis_number]
 
         if max_index < index_size:
@@ -148,24 +154,28 @@ def _maybe_expand_voxel_space(self, voxel_space_obj, max_index, axis_number):
         mat_add_size[axis_number] = layers_to_add
         mat_add = np.zeros(mat_add_size, dtype=np.int8)
 
-        voxel_space_obj.space = np.concatenate((voxel_space_obj.space, mat_add), axis=axis_number)
+        new_space = np.concatenate((space, mat_add), axis=axis_number)
+
+        if is_voxel_space:
+            voxel_space_obj.space = new_space
+            return voxel_space_obj
 
-        return voxel_space_obj
+        return new_space
 
     def _deposit_sphere(self, radius, voxel_space, nozzle_height, target_volume):
-        # `voxel_space` is a VoxelSpace instance; deform and fill it in-place.
-        voxel_space = self.deform_voxel_space_for_big_spheres(voxel_space, radius)
+        # Accept and return either a VoxelSpace object or a raw ndarray.
+        vs = self.deform_voxel_space_for_big_spheres(voxel_space, radius)
 
         lower_indexes, upper_indexes = self.find_sphere_limits(radius, nozzle_height)
 
-        voxel_space = self.fill_voxels(voxel_space, radius, lower_indexes, upper_indexes)
+        vs = self.fill_voxels(vs, radius, lower_indexes, upper_indexes)
 
-        current_volume = GeometryMath.calculate_filled_volume(voxel_space, self.voxel_size)
+        current_volume = GeometryMath.calculate_filled_volume(vs, self.voxel_size)
 
         volume_overshoot = current_volume / target_volume - 1.0
 
-        # Return the computed overshoot and the (possibly mutated) VoxelSpace object
-        return volume_overshoot, voxel_space
+        # Return the computed overshoot and the (possibly mutated) voxel container
+        return volume_overshoot, vs
 
     def _increase_solver_tolerance(self, radius_a, radius_b):
         return radius_b - radius_a < self.voxel_size * 0.5
diff --git a/app/physics/volume.py b/app/physics/volume.py
index 977a788..4ca2005 100644
--- a/app/physics/volume.py
+++ b/app/physics/volume.py
@@ -27,7 +27,7 @@ def get_volumes_for_filament(
         if not consider_acceleration:
             # Simple case: evenly distribute the volume across steps
             volume_per_step = total_volume / number_simulation_steps
-            return [volume_per_step * (i + 1) for i in range(number_simulation_steps)]
+            return [volume_per_step for i in range(number_simulation_steps)]
         
         # Complex case: use acceleration-specific implementation
         from app.physics.acceleration.volume import AccelerationVolume
diff --git a/tests/geometry/voxel_space_test.py b/tests/geometry/voxel_space_test.py
index d6c3895..7bfa546 100644
--- a/tests/geometry/voxel_space_test.py
+++ b/tests/geometry/voxel_space_test.py
@@ -13,6 +13,10 @@ def __init__(self):
         self._number_printed_filaments = 0
         self._filaments_coordinates = list()
         self._default_nozzle_speed = 10.0
+        # Patch: provide a _printer attribute for compatibility
+        class _Printer:
+            feedstock_filament_diameter = 1.75
+        self._printer = _Printer()
 
     def read(self):
         pass
@@ -47,6 +51,7 @@ def __init__(self):
         self.step_size = 0.2
         self.radius_increment = 0.1
         self.solver_tolerance = 0.001
+        self.consider_acceleration = False
 
 
 class FakePrinter:
@@ -114,17 +119,19 @@ def __init__(self):
         self.step_size = 0.2
         self.radius_increment = 0.1
         self.solver_tolerance = 0.001
+        self.consider_acceleration = False
+        self.sphere_z_offset = 0.0
 
 
 class TestPrint:
     def test_should_print(self):
+        printer = FakePrinter()
         test_instruction = Gcode(
-            gcode_path="tests/fixtures/gcode_example.gcode", default_nozzle_speed=40.0
+            gcode_path="tests/fixtures/gcode_example.gcode", default_nozzle_speed=40.0, printer=printer
         )
         test_instruction.read()
 
         simulation_config = FakeSimulationPrint()
-        printer = FakePrinter()
 
         voxel_space = VoxelSpace(
             instruction=test_instruction,

From d5e92f720d43c0afeb439e06458856d3a15a6b3b Mon Sep 17 00:00:00 2001
From: Kyle Grover <kyle@kylegrover.com>
Date: Mon, 17 Nov 2025 19:09:31 -0800
Subject: [PATCH 08/10] add a lightweight preview mode that doesn't do any
 volume conservation available with a preview_mode: true or by passing a
 --preview flag

---
 README.md                   |  7 ++++
 app/configs/arguments.py    |  1 +
 app/configs/simulation.py   |  3 ++
 app/geometry/voxel_space.py | 66 ++++++++++++++++++++++++++++++-------
 volco.py                    | 14 +++++++-
 5 files changed, 78 insertions(+), 13 deletions(-)

diff --git a/README.md b/README.md
index e060843..0f24b2b 100644
--- a/README.md
+++ b/README.md
@@ -134,6 +134,7 @@ VOLCO uses two configuration files: simulation settings and printer settings. Be
 | `solver_tolerance` | Tolerance for volume conservation in the bisection method. | 0.0001 |
 | `consider_acceleration` | Whether to consider acceleration in volume distribution. | false |
 | `stl_ascii` | Whether to export STL in ASCII format (true) or binary (false). | false |
+| `preview_mode` | If true, runs a fast, lightweight preview that traces the G-code path and fills extrusions as capsules. No physics or overlap checks. Intended for quick feedback before running the full simulation. | false |
 
 ### Printer Configuration
 
@@ -156,6 +157,12 @@ VOLCO uses two configuration files: simulation settings and printer settings. Be
 
 - **consider_acceleration**: When true, the simulation accounts for acceleration and deceleration, which can provide more accurate results but increases computation time.
 
+## Preview Mode
+
+Preview mode is designed for fast, lightweight visualization of the print path. It traces the G-code path and fills extrusions as capsules (cylinders with rounded ends), matching the nozzle diameter. This mode skips all physics and overlap checks, making it ideal for quick feedback and rapid iteration before running the full, volume-conserving simulation.
+
+If you want to run VOLCO in maximum speed mode (no volume conservation, no overlap checks), set `max_speed_mode: true` in your simulation config. This will trace the G-code path and fill voxels along it, suitable for fast preview, parallel, or GPU-optimized runs.
+
 ## Finite Element Analysis (FEA)
 
 VOLCO includes a Finite Element Analysis (FEA) module that enables structural analysis of the simulated 3D printed parts. With just one line of code, you can analyze the structural behavior of your VOLCO simulation results:
diff --git a/app/configs/arguments.py b/app/configs/arguments.py
index c18f654..a2c999e 100644
--- a/app/configs/arguments.py
+++ b/app/configs/arguments.py
@@ -8,5 +8,6 @@ def get_options():
         parser.add_argument("--gcode", type=str)
         parser.add_argument("--sim", type=str)
         parser.add_argument("--printer", type=str)
+        parser.add_argument("--preview", action="store_true", help="Enable preview mode (fast, lightweight visualization)")
 
         return parser.parse_args()
diff --git a/app/configs/simulation.py b/app/configs/simulation.py
index 63a5832..304f82e 100644
--- a/app/configs/simulation.py
+++ b/app/configs/simulation.py
@@ -53,3 +53,6 @@ def _load_config_from_dict(self, config):
         # Default acceleration and STL settings
         self.consider_acceleration = config.get("consider_acceleration", False)
         self.stl_ascii = config.get("stl_ascii", False)
+
+        # Preview mode (fast, lightweight visualization)
+        self.preview_mode = config.get("preview_mode", False)
diff --git a/app/geometry/voxel_space.py b/app/geometry/voxel_space.py
index d4c75cf..f0d37bb 100644
--- a/app/geometry/voxel_space.py
+++ b/app/geometry/voxel_space.py
@@ -95,41 +95,35 @@ def initialize_space(self):
         self._filled_voxels_count = 0
 
     def print(self):
+        # Check for preview mode in simulation config
+        if hasattr(self._simulation, "preview_mode") and self._simulation.preview_mode:
+            self.print_preview_mode()
+            return
+
         number_printed_filaments = 0
         number_printed_layers = 0
-
         initial_z_coordinate = 0.0
-
         print(self._instruction.filaments_coordinates)
-
         for filament_coordinates in self._instruction.filaments_coordinates:
             (
                 initial_coordinate,
                 final_coordinate,
             ) = self._find_initial_and_final_filament_coordinates(filament_coordinates)
-
             direction_vector = GeometryMath.direction_vector(
                 initial_coordinate, final_coordinate
             )
-
             printing_speed = filament_coordinates[1][4]
-
             volume = filament_coordinates[2]  # Now this is volume, not E
-
             number_printed_filaments += 1
-
             if final_coordinate[2] > initial_z_coordinate:
                 number_printed_layers += 1
                 initial_z_coordinate = final_coordinate[2]
-
             filament_length = GeometryMath.distance(
                 initial_coordinate, final_coordinate
             )
-
             number_simulation_steps, step_size = self._find_simulation_step_info(
                 filament_length
             )
-
             volumes = Volume.get_volumes_for_filament(
                 number_simulation_steps=number_simulation_steps,
                 total_volume=volume,
@@ -138,7 +132,6 @@ def print(self):
                 printing_speed=printing_speed,
                 printer=self._printer,
             )
-
             self._deposit_filament(
                 number_simulation_steps=number_simulation_steps,
                 step_size=step_size,
@@ -147,6 +140,55 @@ def print(self):
                 volumes=volumes,
             )
 
+    def print_preview_mode(self):
+        """
+        Preview mode: Fast, lightweight visualization. Traces G-code path and fills a capsule (cylinder with rounded endcaps) between last and current point if extruder is on. No physics, no overlap checks. Intended for quick feedback before running the full simulation.
+        """
+        nozzle_radius = self._printer.nozzle_diameter / 2.0
+        voxel_size = self._simulation.voxel_size
+        for filament_coordinates in self._instruction.filaments_coordinates:
+            initial, final = self._find_initial_and_final_filament_coordinates(filament_coordinates)
+            extrusion_start = filament_coordinates[0][3]
+            extrusion_end = filament_coordinates[1][3]
+            # Only deposit if extrusion is positive (extruder is on)
+            if extrusion_end > extrusion_start:
+                # Get bounding box for the capsule
+                min_x = min(initial[0], final[0]) - nozzle_radius
+                max_x = max(initial[0], final[0]) + nozzle_radius
+                min_y = min(initial[1], final[1]) - nozzle_radius
+                max_y = max(initial[1], final[1]) + nozzle_radius
+                min_z = min(initial[2], final[2]) - nozzle_radius
+                max_z = max(initial[2], final[2]) + nozzle_radius
+                i_min = max(0, int(min_x / voxel_size))
+                i_max = min(self.space.shape[0] - 1, int(max_x / voxel_size))
+                j_min = max(0, int(min_y / voxel_size))
+                j_max = min(self.space.shape[1] - 1, int(max_y / voxel_size))
+                k_min = max(0, int(min_z / voxel_size))
+                k_max = min(self.space.shape[2] - 1, int(max_z / voxel_size))
+                # Capsule math: for each voxel, check if its center is within the capsule
+                p1 = np.array(initial)
+                p2 = np.array(final)
+                seg = p2 - p1
+                seg_len2 = np.dot(seg, seg)
+                for i in range(i_min, i_max + 1):
+                    for j in range(j_min, j_max + 1):
+                        for k in range(k_min, k_max + 1):
+                            voxel_center = np.array([
+                                (i + 0.5) * voxel_size,
+                                (j + 0.5) * voxel_size,
+                                (k + 0.5) * voxel_size
+                            ])
+                            # Project voxel_center onto segment
+                            if seg_len2 > 0:
+                                t = np.dot(voxel_center - p1, seg) / seg_len2
+                                t = max(0, min(1, t))
+                                closest = p1 + t * seg
+                            else:
+                                closest = p1
+                            dist2 = np.sum((voxel_center - closest) ** 2)
+                            if dist2 <= nozzle_radius ** 2:
+                                self.space[i, j, k] = 1
+
     def _find_initial_and_final_filament_coordinates(self, filament_coordinates):
         initial = [
             filament_coordinates[0][0] + self.filament_translations["x"],
diff --git a/volco.py b/volco.py
index 8740cae..55361b7 100644
--- a/volco.py
+++ b/volco.py
@@ -114,8 +114,20 @@ def run_simulation(
 if __name__ == "__main__":
     options = Arguments.get_options()
 
+    # Load simulation config from file if provided
+    sim_config = None
+    if options.sim:
+        with open(options.sim, "r") as f:
+            sim_config = json.load(f)
+    # If --preview is set, override preview_mode in config
+    if options.preview:
+        if sim_config is None:
+            sim_config = {}
+        sim_config["preview_mode"] = True
+
     run_simulation(
         gcode_path=options.gcode,
         printer_config_path=options.printer,
-        sim_config_path=options.sim
+        sim_config=sim_config if sim_config is not None else None,
+        sim_config_path=None if sim_config is not None else options.sim
     )

From fed8a7d0674e6c9fcdfe204f0dec88a0ce3d27b8 Mon Sep 17 00:00:00 2001
From: Kyle Grover <kyle@kylegrover.com>
Date: Mon, 17 Nov 2025 19:26:32 -0800
Subject: [PATCH 09/10] improve preview mode voxel processing speed

---
 app/geometry/voxel_space.py | 49 +++++++++++++++++--------------------
 1 file changed, 23 insertions(+), 26 deletions(-)

diff --git a/app/geometry/voxel_space.py b/app/geometry/voxel_space.py
index f0d37bb..5949171 100644
--- a/app/geometry/voxel_space.py
+++ b/app/geometry/voxel_space.py
@@ -146,13 +146,12 @@ def print_preview_mode(self):
         """
         nozzle_radius = self._printer.nozzle_diameter / 2.0
         voxel_size = self._simulation.voxel_size
+        # Use the same logic as the default mode: fill capsule for every filament segment
         for filament_coordinates in self._instruction.filaments_coordinates:
             initial, final = self._find_initial_and_final_filament_coordinates(filament_coordinates)
-            extrusion_start = filament_coordinates[0][3]
-            extrusion_end = filament_coordinates[1][3]
-            # Only deposit if extrusion is positive (extruder is on)
-            if extrusion_end > extrusion_start:
-                # Get bounding box for the capsule
+            # Only fill if the segment has positive volume
+            volume = filament_coordinates[2]
+            if volume > 0:
                 min_x = min(initial[0], final[0]) - nozzle_radius
                 max_x = max(initial[0], final[0]) + nozzle_radius
                 min_y = min(initial[1], final[1]) - nozzle_radius
@@ -165,29 +164,27 @@ def print_preview_mode(self):
                 j_max = min(self.space.shape[1] - 1, int(max_y / voxel_size))
                 k_min = max(0, int(min_z / voxel_size))
                 k_max = min(self.space.shape[2] - 1, int(max_z / voxel_size))
-                # Capsule math: for each voxel, check if its center is within the capsule
+                # Vectorized capsule SDF fill
                 p1 = np.array(initial)
                 p2 = np.array(final)
-                seg = p2 - p1
-                seg_len2 = np.dot(seg, seg)
-                for i in range(i_min, i_max + 1):
-                    for j in range(j_min, j_max + 1):
-                        for k in range(k_min, k_max + 1):
-                            voxel_center = np.array([
-                                (i + 0.5) * voxel_size,
-                                (j + 0.5) * voxel_size,
-                                (k + 0.5) * voxel_size
-                            ])
-                            # Project voxel_center onto segment
-                            if seg_len2 > 0:
-                                t = np.dot(voxel_center - p1, seg) / seg_len2
-                                t = max(0, min(1, t))
-                                closest = p1 + t * seg
-                            else:
-                                closest = p1
-                            dist2 = np.sum((voxel_center - closest) ** 2)
-                            if dist2 <= nozzle_radius ** 2:
-                                self.space[i, j, k] = 1
+                ba = p2 - p1
+                ba_dot_ba = np.dot(ba, ba) if np.dot(ba, ba) > 0 else 1e-8
+                ii, jj, kk = np.meshgrid(
+                    np.arange(i_min, i_max + 1),
+                    np.arange(j_min, j_max + 1),
+                    np.arange(k_min, k_max + 1),
+                    indexing='ij'
+                )
+                voxel_centers = np.stack([
+                    (ii + 0.5) * voxel_size,
+                    (jj + 0.5) * voxel_size,
+                    (kk + 0.5) * voxel_size
+                ], axis=-1)
+                pa = voxel_centers - p1
+                h = np.clip(np.sum(pa * ba, axis=-1) / ba_dot_ba, 0.0, 1.0)
+                capsule_sdf = np.sqrt(np.sum((pa - ba * h[..., None]) ** 2, axis=-1)) - nozzle_radius
+                inside = capsule_sdf <= 0
+                self.space[ii[inside], jj[inside], kk[inside]] = 1
 
     def _find_initial_and_final_filament_coordinates(self, filament_coordinates):
         initial = [

From 54cb0185b98d058a00488adc66bca886ac3e948a Mon Sep 17 00:00:00 2001
From: Kyle Grover <kyle@kylegrover.com>
Date: Mon, 17 Nov 2025 19:51:40 -0800
Subject: [PATCH 10/10] refactor create_box_representation for improved filled
 voxel processing

---
 app/reporter/mesh.py | 129 ++++++++++++++++---------------------------
 1 file changed, 49 insertions(+), 80 deletions(-)

diff --git a/app/reporter/mesh.py b/app/reporter/mesh.py
index 95e74db..d1eacbf 100644
--- a/app/reporter/mesh.py
+++ b/app/reporter/mesh.py
@@ -30,101 +30,69 @@ def create_box_representation(voxel_space, voxel_size):
     trimesh.Trimesh
         A mesh containing only visible faces of filled voxels
     """
-    # Find the indices of filled voxels
-    filled_voxels = np.where(voxel_space > 0)
-    
-    # If no filled voxels, return an empty scene
-    if len(filled_voxels[0]) == 0:
+    # Find filled voxels
+    filled = (voxel_space > 0)
+    if not np.any(filled):
         return trimesh.Scene()
-    
-    # Get the dimensions of the voxel space
+
     max_i, max_j, max_k = voxel_space.shape
-    
-    # Define a unit cube vertices (8 corners)
     unit_cube_vertices = np.array([
-        [-0.5, -0.5, -0.5],  # 0: bottom, back, left
-        [0.5, -0.5, -0.5],   # 1: bottom, back, right
-        [0.5, 0.5, -0.5],    # 2: bottom, front, right
-        [-0.5, 0.5, -0.5],   # 3: bottom, front, left
-        [-0.5, -0.5, 0.5],   # 4: top, back, left
-        [0.5, -0.5, 0.5],    # 5: top, back, right
-        [0.5, 0.5, 0.5],     # 6: top, front, right
-        [-0.5, 0.5, 0.5]     # 7: top, front, left
+        [-0.5, -0.5, -0.5], [0.5, -0.5, -0.5], [0.5, 0.5, -0.5], [-0.5, 0.5, -0.5],
+        [-0.5, -0.5, 0.5], [0.5, -0.5, 0.5], [0.5, 0.5, 0.5], [-0.5, 0.5, 0.5]
     ])
-    
-    # Define the faces for each side of the cube (2 triangles per face)
-    # Each face is associated with a direction (negative or positive x, y, z)
     face_definitions = [
-        # Face indices, Direction to check (di, dj, dk)
-        ([[0, 2, 1], [0, 3, 2]], (0, 0, -1)),  # bottom face (negative z)
-        ([[4, 5, 6], [4, 6, 7]], (0, 0, 1)),   # top face (positive z)
-        ([[0, 1, 5], [0, 5, 4]], (0, -1, 0)),  # back face (negative y)
-        ([[2, 3, 7], [2, 7, 6]], (0, 1, 0)),   # front face (positive y)
-        ([[0, 4, 7], [0, 7, 3]], (-1, 0, 0)),  # left face (negative x)
-        ([[1, 2, 6], [1, 6, 5]], (1, 0, 0))    # right face (positive x)
+        ([[0, 2, 1], [0, 3, 2]], (0, 0, -1)),
+        ([[4, 5, 6], [4, 6, 7]], (0, 0, 1)),
+        ([[0, 1, 5], [0, 5, 4]], (0, -1, 0)),
+        ([[2, 3, 7], [2, 7, 6]], (0, 1, 0)),
+        ([[0, 4, 7], [0, 7, 3]], (-1, 0, 0)),
+        ([[1, 2, 6], [1, 6, 5]], (1, 0, 0))
     ]
-    
-    # First pass: determine which voxels need to have vertices added
-    # and which faces need to be rendered
-    voxels_to_render = {}  # Maps voxel index to list of faces to render
-    
-    for idx, (i, j, k) in enumerate(zip(*filled_voxels)):
-        faces_to_render = []
-        
-        # Check each face to see if it should be rendered
-        for face_idx, (face_triangles, (di, dj, dk)) in enumerate(face_definitions):
-            # Check if this face is at the boundary or adjacent to an empty voxel
-            ni, nj, nk = i + di, j + dj, k + dk
-            
-            # If the neighbor is outside the voxel space or is empty, render this face
-            if (ni < 0 or ni >= max_i or
-                nj < 0 or nj >= max_j or
-                nk < 0 or nk >= max_k or
-                voxel_space[ni, nj, nk] == 0):
-                
-                faces_to_render.append(face_idx)
-        
-        # If this voxel has faces to render, add it to the dictionary
-        if faces_to_render:
-            voxels_to_render[(i, j, k)] = faces_to_render
-    
-    # Second pass: create vertices and faces
+
+    # Vectorized neighbor check for all filled voxels
+    filled_indices = np.array(np.where(filled)).T
     all_vertices = []
     all_faces = []
-    voxel_to_vertex_idx = {}  # Maps voxel coordinates to starting vertex index
-    
-    for idx, (i, j, k) in enumerate(voxels_to_render.keys()):
-        # Calculate the center position of the voxel
-        center = np.array([
-            (i) * voxel_size,
-            (j) * voxel_size,
-            (k) * voxel_size
-        ])
-        
-        # Scale and translate the unit cube vertices for this voxel
-        voxel_vertices = unit_cube_vertices * voxel_size + center
-        
-        # Add vertices for this voxel
-        vertex_start = len(all_vertices)
-        all_vertices.extend(voxel_vertices)
-        voxel_to_vertex_idx[(i, j, k)] = vertex_start
-        
-        # Add faces for this voxel
-        for face_idx in voxels_to_render[(i, j, k)]:
-            face_triangles = face_definitions[face_idx][0]
+    vertex_count = 0
+    for face_idx, (face_triangles, (di, dj, dk)) in enumerate(face_definitions):
+        # Use numpy slicing for neighbor checks
+        if di != 0:
+            if di > 0:
+                mask = np.zeros_like(filled)
+                mask[:-1, :, :] = filled[:-1, :, :] & (~filled[1:, :, :])
+            else:
+                mask = np.zeros_like(filled)
+                mask[1:, :, :] = filled[1:, :, :] & (~filled[:-1, :, :])
+        elif dj != 0:
+            if dj > 0:
+                mask = np.zeros_like(filled)
+                mask[:, :-1, :] = filled[:, :-1, :] & (~filled[:, 1:, :])
+            else:
+                mask = np.zeros_like(filled)
+                mask[:, 1:, :] = filled[:, 1:, :] & (~filled[:, :-1, :])
+        elif dk != 0:
+            if dk > 0:
+                mask = np.zeros_like(filled)
+                mask[:, :, :-1] = filled[:, :, :-1] & (~filled[:, :, 1:])
+            else:
+                mask = np.zeros_like(filled)
+                mask[:, :, 1:] = filled[:, :, 1:] & (~filled[:, :, :-1])
+        face_voxels = np.array(np.where(mask)).T
+        for idx in face_voxels:
+            i, j, k = idx
+            center = np.array([i, j, k]) * voxel_size
+            voxel_vertices = unit_cube_vertices * voxel_size + center
+            v_start = vertex_count
+            all_vertices.extend(voxel_vertices)
+            vertex_count += 8
             for triangle in face_triangles:
-                all_faces.append([t + vertex_start for t in triangle])
-    
-    # Convert lists to numpy arrays
+                all_faces.append([v_start + t for t in triangle])
     if all_vertices and all_faces:
         vertices_array = np.array(all_vertices)
         faces_array = np.array(all_faces)
-        
-        # Create a mesh from the vertices and faces
         mesh = trimesh.Trimesh(vertices=vertices_array, faces=faces_array)
         return mesh
     else:
-        # If no faces were added, return an empty scene
         return trimesh.Scene()
 
 def generate_mesh_from_voxels(voxel_space, voxel_size):
@@ -212,6 +180,7 @@ def export_mesh_to_stl(mesh, file_path, ascii_format=True):
         
         # Use the generic export function which works for both Scene and Trimesh objects
         trimesh.exchange.export.export_mesh(mesh, file_path, **export_options)
+        # mesh.export(file_path, **export_options)
             
         logger.info(f"[Mesh]: STL exported in {format_type} format!")
         return file_path