marching_cubes.py

# Author: Joelene Hales, 2024

import numpy as np
import glm
import os
import sys


# Scalar fields
def wavy(x, y, z):
    """ Scalar field of a wavy surface, defined by the function 
        f(x, y, z) =  y - sin(x)cos(z) 
    """
    return y - np.sin(x) * np.cos(z)

def hyperboloid(x, y, z):
    """ Scalar field of a hyperboloid, defined by the function
        f(x, y, z) =  x^2 - y^2 - z^2 - z 
    """
    return x**2 - y**2 - z**2 - z


# Helper functions used in the marching cubes algorithm
def _lookup_configuration(case):
    """ Lookup which of the cube's edge midpoints define the vertices of the
    triangle(s) in a case of marching cubes.
    
    Parameters
    ----------
    case : int
        Case of marching cubes encoded as an 8-bit number

    Returns
    -------
    configuration : list[int]
        Edge indices corresponding to the given case of marching cubes
    
    """

    # Define lookup table. Each index of the table corresponds to one case of
    # marching cubes. Each list element corresponds to one edge of the cube.
    lookup_table = [
        [],
        [0, 8, 3],    # Index 1: back bottom left corner of cube
        [0, 1, 9],    # Index 2: back bottom right corner of cube 
        [1, 8, 3, 9, 8, 1],
        [1, 2, 10],   # Index 4: front bottom right corner of cube
        [0, 8, 3, 1, 2, 10],
        [9, 2, 10, 0, 2, 9],
        [2, 8, 3, 2, 10, 8, 10, 9, 8],
        [3, 11, 2],   # Index 8: front bottom left corner of cube
        [0, 11, 2, 8, 11, 0],
        [1, 9, 0, 2, 3, 11],
        [1, 11, 2, 1, 9, 11, 9, 8, 11],
        [3, 10, 1, 11, 10, 3],
        [0, 10, 1, 0, 8, 10, 8, 11, 10],
        [3, 9, 0, 3, 11, 9, 11, 10, 9],
        [9, 8, 10, 10, 8, 11],
        [4, 7, 8],   # Index 16: bottom top left corner of cube
        [4, 3, 0, 7, 3, 4],
        [0, 1, 9, 8, 4, 7],
        [4, 1, 9, 4, 7, 1, 7, 3, 1],
        [1, 2, 10, 8, 4, 7],
        [3, 4, 7, 3, 0, 4, 1, 2, 10],
        [9, 2, 10, 9, 0, 2, 8, 4, 7],
        [2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4],
        [8, 4, 7, 3, 11, 2],
        [11, 4, 7, 11, 2, 4, 2, 0, 4],
        [9, 0, 1, 8, 4, 7, 2, 3, 11],
        [4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1],
        [3, 10, 1, 3, 11, 10, 7, 8, 4],
        [1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4],
        [4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3],
        [4, 7, 11, 4, 11, 9, 9, 11, 10],
        [9, 5, 4],   # Index 32: bottom top right corner of cube
        [9, 5, 4, 0, 8, 3],
        [0, 5, 4, 1, 5, 0],
        [8, 5, 4, 8, 3, 5, 3, 1, 5],
        [1, 2, 10, 9, 5, 4],
        [3, 0, 8, 1, 2, 10, 4, 9, 5],
        [5, 2, 10, 5, 4, 2, 4, 0, 2],
        [2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8],
        [9, 5, 4, 2, 3, 11],
        [0, 11, 2, 0, 8, 11, 4, 9, 5],
        [0, 5, 4, 0, 1, 5, 2, 3, 11],
        [2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5],
        [10, 3, 11, 10, 1, 3, 9, 5, 4],
        [4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10],
        [5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3],
        [5, 4, 8, 5, 8, 10, 10, 8, 11],
        [9, 7, 8, 5, 7, 9],
        [9, 3, 0, 9, 5, 3, 5, 7, 3],
        [0, 7, 8, 0, 1, 7, 1, 5, 7],
        [1, 5, 3, 3, 5, 7],
        [9, 7, 8, 9, 5, 7, 10, 1, 2],
        [10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3],
        [8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2],
        [2, 10, 5, 2, 5, 3, 3, 5, 7],
        [7, 9, 5, 7, 8, 9, 3, 11, 2],
        [9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11],
        [2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7],
        [11, 2, 1, 11, 1, 7, 7, 1, 5],
        [9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11],
        [5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0],
        [11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0],
        [11, 10, 5, 7, 11, 5],
        [10, 6, 5],   # Index 64: front top right corner of cube
        [0, 8, 3, 5, 10, 6],
        [9, 0, 1, 5, 10, 6],
        [1, 8, 3, 1, 9, 8, 5, 10, 6],
        [1, 6, 5, 2, 6, 1],
        [1, 6, 5, 1, 2, 6, 3, 0, 8],
        [9, 6, 5, 9, 0, 6, 0, 2, 6],
        [5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8],
        [2, 3, 11, 10, 6, 5],
        [11, 0, 8, 11, 2, 0, 10, 6, 5],
        [0, 1, 9, 2, 3, 11, 5, 10, 6],
        [5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11],
        [6, 3, 11, 6, 5, 3, 5, 1, 3],
        [0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6],
        [3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9],
        [6, 5, 9, 6, 9, 11, 11, 9, 8],
        [5, 10, 6, 4, 7, 8],
        [4, 3, 0, 4, 7, 3, 6, 5, 10],
        [1, 9, 0, 5, 10, 6, 8, 4, 7],
        [10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4],
        [6, 1, 2, 6, 5, 1, 4, 7, 8],
        [1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7],
        [8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6],
        [7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9],
        [3, 11, 2, 7, 8, 4, 10, 6, 5],
        [5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11],
        [0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6],
        [9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6],
        [8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6],
        [5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11],
        [0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7],
        [6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9],
        [10, 4, 9, 6, 4, 10],
        [4, 10, 6, 4, 9, 10, 0, 8, 3],
        [10, 0, 1, 10, 6, 0, 6, 4, 0],
        [8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10],
        [1, 4, 9, 1, 2, 4, 2, 6, 4],
        [3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4],
        [0, 2, 4, 4, 2, 6],
        [8, 3, 2, 8, 2, 4, 4, 2, 6],
        [10, 4, 9, 10, 6, 4, 11, 2, 3],
        [0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6],
        [3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10],
        [6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1],
        [9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3],
        [8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1],
        [3, 11, 6, 3, 6, 0, 0, 6, 4],
        [6, 4, 8, 11, 6, 8],
        [7, 10, 6, 7, 8, 10, 8, 9, 10],
        [0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10],
        [10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0],
        [10, 6, 7, 10, 7, 1, 1, 7, 3],
        [1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7],
        [2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9],
        [7, 8, 0, 7, 0, 6, 6, 0, 2],
        [7, 3, 2, 6, 7, 2],
        [2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7],
        [2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7],
        [1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11],
        [11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1],
        [8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6],
        [0, 9, 1, 11, 6, 7],
        [7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0],
        [7, 11, 6],
        [7, 6, 11],   # Index 128: front top left corner of cube
        [3, 0, 8, 11, 7, 6],
        [0, 1, 9, 11, 7, 6],
        [8, 1, 9, 8, 3, 1, 11, 7, 6],
        [10, 1, 2, 6, 11, 7],
        [1, 2, 10, 3, 0, 8, 6, 11, 7],
        [2, 9, 0, 2, 10, 9, 6, 11, 7],
        [6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8],
        [7, 2, 3, 6, 2, 7],
        [7, 0, 8, 7, 6, 0, 6, 2, 0],
        [2, 7, 6, 2, 3, 7, 0, 1, 9],
        [1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6],
        [10, 7, 6, 10, 1, 7, 1, 3, 7],
        [10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8],
        [0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7],
        [7, 6, 10, 7, 10, 8, 8, 10, 9],
        [6, 8, 4, 11, 8, 6],
        [3, 6, 11, 3, 0, 6, 0, 4, 6],
        [8, 6, 11, 8, 4, 6, 9, 0, 1],
        [9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6],
        [6, 8, 4, 6, 11, 8, 2, 10, 1],
        [1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6],
        [4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9],
        [10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3],
        [8, 2, 3, 8, 4, 2, 4, 6, 2],
        [0, 4, 2, 4, 6, 2],
        [1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8],
        [1, 9, 4, 1, 4, 2, 2, 4, 6],
        [8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1],
        [10, 1, 0, 10, 0, 6, 6, 0, 4],
        [4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3],
        [10, 9, 4, 6, 10, 4],
        [4, 9, 5, 7, 6, 11],
        [0, 8, 3, 4, 9, 5, 11, 7, 6],
        [5, 0, 1, 5, 4, 0, 7, 6, 11],
        [11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5],
        [9, 5, 4, 10, 1, 2, 7, 6, 11],
        [6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5],
        [7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2],
        [3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6],
        [7, 2, 3, 7, 6, 2, 5, 4, 9],
        [9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7],
        [3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0],
        [6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8],
        [9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7],
        [1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4],
        [4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10],
        [7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10],
        [6, 9, 5, 6, 11, 9, 11, 8, 9],
        [3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5],
        [0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11],
        [6, 11, 3, 6, 3, 5, 5, 3, 1],
        [1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6],
        [0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10],
        [11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5],
        [6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3],
        [5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2],
        [9, 5, 6, 9, 6, 0, 0, 6, 2],
        [1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8],
        [1, 5, 6, 2, 1, 6],
        [1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6],
        [10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0],
        [0, 3, 8, 5, 6, 10],
        [10, 5, 6],
        [11, 5, 10, 7, 5, 11],
        [11, 5, 10, 11, 7, 5, 8, 3, 0],
        [5, 11, 7, 5, 10, 11, 1, 9, 0],
        [10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1],
        [11, 1, 2, 11, 7, 1, 7, 5, 1],
        [0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11],
        [9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7],
        [7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2],
        [2, 5, 10, 2, 3, 5, 3, 7, 5],
        [8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5],
        [9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2],
        [9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2],
        [1, 3, 5, 3, 7, 5],
        [0, 8, 7, 0, 7, 1, 1, 7, 5],
        [9, 0, 3, 9, 3, 5, 5, 3, 7],
        [9, 8, 7, 5, 9, 7],
        [5, 8, 4, 5, 10, 8, 10, 11, 8],
        [5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0],
        [0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5],
        [10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4],
        [2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8],
        [0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11],
        [0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5],
        [9, 4, 5, 2, 11, 3],
        [2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4],
        [5, 10, 2, 5, 2, 4, 4, 2, 0],
        [3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9],
        [5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2],
        [8, 4, 5, 8, 5, 3, 3, 5, 1],
        [0, 4, 5, 1, 0, 5],
        [8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5],
        [9, 4, 5],
        [4, 11, 7, 4, 9, 11, 9, 10, 11],
        [0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11],
        [1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11],
        [3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4],
        [4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2],
        [9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3],
        [11, 7, 4, 11, 4, 2, 2, 4, 0],
        [11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4],
        [2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9],
        [9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7],
        [3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10],
        [1, 10, 2, 8, 7, 4],
        [4, 9, 1, 4, 1, 7, 7, 1, 3],
        [4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1],
        [4, 0, 3, 7, 4, 3],
        [4, 8, 7],
        [9, 10, 8, 10, 11, 8],
        [3, 0, 9, 3, 9, 11, 11, 9, 10],
        [0, 1, 10, 0, 10, 8, 8, 10, 11],
        [3, 1, 10, 11, 3, 10],
        [1, 2, 11, 1, 11, 9, 9, 11, 8],
        [3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9],
        [0, 2, 11, 8, 0, 11],
        [3, 2, 11],
        [2, 3, 8, 2, 8, 10, 10, 8, 9],
        [9, 10, 2, 0, 9, 2],
        [2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8],
        [1, 10, 2],
        [1, 3, 8, 9, 1, 8],
        [0, 9, 1],
        [0, 3, 8],
        []
    ]

    configuration = lookup_table[case]

    return configuration


def _get_vertex_positions(configuration, coordinate, stepsize):
    """ Get the vertex positions corresponding to a list of edge indices.
    
    Parameters
    ----------
    configuration : list[int]
        List of edge indices
    coordinate : tuple(float, float, float)
        Bottom-left corner of cube
    stepsize : float
        Size length of cube

    Returns
    -------
    configuration_vertices : list[int]
        List of triangle vertex coordinates

    """

    x, y, z = coordinate  # Unpack coordinates

    # Relative position on the cube of each edge's midpoint
    edge_midpoints = [
        (0.5, 0.0, 0.0),
        (1.0, 0.0, 0.5),
        (0.5, 0.0, 1.0),
        (0.0, 0.0, 0.5),
        (0.5, 1.0, 0.0),
        (1.0, 1.0, 0.5),
        (0.5, 1.0, 1.0),
        (0.0, 1.0, 0.5),
        (0.0, 0.5, 0.0),
        (1.0, 0.5, 0.0),
        (1.0, 0.5, 1.0),
        (0.0, 0.5, 1.0)
    ]

    # Calculate vertex coordinate for each edge index in the configuration
    configuration_vertices = []
    for index in configuration:

        edge_x, edge_y, edge_z = edge_midpoints[index]

        configuration_vertices.append(x + stepsize * edge_x)
        configuration_vertices.append(y + stepsize * edge_y)
        configuration_vertices.append(z + stepsize * edge_z)

    return configuration_vertices


# Main functions
def marching_cubes(scalar_field, isovalue, volume_min, volume_max, stepsize):
    """ Generates a triangle mesh of an object from a scalar field using marching cubes.
    
    Parameters
    ----------
    scalar_field : function
        Scalar field f(x, y, z) defining the object to generate a triangle mesh of
    isovalue : float
        Boundary value used to define which points are considered to be inside
        of the object. If f(x, y, z) < isovalue, then the point (x, y, z) is
        considered to be inside of the object. If f(x, y, z) > isovalue, then
        the point (x, y, z) is considered to be outside of the object.
    volume_min : int
        Minimum value of x, y, and z to use in the algorithm.
    volume_max : int
        Maximum value of x, y, and z to use in the algorithm.
    stepsize : float
        Side length of each cube in the algorithm.

    Returns
    -------
    mesh_vertices : list[float]
        Vertices in the triangle mesh for the object generated by the algorithm
    
    """

    mesh_vertices = []  # Vertices in the triangle mesh for the object generated by the algorithm

    # Indices in the lookup table corresponding each cube corner. 
    # Used to build the bitmask for each case of marching cubes
    FRONT_BOTTOM_LEFT = 8  
    FRONT_BOTTOM_RIGHT = 4
    FRONT_TOP_LEFT = 128
    FRONT_TOP_RIGHT = 64

    BACK_BOTTOM_LEFT = 1
    BACK_BOTTOM_RIGHT = 2
    BACK_TOP_LEFT = 16
    BACK_TOP_RIGHT = 32

    # Divide the volume into cubes
    for y in np.arange(volume_min, volume_max, stepsize):
        for x in np.arange(volume_min, volume_max, stepsize):
            for z in np.arange(volume_min, volume_max, stepsize):

                # Test if each of the 8 corners of the cube fall inside or
                # outside of the object, and build up bitmask
                case = 0
                if (scalar_field(x, y, z) < isovalue):
                    case |= BACK_BOTTOM_LEFT
                if (scalar_field(x+stepsize, y, z) < isovalue):
                    case |= BACK_BOTTOM_RIGHT
                if (scalar_field(x, y+stepsize, z) < isovalue):
                    case |= BACK_TOP_LEFT
                if (scalar_field(x+stepsize, y+stepsize, z) < isovalue):
                    case |= BACK_TOP_RIGHT
                if (scalar_field(x, y, z+stepsize) < isovalue):
                    case |= FRONT_BOTTOM_LEFT
                if (scalar_field(x+stepsize, y, z+stepsize) < isovalue):
                    case |= FRONT_BOTTOM_RIGHT
                if (scalar_field(x, y+stepsize, z+stepsize) < isovalue):
                    case |= FRONT_TOP_LEFT
                if (scalar_field(x+stepsize, y+stepsize, z+stepsize) < isovalue):
                    case |= FRONT_TOP_RIGHT

                edge_indices = _lookup_configuration(case)

                cube_vertices = _get_vertex_positions(edge_indices, (x, y, z), stepsize)

                for vertex in cube_vertices:
                    mesh_vertices.append(vertex)

    return mesh_vertices


def compute_normals(vertices):
    """ Computes the normal vectors for each vertex in a triangle mesh.

    Parameters
    ----------
    vertices : list[float]
        Vertices in the triangle mesh
    
    Returns
    -------
    normals : list[float]
        Normal vectors for each vertex in the triangle mesh
    
    """

    normals = []  # Normal vectors for each vertex in the triangle mesh

    # Iterate over each triangle in the triangle mesh
    for index in np.arange(0, len(vertices), 9):
        
        # Define each vertex of the triangle
        vertexA = glm.vec3(vertices[index], vertices[index+1], vertices[index+2])
        vertexB = glm.vec3(vertices[index+3], vertices[index+4], vertices[index+5])
        vertexC = glm.vec3(vertices[index+6], vertices[index+7], vertices[index+8])

        # Define triangle edges
        edgeA = vertexB - vertexA
        edgeB = vertexC - vertexA

        # Compute cross product and normalize for normal vector
        normal = glm.normalize(glm.cross(edgeA, edgeB))

        # Add normal to the list 3 times (once for each vertex of the triangle)
        for i in range(3):
            normals.append(normal.x)
            normals.append(normal.y)
            normals.append(normal.z)

    return normals


def writePLY(vertices, normals, filename, comment=None):
    """ Writes a triangle mesh to a PLY file.

    Parameters
    ----------
    vertices : list[int]
        Vertices in the triangle mesh
    normals : list[int]
        Normal vectors for each vertex in the triangle mesh
    filename : str
        Filename of PLY file to write triangle mesh to
    comment : str
        Comment line to write in PLY file
    
    """

    # Get filepath to file
    directory = os.path.dirname(os.path.abspath(__file__))
    filepath = os.path.join(directory, filename)

    num_vertices = int(len(vertices)/3)   # Each vertex has 3 coordinates (x, y, z)
    num_faces = int(num_vertices/3)       # 3 vertices make up 1 triangle face

    # Write data to file
    with open(filepath, "w") as file:
        
        # Header
        file.write("ply\n")
        file.write("format ascii 1.0\n")
        file.write("comment {}\n".format(comment))

        # Vertex metadata
        file.write("element vertex {}\n".format(num_vertices))  
        file.write("property float x\nproperty float y\nproperty float z\n")
        file.write("property float nx\nproperty float ny\nproperty float nz\n")
        
        # Triangle face metadata
        file.write("element face {}\n".format(num_faces))  
        file.write("property list uchar uint vertex_indices\n")

        file.write("end_header\n")

        # Vertex data
        for i in np.arange(0, len(vertices), 3):  # Iterate over each vertex
            file.write("{} {} {} {} {} {}\n".format(vertices[i], vertices[i+1], vertices[i+2],
                                                    normals[i], normals[i+1], normals[i+2]))

        # Triangle face indices
        for i in np.arange(0, num_vertices, 3):  # Iterate over groups of vertex indices forming each face
            file.write("3 {} {} {}\n".format(i, i+1, i+2))


# Validate input and unpack arguments
if (len(sys.argv) == 7):
    _, filename, scalar_field_name, isovalue, volume_min, volume_max, stepsize = sys.argv
else:  # Invalid input
    raise TypeError("Expected 6 arguments but {} were given.".format(len(sys.argv) - 1))


# Define scalar field function to use from input
if (scalar_field_name == "wavy"):
    scalar_field = wavy
elif (scalar_field_name == "hyperboloid"):
    scalar_field = hyperboloid
else:  # Invalid input
    raise ValueError("Undefined scalar field: '{}'. Please enter one of the following options:\n'wavy', 'hyperboloid'".format(scalar_field_name))


# Run marching cubes algorithm and save result
vertices = marching_cubes(scalar_field, float(isovalue), int(volume_min), int(volume_max), float(stepsize))  # Generate vertex positions
normals = compute_normals(vertices)    # Compute normals

# Save to file
comment = "Scalar field: '{}', Isovalue: {}, Min. volume: {}, Max. volume: {}, Stepsize: {}".format(scalar_field_name, isovalue, volume_min, volume_max, stepsize)
writePLY(vertices, normals, filename, comment)

print("Marching cubes computation complete!")