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 efcf00a..a3e181f 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,24 +37,56 @@ def deposit_sphere( ), ) - return voxel_space + return voxel_space_out + + def fill_voxels(self, voxel_space_obj, radius, lower_indexes, upper_indexes): + # 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 - def fill_voxels(self, voxel_space, radius, lower_indexes, upper_indexes): empty_voxels = GeometryMath.find_empty_voxels_in_space( - voxel_space, lower_indexes, upper_indexes + 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 - - return voxel_space + if not empty_voxels: + 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 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 + + 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) + current_vals = space[xi, yj, zk] + new_mask = current_vals == 0 + n_new = int(np.count_nonzero(new_mask)) + if n_new > 0: + # 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 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) @@ -85,56 +120,62 @@ def find_sphere_limits(self, radius, nozzle_height): """ def deform_voxel_space_for_big_spheres(self, voxel_space, radius): + # 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, max_index, axis_number): - size = voxel_space.shape + def _maybe_expand_voxel_space(self, voxel_space_obj, max_index, axis_number): + # 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: - return voxel_space + return voxel_space_obj number_to_be_added = max_index - index_size + 1 - new_size = list(size) - new_size[axis_number] = 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) - mat_add = np.zeros(new_size, dtype=np.int8) + # Create only the additional block to concatenate + mat_add_size = list(size) + 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) + new_space = np.concatenate((space, mat_add), axis=axis_number) - return voxel_space + if is_voxel_space: + voxel_space_obj.space = new_space + return voxel_space_obj - def _deposit_sphere(self, radius, voxel_space, nozzle_height, target_volume): - copy_voxel_space = voxel_space.copy() + return new_space - copy_voxel_space = self.deform_voxel_space_for_big_spheres( - copy_voxel_space, radius - ) + def _deposit_sphere(self, radius, voxel_space, nozzle_height, target_volume): + # 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) - copy_voxel_space = self.fill_voxels( - copy_voxel_space, radius, lower_indexes, upper_indexes - ) + vs = self.fill_voxels(vs, radius, lower_indexes, upper_indexes) - current_volume = GeometryMath.calculate_filled_volume( - copy_voxel_space, self.voxel_size - ) + current_volume = GeometryMath.calculate_filled_volume(vs, self.voxel_size) volume_overshoot = current_volume / target_volume - 1.0 - return volume_overshoot, copy_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/geometry/voxel_space.py b/app/geometry/voxel_space.py index 819d576..d4c75cf 100644 --- a/app/geometry/voxel_space.py +++ b/app/geometry/voxel_space.py @@ -45,12 +45,54 @@ 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. + # 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, ) + # reset the running counter when allocating space + self._filled_voxels_count = 0 def print(self): number_printed_filaments = 0 @@ -135,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): @@ -183,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/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/app/solvers/bisection_method.py b/app/solvers/bisection_method.py index 953ce58..80b5e09 100644 --- a/app/solvers/bisection_method.py +++ b/app/solvers/bisection_method.py @@ -23,35 +23,43 @@ 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, 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 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") or hasattr(out, "space"): + 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" 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,