Merge pull request #10 from swiftcoder/0.1-alpha.0

Extended Marching Cubes and Dual Contouring

Author: swiftcoder

Cargo.toml

@@ -6,7 +6,7 @@ license = "Apache-2.0"
 name = "isosurface"
 readme = "README.md"
 repository = "https://github.com/swiftcoder/isosurface"
-version = "0.0.4"
+version = "0.1.0-alpha.0"
 
 [dev-dependencies]
 cgmath = "^0.17"
@@ -15,12 +15,12 @@ glium = "^0.26"
 glium_text_rusttype = "^0.3"
 
 [profile.dev]
-opt-level = 2
+opt-level = 3
 
 [profile.release]
-opt-level = 3
 lto = true
+opt-level = 3
 
 [[bench]]
-name = "isosurface"
 harness = false
+name = "isosurface"

README.md

@@ -1,22 +1,31 @@
+[![Crates.io](https://img.shields.io/crates/v/isosurface.svg)](https://crates.io/crates/isosurface)
+[![Docs.rs](https://docs.rs/isosurface/badge.svg)](https://docs.rs/isosurface)
+[![Github license](https://img.shields.io/github/license/swiftcoder/isosurface.svg)](https://github.com/swiftcoder/isosurface/blob/trunk/LICENSE.md)
+
 # Isosurface
-Isosurface extraction algorithms for Rust. Currently only a few techniques are implemented.
+Isosurface extraction algorithms implemented in Rust. The classical Marching Cubes and Dual Contouring techniques are included, along with more modern variations on the theme.
+
+In the interest of education, the documentation for each extraction algorithm links to the relevant academic papers.
+
+## Example programs
+`cargo run --example sampler` will execute the sampler, which allows you to compare a variety of algorithms and implicit surfaces.
 
-This crate intentionally has no dependencies to keep the footprint of the library small. The examples rely on the `glium`, `glium_text_rusttype`, and `cgmath` crates.
+ `cargo run --example deferred_rasterisation` will execute a demonstration of GPU-side deferred rasterisation from point clouds. This is a technique pioneered by Gavan Woolery, of [Voxel Quest](https://www.voxelquest.com) fame.
 
-# Marching Cubes
-The Marching Cubes implementation produces perfectly indexed meshes with few duplicate vertices, through the use of a (fairly involved) index caching system. The complexity of the cache could no doubt be reduced through some clever arithmetic, but it is not currently a bottleneck.
+## Dependencies
+This library intentionally has no dependencies. While that requires some redevelopment of common code (i.e. the Vec3 type), it keeps the footprint of the library small, and compile times low for consuming crates. The examples do however rely on the `glium`, `glium_text_rusttype`, and `cgmath` crates, to avoid reinventing the world.
 
-The implementation has been optimised for performance, with memory use kept as a low as possible considering. For an NxNxN voxel chunk, it will allocate roughly NxN of f32 storage for isosurface values, and Nx(N+1) of u32 storage for the index cache.
+## 32-bit indices
+For simplicity vertex indices have been fixed at 32-bits, because for chunks of 32x32x32 and larger you'll often end up with more than 65k vertices. If you are targeting a mobile platform that supports only 16-bit indices, you'll need to keep to smaller chunk sizes, or split the mesh on the output side.
 
-Indices are 32-bit because for chunks of 32x32 and larger you'll typically end up with more than 65k vertices. If you are targeting a mobile platform that supports only 16-bit indices, you'll need to use smaller chunk sizes, and truncate on the output side.
+## Why are optimisations enabled in debug builds?
+Without optimisations enabled, debug builds are around 70x slower. The implementation relies on a lot of nested for loops over integer ranges, and the range iterators themselves entirely dominate the CPU profiles in unoptimised builds.
 
-# Linear Hashed Marching Cubes
-A very efficient algorithm using interleaved integer coordinates to represent octree cells, and storing them in a hash table. Results in better mesh quality than regular marching cubes, and is significantly faster. Memory usage is less predictable, but shouldn't be significantly higher than standard marching cubes.
- 
-# Point Clouds and Deferred Rasterisation
-Point cloud extraction is typically not all that useful, given that point clouds don't contain any data about the actual surface. However, Gavan Woolery (gavanw@) posted an interesting image of reconstructing surface data in image space on the GPU, so I've added a simple example of that.  
+While this can be worked around by converting the `for 0..8` style of loop to a while loop with manual counter, the result is quite unpleasant, and distinctly not in the spirit of rust. I'd rather leave optimisations enabled, and wait for the compiler to become better at handling iterators in debug builds.
 
-# Why are optimisations enabled in debug builds?
-Without optimisations enabled, debug builds are 70x slower (1 minute to extract a 256^3 volume, versus ~800 milliseconds). 
+If you take a dependency on this crate and run into the same issue, you can tell Cargo to compile just this one crate in release mode, by adding the following to your `Cargo.toml`:
 
-The marching cubes implementation relies on a lot of nested for loops over integer ranges, and the range iterators themselves entirely dominate the CPU profiles in unoptimised builds. While this could likely be worked around by converting the `for 0..8` style of loop to a while loop with manual counter, that seems ugly and distinctly not in the spirit of rust. I'd rather leave optimisations enabled, and wait for the compiler to become better at handling iterators.
+```
+[profile.dev.package.isosurface]
+opt-level = 3
+```

benches/isosurface.rs

@@ -1,41 +1,45 @@
+// Copyright 2021 Tristam MacDonald
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
 use criterion::{criterion_group, criterion_main, Criterion};
 use isosurface::{
-    linear_hashed_marching_cubes::LinearHashedMarchingCubes, marching_cubes::MarchingCubes,
-    source::Source,
+    distance::Signed, extractor::IndexedVertices, implicit::Torus, sampler::Sampler,
+    LinearHashedMarchingCubes, MarchingCubes,
 };
 
-fn torus(x: f32, y: f32, z: f32) -> f32 {
-    const R1: f32 = 1.0 / 4.0;
-    const R2: f32 = 1.0 / 10.0;
-    let q_x = ((x * x + y * y).sqrt()).abs() - R1;
-    let len = (q_x * q_x + z * z).sqrt();
-    len - R2
-}
-
-pub struct Torus {}
-
-impl Source for Torus {
-    fn sample(&self, x: f32, y: f32, z: f32) -> f32 {
-        torus(x - 0.5, y - 0.5, z - 0.5)
-    }
-}
-
 fn marching_cubes() {
-    let torus = Torus {};
+    let torus = Torus::new(0.25, 0.1);
+    let sampler = Sampler::new(&torus);
+
     let mut vertices = vec![];
     let mut indices = vec![];
+    let mut extractor = IndexedVertices::new(&mut vertices, &mut indices);
 
-    let mut marching_cubes = MarchingCubes::new(256);
-    marching_cubes.extract(&torus, &mut vertices, &mut indices);
+    let mut marching_cubes = MarchingCubes::<Signed>::new(128);
+    marching_cubes.extract(&sampler, &mut extractor);
 }
 
 fn linear_hashed_marching_cubes() {
-    let torus = Torus {};
+    let torus = Torus::new(0.25, 0.1);
+    let sampler = Sampler::new(&torus);
+
     let mut vertices = vec![];
     let mut indices = vec![];
+    let mut extractor = IndexedVertices::new(&mut vertices, &mut indices);
 
-    let mut marching_cubes = LinearHashedMarchingCubes::new(8);
-    marching_cubes.extract(&torus, &mut vertices, &mut indices);
+    let mut marching_cubes = LinearHashedMarchingCubes::new(7);
+    marching_cubes.extract(&sampler, &mut extractor);
 }
 
 fn marching_cubes_benchmark(c: &mut Criterion) {

examples/common/mod.rs

@@ -1,4 +1,4 @@
-// Copyright 2018 Tristam MacDonald
+// Copyright 2021 Tristam MacDonald
 //
 // Licensed under the Apache License, Version 2.0 (the "License");
 // you may not use this file except in compliance with the License.
@@ -18,10 +18,10 @@ pub mod text;
 use std::mem;
 use std::slice;
 
-/// This is used to reinterpret slices of floats as slices of repr(C) structs, without any
-/// copying. It is optimal, but it is also punching holes in the type system. I hope that Rust
-/// provides safe functionality to handle this in the future. In the meantime, reproduce
-/// this workaround at your own risk.
+/// This is used to reinterpret slices of floats as slices of repr(C) structs,
+/// without any copying. It is optimal, but it is also punching holes in the
+/// type system. I hope that Rust provides safe functionality to handle this in
+/// the future. In the meantime, reproduce this workaround at your own risk.
 pub fn reinterpret_cast_slice<S, T>(input: &[S]) -> &[T] {
     let length_in_bytes = input.len() * mem::size_of::<S>();
     let desired_length = length_in_bytes / mem::size_of::<T>();

examples/common/sources.rs

@@ -1,4 +1,4 @@
-// Copyright 2018 Tristam MacDonald
+// Copyright 2021 Tristam MacDonald
 //
 // Licensed under the Apache License, Version 2.0 (the "License");
 // you may not use this file except in compliance with the License.
@@ -13,68 +13,45 @@
 // limitations under the License.
 
 //! Isosurface definitions for use in multiple examples
+use isosurface::{
+    distance::{Directed, Signed},
+    math::Vec3,
+    source::{HermiteSource, ScalarSource, VectorSource},
+};
 
-use isosurface::source::Source;
+pub trait AllSources: ScalarSource + VectorSource + HermiteSource {}
 
-/// The distance-field equation for a torus
-fn torus(x: f32, y: f32, z: f32) -> f32 {
-    const R1: f32 = 1.0 / 4.0;
-    const R2: f32 = 1.0 / 10.0;
-    let q_x = ((x * x + y * y).sqrt()).abs() - R1;
-    let len = (q_x * q_x + z * z).sqrt();
-    len - R2
-}
-
-pub struct Torus {}
-
-impl Source for Torus {
-    fn sample(&self, x: f32, y: f32, z: f32) -> f32 {
-        torus(x - 0.5, y - 0.5, z - 0.5)
-    }
-}
+impl<S: ScalarSource + VectorSource + HermiteSource> AllSources for S {}
 
-fn abs(x: f32, y: f32, z: f32) -> (f32, f32, f32) {
-    (
-        if x > 0.0 { x } else { -x },
-        if y > 0.0 { y } else { -y },
-        if z > 0.0 { z } else { -z },
-    )
+pub struct DemoSource<'a> {
+    pub source: Box<dyn 'a + AllSources>,
 }
 
-fn max(px: f32, py: f32, pz: f32, qx: f32, qy: f32, qz: f32) -> (f32, f32, f32) {
-    (
-        if px > qx { px } else { qx },
-        if py > qy { py } else { qy },
-        if pz > qz { pz } else { qz },
-    )
-}
-
-/// The distance field equation for a cube
-fn cube(px: f32, py: f32, pz: f32, bx: f32, by: f32, bz: f32) -> f32 {
-    let (ax, ay, az) = abs(px, py, pz);
-    let (dx, dy, dz) = (ax - bx, ay - by, az - bz);
-    let (mx, my, mz) = max(dx, dy, dz, 0.0, 0.0, 0.0);
-    let l = (mx * mx + my * my + mz * mz).sqrt();
-    dx.max(dz.max(dy)).min(0.0) + l
+impl<'a> DemoSource<'a> {
+    pub fn new<S: 'a + AllSources>(source: S) -> Self {
+        Self {
+            source: Box::new(source),
+        }
+    }
 }
 
-/// The distance field equation for a sphere
-fn sphere(x: f32, y: f32, z: f32, r: f32) -> f32 {
-    (x * x + y * y + z * z).sqrt() - r
+impl<'a> ScalarSource for DemoSource<'a> {
+    fn sample_scalar(&self, p: Vec3) -> Signed {
+        let q = p - Vec3::from_scalar(0.5);
+        self.source.sample_scalar(q)
+    }
 }
 
-/// Subtract one distance field from another (i.e. CSG difference operation)
-fn subtract(d1: f32, d2: f32) -> f32 {
-    d2.max(-d1)
+impl<'a> VectorSource for DemoSource<'a> {
+    fn sample_vector(&self, p: Vec3) -> Directed {
+        let q = p - Vec3::from_scalar(0.5);
+        self.source.sample_vector(q)
+    }
 }
 
-pub struct CubeSphere {}
-
-impl Source for CubeSphere {
-    fn sample(&self, x: f32, y: f32, z: f32) -> f32 {
-        subtract(
-            sphere(x - 0.5, y - 0.5, z - 0.5, 0.25),
-            cube(x - 0.5, y - 0.5, z - 0.5, 0.2, 0.2, 0.2),
-        )
+impl<'a> HermiteSource for DemoSource<'a> {
+    fn sample_normal(&self, p: Vec3) -> Vec3 {
+        let q = p - Vec3::from_scalar(0.5);
+        self.source.sample_normal(q)
     }
 }

examples/common/text.rs

@@ -1,4 +1,4 @@
-// Copyright 2018 Tristam MacDonald
+// Copyright 2021 Tristam MacDonald
 //
 // Licensed under the Apache License, Version 2.0 (the "License");
 // you may not use this file except in compliance with the License.
@@ -13,11 +13,11 @@
 // limitations under the License.
 
 //! Conveniences to make the `glium_text` API easier to use in samples.
-
 use cgmath::{Matrix4, Vector3};
 
-/// Produce a transform matrix that will display text at offset column `x`, row `y`, in a
-/// display-filling coordinate space N characters wide and N*aspect rows high.
+/// Produce a transform matrix that will display text at offset column `x`, row
+/// `y`, in a display-filling coordinate space N characters wide and N*aspect
+/// rows high.
 pub fn layout_text(characters_per_row: f32, aspect: f32, x: f32, y: f32) -> Matrix4<f32> {
     let inv_scale = 2.0 / characters_per_row;
     Matrix4::from_translation(Vector3::new(-1.0, -1.0, 0.0))

examples/deferred_rasterisation.rs

@@ -1,4 +1,4 @@
-// Copyright 2018 Tristam MacDonald
+// Copyright 2021 Tristam MacDonald
 //
 // Licensed under the Apache License, Version 2.0 (the "License");
 // you may not use this file except in compliance with the License.
@@ -19,21 +19,22 @@ extern crate isosurface;
 
 mod common;
 
-use crate::common::reinterpret_cast_slice;
-use crate::common::sources::Torus;
+use crate::common::{reinterpret_cast_slice, sources::DemoSource};
 use cgmath::{vec3, Matrix4, Point3, SquareMatrix};
-use glium::glutin;
-use glium::glutin::{
-    dpi::LogicalSize,
-    event::{ElementState, KeyboardInput, VirtualKeyCode, WindowEvent},
-    Api, GlProfile, GlRequest,
+use glium::{
+    glutin::{
+        self,
+        dpi::LogicalSize,
+        event::{ElementState, KeyboardInput, VirtualKeyCode, WindowEvent},
+        Api, GlProfile, GlRequest,
+    },
+    texture::{DepthFormat, DepthTexture2d, MipmapsOption, Texture2d, UncompressedFloatFormat},
+    Surface,
 };
-use glium::texture::{
-    DepthFormat, DepthTexture2d, MipmapsOption, Texture2d, UncompressedFloatFormat,
+use isosurface::{
+    distance::Signed, extractor::OnlyInterleavedNormals, implicit::Torus, sampler::Sampler,
+    source::CentralDifference, PointCloud,
 };
-use glium::Surface;
-use isosurface::point_cloud::PointCloud;
-use isosurface::source::CentralDifference;
 
 #[derive(Copy, Clone)]
 #[repr(C)]
@@ -52,9 +53,9 @@ struct VertexWithNormal {
 
 implement_vertex!(VertexWithNormal, position, normal);
 
-// This technique is derived from an image tweeted by Gavan Woolery (gavanw@). it needs some
-// refinement, but I think I've captured an approximation of his rendering technique.
-// https://twitter.com/gavanw/status/717265068086308865
+// This technique is derived from an image tweeted by Gavan Woolery (gavanw@).
+// it needs some refinement, but I think I've captured an approximation of his
+// rendering technique. https://twitter.com/gavanw/status/717265068086308865
 
 fn main() {
     let events_loop = glutin::event_loop::EventLoop::new();
@@ -76,13 +77,15 @@ fn main() {
 
     let subdivisions = 64;
 
-    let torus = Torus {};
+    let torus = DemoSource::new(Torus::new(0.25, 0.1));
     let central_difference = CentralDifference::new(torus);
+    let sampler = Sampler::new(&central_difference);
 
     let mut vertices = vec![];
-    let mut marcher = PointCloud::new(subdivisions);
+    let mut extractor = OnlyInterleavedNormals::new(&mut vertices, &sampler);
+    let mut marcher = PointCloud::<Signed>::new(subdivisions);
 
-    marcher.extract_midpoints_with_normals(&central_difference, &mut vertices);
+    marcher.extract(&sampler, &mut extractor);
 
     let vertex_buffer: glium::VertexBuffer<VertexWithNormal> = {
         glium::VertexBuffer::new(&display, reinterpret_cast_slice(&vertices))
@@ -226,7 +229,8 @@ fn main() {
     );
     let model = Matrix4::identity();
 
-    // We need two textures to ping-pong between, and one of them needs an attached depth buffer for the initial pass
+    // We need two textures to ping-pong between, and one of them needs an attached
+    // depth buffer for the initial pass
     let position1 = Texture2d::empty_with_format(
         &display,
         UncompressedFloatFormat::F32F32F32F32,
@@ -358,8 +362,8 @@ fn main() {
                 .expect("failed to draw to surface");
         }
 
-        // pass 1 through N-1, ping-pong render both buffers in turn, spreading the points across
-        // the faces of their respective cubes
+        // pass 1 through N-1, ping-pong render both buffers in turn, spreading the
+        // points across the faces of their respective cubes
         for i in 0..3 {
             let framebuffer = if i % 2 == 0 {
                 &mut framebuffer2
@@ -390,7 +394,8 @@ fn main() {
                 .expect("failed to draw to surface");
         }
 
-        // final pass, composite the last buffer to the screen, performing lighting in the process
+        // final pass, composite the last buffer to the screen, performing lighting in
+        // the process
         {
             let mut surface = display.draw();
             surface.clear_color_and_depth((0.306, 0.267, 0.698, 0.0), 1.0);