Skip to content

⚡ Bolt: Optimized Marching Cubes Isosurface Extraction #381

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
fderuiter wants to merge 1 commit into
base: main
Choose a base branch
from
Draft

⚡ Bolt: Optimized Marching Cubes Isosurface Extraction #381

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
68 changes: 41 additions & 27 deletions math_explorer/src/applied/isosurface/marching_cubes.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -143,15 +143,15 @@ pub fn extract_isosurface(grid: &VoxelGrid, threshold: f32) -> Result<Mesh, Isos
// Cache for the "Right Face" gradients of the previous iteration (x-1).
// Corresponds to vertices 1, 2, 5, 6 of (x-1), which become 0, 3, 4, 7 of (x).
let mut cached_gradients: Option<[Point3D; 4]> = None;
// Cache for "Right Face" values of previous iteration.
let mut cached_values: Option<[f32; 4]> = None;

for x in 0..grid.width - 1 {
let base_idx = zy_base + x;
let x_pos = grid.origin.x + (x as f32) * grid.voxel_size.x;

// 1. Determine the index of the case (0-255)
let mut cube_index = 0;
let mut corner_values = [0.0; 8];
let mut corner_pos = [Point3D::new(0.0, 0.0, 0.0); 8];
let mut corner_normals = [Point3D::new(0.0, 0.0, 0.0); 8];

// Direct access for corner values to avoid redundant index calculation
Expand All @@ -165,44 +165,53 @@ pub fn extract_isosurface(grid: &VoxelGrid, threshold: f32) -> Result<Mesh, Isos
// Max = (D-2)*Sz + (H-2)*Sy + (W-2) + 1 + Sy + Sz
// = (D-1)*Sz + (H-1)*Sy + W-1
// Which is exactly the last element. So indices are valid.
let v0 = unsafe { *data.get_unchecked(base_idx) };
let v1 = unsafe { *data.get_unchecked(base_idx + 1) };
let v2 = unsafe { *data.get_unchecked(base_idx + 1 + stride_y) };
let v3 = unsafe { *data.get_unchecked(base_idx + stride_y) };
let v4 = unsafe { *data.get_unchecked(base_idx + stride_z) };
let v5 = unsafe { *data.get_unchecked(base_idx + 1 + stride_z) };
let v6 = unsafe { *data.get_unchecked(base_idx + 1 + stride_y + stride_z) };
let v7 = unsafe { *data.get_unchecked(base_idx + stride_y + stride_z) };

corner_values[0] = v0;

// Profiler Optimization: Reuse values from previous iteration (cache)
// Left Face (0, 3, 4, 7) comes from previous Right Face (1, 2, 5, 6)
let (v0, v3, v4, v7) = if let Some(vals) = cached_values {
(vals[0], vals[1], vals[2], vals[3])
} else {
unsafe {
(
*data.get_unchecked(base_idx),
*data.get_unchecked(base_idx + stride_y),
*data.get_unchecked(base_idx + stride_z),
*data.get_unchecked(base_idx + stride_y + stride_z),
)
}
};

// Right Face (1, 2, 5, 6) - Always new reads
let next_x_idx = base_idx + 1;
let v1 = unsafe { *data.get_unchecked(next_x_idx) };
let v2 = unsafe { *data.get_unchecked(next_x_idx + stride_y) };
let v5 = unsafe { *data.get_unchecked(next_x_idx + stride_z) };
let v6 = unsafe { *data.get_unchecked(next_x_idx + stride_y + stride_z) };

// Update cache for next iteration
cached_values = Some([v1, v2, v5, v6]);

if v0 < threshold {
cube_index |= 1;
}
corner_values[1] = v1;
if v1 < threshold {
cube_index |= 2;
}
corner_values[2] = v2;
if v2 < threshold {
cube_index |= 4;
}
corner_values[3] = v3;
if v3 < threshold {
cube_index |= 8;
}
corner_values[4] = v4;
if v4 < threshold {
cube_index |= 16;
}
corner_values[5] = v5;
if v5 < threshold {
cube_index |= 32;
}
corner_values[6] = v6;
if v6 < threshold {
cube_index |= 64;
}
corner_values[7] = v7;
if v7 < threshold {
cube_index |= 128;
}
Expand All @@ -214,19 +223,24 @@ pub fn extract_isosurface(grid: &VoxelGrid, threshold: f32) -> Result<Mesh, Isos
continue;
}

// Construct full array only when needed
let corner_values = [v0, v1, v2, v3, v4, v5, v6, v7];

// Only compute positions if needed
let next_x_pos = x_pos + grid.voxel_size.x;
let next_y_pos = y_pos + grid.voxel_size.y;
let next_z_pos = z_pos + grid.voxel_size.z;

corner_pos[0] = Point3D::new(x_pos, y_pos, z_pos);
corner_pos[1] = Point3D::new(next_x_pos, y_pos, z_pos);
corner_pos[2] = Point3D::new(next_x_pos, next_y_pos, z_pos);
corner_pos[3] = Point3D::new(x_pos, next_y_pos, z_pos);
corner_pos[4] = Point3D::new(x_pos, y_pos, next_z_pos);
corner_pos[5] = Point3D::new(next_x_pos, y_pos, next_z_pos);
corner_pos[6] = Point3D::new(next_x_pos, next_y_pos, next_z_pos);
corner_pos[7] = Point3D::new(x_pos, next_y_pos, next_z_pos);
let corner_pos = [
Point3D::new(x_pos, y_pos, z_pos),
Point3D::new(next_x_pos, y_pos, z_pos),
Point3D::new(next_x_pos, next_y_pos, z_pos),
Point3D::new(x_pos, next_y_pos, z_pos),
Point3D::new(x_pos, y_pos, next_z_pos),
Point3D::new(next_x_pos, y_pos, next_z_pos),
Point3D::new(next_x_pos, next_y_pos, next_z_pos),
Point3D::new(x_pos, next_y_pos, next_z_pos),
];

// Profiler Optimization: Lazy Gradient Computation & Sliding Window

Expand Down