marching-cubes 0.1.2

use super::container::*;
use super::tables::*;

/// A struct representing the Marching Cubes algorithm.
///
/// The Marching Cubes algorithm is used to extract a triangular mesh representation of a three-dimensional surface from a scalar field.
///
/// # Examples
///
/// ```
/// use iso::{MarchingCubes, GridCell, Triangle};
///
/// let grid = GridCell {
///     positions: [
///         [0.0, 0.0, 0.0],
///         [1.0, 0.0, 0.0],
///         [1.0, 1.0, 0.0],
///         [0.0, 1.0, 0.0],
///         [0.0, 0.0, 1.0],
///         [1.0, 0.0, 1.0],
///         [1.0, 1.0, 1.0],
///         [0.0, 1.0, 1.0]
///     ],
///     value: [0.0, 0.5, 0.5, 1.0, 0.0, 1.0, 1.0, 0.0],
/// };
///
/// let isolevel = 0.5;
/// let mc = MarchingCubes::new(isolevel, grid);
/// let mut triangles = vec![];
/// let triangle_count = mc.polygonise(&mut triangles);
///
/// assert_eq!(triangle_count, 4);
/// ```
#[derive(Debug, Clone, Copy)]
pub struct MarchingCubes {
    iso_value: f32,
    grid: GridCell,
}

#[derive(Debug, Clone, Copy)]
pub struct Triangle {
    pub positions: [Vec3; 3],
}

#[derive(Debug, Clone, Copy)]
pub struct GridCell {
    pub positions: [Vec3; 8],
    pub value: [f32; 8],
}

impl MarchingCubes {
    /// Creates a new instance of the `MarchingCubes` struct.
    ///
    /// # Parameters
    ///
    /// * `iso_value`: A float value representing the isolevel to use when determining which parts of the grid cell are inside and outside the surface.
    /// * `grid`: A `GridCell` struct representing the cell to be polygonised. It contains the positions of the cell's vertices and their values.
    ///
    /// # Returns
    ///
    /// A new instance of the `MarchingCubes` struct.
    pub fn new(iso_value: f32, grid: GridCell) -> Self {
        Self {
            iso_value,
            grid
        }
    }

    /// Polygonises a grid cell based on the stored isolevel, populating the given `triangles` slice with the resulting triangles.
    /// Look at [polygonise] for more information regarding this function.
    ///
    /// # Parameters
    ///
    /// * `triangles`: A mutable slice of `Triangle` structs to be filled with the resulting triangles.
    ///
    /// # Returns
    ///
    /// The number of triangles generated, as an `i32`.
    pub fn polygonise(self, triangles: &mut [Triangle]) -> i32 {
        polygonise(self.grid, self.iso_value, triangles)
    }
}

/// Polygonises a grid cell based on the given isolevel, populating the given
/// `triangles` slice with the resulting triangles.
///
/// # Parameters
///
/// * `grid_cell`: A `GridCell` struct representing the cell to be polygonised. It contains the positions of the cell's vertices and their values.
/// * `isolevel`: A float value representing the isolevel to use when determining
///               which parts of the cell are inside and outside the surface. The isolevel separates the grid cell into parts that are either inside or outside the surface being represented.
/// * `triangles`: A mutable slice of `Triangle` structs to be filled with the
///                resulting triangles. The triangles are generated by triangulating the parts of the grid cell that are inside the surface.
///
/// # Returns
///
/// The number of triangles generated, as an `i32`.
///
/// # Example
///
/// ```
/// use iso::{GridCell, Triangle, polygonise};
///
/// let grid_cell = GridCell {
///     positions: [
///         [0.0, 0.0, 0.0],
///         [1.0, 0.0, 0.0],
///         [1.0, 1.0, 0.0],
///         [0.0, 1.0, 0.0],
///         [0.0, 0.0, 1.0],
///         [1.0, 0.0, 1.0],
///         [1.0, 1.0, 1.0],
///         [0.0, 1.0, 1.0]
///     ],
///     value: [0.0, 0.5, 0.5, 1.0, 0.0, 1.0, 1.0, 0.0],
/// };
/// let mut triangles = vec![];
///
/// let isolevel = 0.5;
/// let result = polygonise(grid_cell, isolevel, &mut triangles);
///
/// assert_eq!(result, 4);
/// ```
pub fn polygonise(grid_cell: GridCell, isolevel: f32, triangles: &mut [Triangle]) -> i32 {
    let mut cube_index: usize;
    let mut vertices_list: [Vec3; 12] = [empty_vec3(); 12];

    cube_index = 0;

    if grid_cell.value[0] < isolevel {
        cube_index |= 1;
    }

    if grid_cell.value[1] < isolevel {
        cube_index |= 2;
    }

    if grid_cell.value[2] < isolevel {
        cube_index |= 4;
    }
    if grid_cell.value[3] < isolevel {
        cube_index |= 8;
    }
    if grid_cell.value[4] < isolevel {
        cube_index |= 16;
    }
    if grid_cell.value[5] < isolevel {
        cube_index |= 32;
    }
    if grid_cell.value[6] < isolevel {
        cube_index |= 64;
    }
    if grid_cell.value[7] < isolevel {
        cube_index |= 128;
    }

    if EDGE_TABLE[cube_index] == 0 {
        return 0i32;
    }

    if EDGE_TABLE[cube_index] & 1 == 0 {
        vertices_list[0] = interpolate_vertex(
            isolevel,
            grid_cell.positions[0],
            grid_cell.positions[1],
            vector2(grid_cell.value[0], grid_cell.value[1]),
        )
    }

    if EDGE_TABLE[cube_index] & 2 == 0 {
        vertices_list[1] = interpolate_vertex(
            isolevel,
            grid_cell.positions[1],
            grid_cell.positions[2],
            vector2(grid_cell.value[1], grid_cell.value[2]),
        )
    }

    if EDGE_TABLE[cube_index] & 4 == 0 {
        vertices_list[2] = interpolate_vertex(
            isolevel,
            grid_cell.positions[2],
            grid_cell.positions[3],
            vector2(grid_cell.value[2], grid_cell.value[3]),
        )
    }

    if EDGE_TABLE[cube_index] & 8 == 0 {
        vertices_list[3] = interpolate_vertex(
            isolevel,
            grid_cell.positions[3],
            grid_cell.positions[0],
            vector2(grid_cell.value[3], grid_cell.value[0]),
        )
    }

    if EDGE_TABLE[cube_index] & 16 == 0 {
        vertices_list[4] = interpolate_vertex(
            isolevel,
            grid_cell.positions[4],
            grid_cell.positions[5],
            vector2(grid_cell.value[4], grid_cell.value[5]),
        )
    }

    if EDGE_TABLE[cube_index] & 32 == 0 {
        vertices_list[5] = interpolate_vertex(
            isolevel,
            grid_cell.positions[5],
            grid_cell.positions[6],
            vector2(grid_cell.value[5], grid_cell.value[6]),
        )
    }

    if EDGE_TABLE[cube_index] & 64 == 0 {
        vertices_list[6] = interpolate_vertex(
            isolevel,
            grid_cell.positions[6],
            grid_cell.positions[7],
            vector2(grid_cell.value[6], grid_cell.value[7]),
        )
    }

    if EDGE_TABLE[cube_index] & 128 == 0 {
        vertices_list[7] = interpolate_vertex(
            isolevel,
            grid_cell.positions[7],
            grid_cell.positions[4],
            vector2(grid_cell.value[7], grid_cell.value[4]),
        )
    }

    if EDGE_TABLE[cube_index] & 256 == 0 {
        vertices_list[8] = interpolate_vertex(
            isolevel,
            grid_cell.positions[0],
            grid_cell.positions[4],
            vector2(grid_cell.value[0], grid_cell.value[4]),
        )
    }
    if EDGE_TABLE[cube_index] & 512 == 0 {
        vertices_list[9] = interpolate_vertex(
            isolevel,
            grid_cell.positions[1],
            grid_cell.positions[5],
            vector2(grid_cell.value[1], grid_cell.value[5]),
        )
    }

    if EDGE_TABLE[cube_index] & 1024 == 0 {
        vertices_list[10] = interpolate_vertex(
            isolevel,
            grid_cell.positions[2],
            grid_cell.positions[6],
            vector2(grid_cell.value[2], grid_cell.value[6]),
        )
    }

    if EDGE_TABLE[cube_index] & 2048 == 0 {
        vertices_list[11] = interpolate_vertex(
            isolevel,
            grid_cell.positions[3],
            grid_cell.positions[7],
            vector2(grid_cell.value[3], grid_cell.value[7]),
        )
    }

    let mut triangle_num = 0;

    for mut i in 0.. {
        let tri = TRI_TABLE[cube_index][i];

        if tri == -1 {
            break;
        }

        i += 3;

        (*triangles)[triangle_num].positions[0] = vertices_list[TRI_TABLE[cube_index][i] as usize];
        (*triangles)[triangle_num].positions[1] =
            vertices_list[TRI_TABLE[cube_index][i + 1] as usize];
        (*triangles)[triangle_num].positions[2] =
            vertices_list[TRI_TABLE[cube_index][i + 2] as usize];

        triangle_num += 1;
    }

    0i32
}

/// Interpolates a vertex using `isolevel`, `point1`, `point2`, and `alpha_points`
///
/// The `isolevel` is used to determine the position of the vertex between `point1` and `point2`.
/// `alpha_points` is used to determine the relative values of `isolevel` and `point1` and `point2`.
/// `point1` and `point2` represent the two points that the vertex will be interpolated between.
///
/// # Parameters
///
/// * `isolevel` - The level of iso to interpolate the vertex at
/// * `point1` - The first point that the vertex will be interpolated between
/// * `point2` - The second point that the vertex will be interpolated between
/// * `alpha_points` - A Vec2 that contains the relative values of `isolevel` and `point1` and `point2`
///
/// # Returns
///
/// A `Vec3` representing the interpolated vertex
///
/// # Example
///
/// ```rust
/// let point1 = vec3(0.0, 0.0, 0.0);
/// let point2 = vec3(1.0, 1.0, 1.0);
/// let alpha_points = vec2(0.5, 0.5);
/// let isolevel = 0.5;
///
/// let result = interpolate_vertex(isolevel, point1, point2, alpha_points);
///
/// assert_eq!(result, vec3(0.5, 0.5, 0.5));
/// ```
fn interpolate_vertex(isolevel: f32, point1: Vec3, point2: Vec3, alpha_points: Vec2) -> Vec3 {
    const ISO_THRESHOLD: f32 = 0.00001;

    let mut point = empty_vec3();
    let factor = (isolevel - alpha_points[0]) / (alpha_points[1] - alpha_points[0]);

    if (isolevel - alpha_points[0]).abs() < ISO_THRESHOLD
        || (alpha_points[0] - alpha_points[1]) < ISO_THRESHOLD
    {
        return point1;
    } else if (isolevel - alpha_points[1]).abs() < ISO_THRESHOLD {
        return point2;
    }

    point[0] = point1[0] + factor * (point2[0] - point1[0]);
    point[1] = point1[1] + factor * (point2[1] - point1[1]);
    point[2] = point1[2] + factor * (point2[2] - point1[2]);
    point
}