Fixed bug where last planes on each axis of SDF were not being used

[chris-marrero]

Say we had an SDF with shape (10, 10, 10). It would have 1000 elements. Previously, surface_nets would assert!((shape.linearize(max) as usize) < sdf.len()). Calling surface_nets with max = (10, 10, 10) would result in shape.linearize((10, 10, 10)) = 10 + (10 * 10) + (10 * 10 * 10) = 1110. You must call it with max = (9, 9, 9) in order to get shape.linearize((9, 9, 9)) = 9 + (9 * 10) + (9 * 10 * 10) = 999. The problem with this is that in estimate_surface, the indices are [min, max) and do not include the max values. So, we end up excluding the max values twice and lose the final planes on each axis of the SDF.

I propose modifying the assert to subtract 1 from each value max in order to match the size of the loop in estimate_surface. I also added an assert to ensure we don't have negative indices for max and bumped the version number. Maybe it should be 0.3 instead of 0.2.1 as it will likely cause changes in results.

I do have one question about the assert on line 138. Surely shape.linearize(min) MUST be less than shape.linearize(max) to produce any result. If so, we could modify it to a < and assert shape.linearize(min) >= 0.

[bonsairobo]

I'm pretty sure the original code is correct, because I've never had issues with missing geometry even when putting together multiple chunks.

I think maybe you are just misinterpreting the parameters. The loop in estimate_surface still samples from max because there are 8 samples per iteration of the loop. So [min, max] is exactly the extent of voxels that will be sampled. The only requirement is that the linear indices don't exceed the length of the SDF buffer, which is why min and max are checked.

I also added an assert to ensure we don't have negative indices for max and bumped the version number.

It is not possible for max to have a negative linear index because it is specified in [u32; 3] coordinates, which always yield a nonnegative linear index.

If you are concerned about missing the back planes, then you need to add more samples to your buffer. Remember that surface nets is a dual method, meaning it only creates one vertex per cell of 8 samples.

[chris-marrero]

Maybe I'm misunderstanding how your implementation works, but here's a minimal example that shows why I made the change I made. I simply create a 3x3 SDF and set the inside node to negative. From what I understand, this should create a cube. However, it will only create 3 faces. If I run it with my PR and call surface_nets(max: [3, 3, 3]), I get all 6 faces.

Also, you're totally right about the the negative indices.

use std::f32::consts::FRAC_PI_4;

use bevy::{
    prelude::*,
    render::{mesh::Indices, render_resource::PrimitiveTopology},
};
use fast_surface_nets::{
    ndshape::{ConstShape, ConstShape3u32},
    surface_nets, SurfaceNetsBuffer,
};

fn setup(mut commands: Commands, mut meshes: ResMut<Assets<Mesh>>) {
    type SdfShape = ConstShape3u32<3, 3, 3>;
    let mut sdf = [1.0; SdfShape::USIZE];
    sdf[SdfShape::linearize([1, 1, 1]) as usize] = -1.0;

    let mut output = SurfaceNetsBuffer::default();
    surface_nets(&sdf, &SdfShape {}, [0, 0, 0], [2, 2, 2], &mut output);
    output.positions.iter_mut().for_each(|p| {
        *p = Quat::from_euler(EulerRot::XYZ, FRAC_PI_4, 0.0, FRAC_PI_4)
            .mul_vec3(Vec3::new(p[0] - 1.0, p[1] - 1.0, p[2] - 1.0))
            .into();
    });

    let mut mesh = Mesh::new(PrimitiveTopology::TriangleList);
    mesh.insert_attribute(Mesh::ATTRIBUTE_POSITION, output.positions);
    mesh.set_indices(Some(Indices::U32(output.indices)));

    commands.spawn(PbrBundle {
        mesh: meshes.add(mesh),
        ..Default::default()
    });

    commands.spawn(Camera3dBundle {
        transform: Transform::from_xyz(2.0, 2.0, 2.0).looking_at(Vec3::ZERO, Vec3::Y),
        ..default()
    });
}

fn update(mut query: Query<(&mut Transform, With<Handle<Mesh>>)>) {
    for (mut t, _) in query.iter_mut() {
        t.rotate_axis(Vec3::new(0.0, 1.0, 0.0), 0.01);
    }
}

fn main() {
    App::new()
        .add_plugins(DefaultPlugins)
        .add_startup_system(setup)
        .add_system(update)
        .run()
}

[bonsairobo]

@chris-marrero You are right! I forgot that there are special cases to skip the maximum planes in a chunk, because I wanted to make this usable with multiple chunks without having neighbors overlap.

fast-surface-nets-rs/src/lib.rs

Line 306 in d1570d4

if y != miny && z != minz && x != maxx - 1 {

If you'd like, we can make this feature optional.

[chris-marrero]

I noticed that if I removed the unsafe code (which I assume is to allow for sparse datasets), my original code would access out of bounds. My initial solution merely changed maxx, maxy, maxz in make_all_quads instead of changing the logic itself. So, my new commit conditionally ignores the final plane check for plane generation. I named the feature eval-max-plane, but I am open to other names.