// David Eberly, Geometric Tools, Redmond WA 98052
// Copyright (c) 1998-2020
// Distributed under the Boost Software License, Version 1.0.
// https://www.boost.org/LICENSE_1_0.txt
// https://www.geometrictools.com/License/Boost/LICENSE_1_0.txt
// Version: 4.0.2019.08.13

#pragma once

#include <Mathematics/SurfaceExtractor.h>
#include <set>

// The level set extraction algorithm implemented here is described
// in the document
// https://www.geometrictools.com/Documentation/ExtractLevelSurfaces.pdf

namespace WwiseGTE
{
    // The image type T must be one of the integer types:  int8_t, int16_t,
    // int32_t, uint8_t, uint16_t or uint32_t.  Internal integer computations
    // are performed using int64_t.  The type Real is for extraction to
    // floating-point vertices.
    template <typename T, typename Real>
    class SurfaceExtractorTetrahedra : public SurfaceExtractor<T, Real>
    {
    public:
        // Convenience type definitions.
        typedef typename SurfaceExtractor<T, Real>::Vertex Vertex;
        typedef typename SurfaceExtractor<T, Real>::Triangle Triangle;

        // The input is a 3D image with lexicographically ordered voxels
        // (x,y,z) stored in a linear array.  Voxel (x,y,z) is stored in the
        // array at location index = x + xBound * (y + yBound * z).  The
        // inputs xBound, yBound and zBound must each be 2 or larger so that
        // there is at least one image cube to process.  The inputVoxels must
        // be nonnull and point to contiguous storage that contains at least
        // xBound * yBound * zBound elements.
        SurfaceExtractorTetrahedra(int xBound, int yBound, int zBound, T const* inputVoxels)
            :
            SurfaceExtractor<T, Real>(xBound, yBound, zBound, inputVoxels),
            mNextVertex(0)
        {
        }

        // Extract level surfaces and return rational vertices.  Use the
        // base-class Extract if you want real-valued vertices.
        virtual void Extract(T level, std::vector<Vertex>& vertices,
            std::vector<Triangle>& triangles) override
        {
            // Adjust the image so that the level set is F(x,y,z) = 0.
            int64_t levelI64 = static_cast<int64_t>(level);
            for (size_t i = 0; i < this->mVoxels.size(); ++i)
            {
                int64_t inputI64 = static_cast<int64_t>(this->mInputVoxels[i]);
                this->mVoxels[i] = inputI64 - levelI64;
            }

            mVMap.clear();
            mESet.clear();
            mTSet.clear();
            mNextVertex = 0;
            vertices.clear();
            triangles.clear();
            for (int z = 0, zp = 1; zp < this->mZBound; ++z, ++zp)
            {
                int zParity = (z & 1);
                for (int y = 0, yp = 1; yp < this->mYBound; ++y, ++yp)
                {
                    int yParity = (y & 1);
                    for (int x = 0, xp = 1; xp < this->mXBound; ++x, ++xp)
                    {
                        int xParity = (x & 1);

                        int i000 = x + this->mXBound * (y + this->mYBound * z);
                        int i100 = i000 + 1;
                        int i010 = i000 + this->mXBound;
                        int i110 = i010 + 1;
                        int i001 = i000 + this->mXYBound;
                        int i101 = i001 + 1;
                        int i011 = i001 + this->mXBound;
                        int i111 = i011 + 1;
                        int64_t f000 = static_cast<int64_t>(this->mVoxels[i000]);
                        int64_t f100 = static_cast<int64_t>(this->mVoxels[i100]);
                        int64_t f010 = static_cast<int64_t>(this->mVoxels[i010]);
                        int64_t f110 = static_cast<int64_t>(this->mVoxels[i110]);
                        int64_t f001 = static_cast<int64_t>(this->mVoxels[i001]);
                        int64_t f101 = static_cast<int64_t>(this->mVoxels[i101]);
                        int64_t f011 = static_cast<int64_t>(this->mVoxels[i011]);
                        int64_t f111 = static_cast<int64_t>(this->mVoxels[i111]);

                        if (xParity ^ yParity ^ zParity)
                        {
                            // 1205
                            ProcessTetrahedron(
                                xp, y, z, f100,
                                xp, yp, z, f110,
                                x, y, z, f000,
                                xp, y, zp, f101);

                            // 3027
                            ProcessTetrahedron(
                                x, yp, z, f010,
                                x, y, z, f000,
                                xp, yp, z, f110,
                                x, yp, zp, f011);

                            // 4750
                            ProcessTetrahedron(
                                x, y, zp, f001,
                                x, yp, zp, f011,
                                xp, y, zp, f101,
                                x, y, z, f000);

                            // 6572
                            ProcessTetrahedron(
                                xp, yp, zp, f111,
                                xp, y, zp, f101,
                                x, yp, zp, f011,
                                xp, yp, z, f110);

                            // 0752
                            ProcessTetrahedron(
                                x, y, z, f000,
                                x, yp, zp, f011,
                                xp, y, zp, f101,
                                xp, yp, z, f110);
                        }
                        else
                        {
                            // 0134
                            ProcessTetrahedron(
                                x, y, z, f000,
                                xp, y, z, f100,
                                x, yp, z, f010,
                                x, y, zp, f001);

                            // 2316
                            ProcessTetrahedron(
                                xp, yp, z, f110,
                                x, yp, z, f010,
                                xp, y, z, f100,
                                xp, yp, zp, f111);

                            // 5461
                            ProcessTetrahedron(
                                xp, y, zp, f101,
                                x, y, zp, f001,
                                xp, yp, zp, f111,
                                xp, y, z, f100);

                            // 7643
                            ProcessTetrahedron(
                                x, yp, zp, f011,
                                xp, yp, zp, f111,
                                x, y, zp, f001,
                                x, yp, z, f010);

                            // 6314
                            ProcessTetrahedron(
                                xp, yp, zp, f111,
                                x, yp, z, f010,
                                xp, y, z, f100,
                                x, y, zp, f001);
                        }
                    }
                }
            }

            // Pack vertices into an array.
            vertices.resize(mVMap.size());
            for (auto const& element : mVMap)
            {
                vertices[element.second] = element.first;
            }

            // Pack edges into an array (computed, but not reported to
            // caller).
            std::vector<Edge> edges(mESet.size());
            size_t i = 0;
            for (auto const& element : mESet)
            {
                edges[i++] = element;
            }

            // Pack triangles into an array.
            triangles.resize(mTSet.size());
            i = 0;
            for (auto const& element : mTSet)
            {
                triangles[i++] = element;
            }
        }

    protected:
        struct Edge
        {
            Edge() = default;

            Edge(int v0, int v1)
            {
                if (v0 < v1)
                {
                    v[0] = v0;
                    v[1] = v1;
                }
                else
                {
                    v[0] = v1;
                    v[1] = v0;
                }
            }

            bool operator==(Edge const& other) const
            {
                return v[0] == other.v[0] && v[1] == other.v[1];
            }

            bool operator<(Edge const& other) const
            {
                for (int i = 0; i < 2; ++i)
                {
                    if (v[i] < other.v[i])
                    {
                        return true;
                    }
                    if (v[i] > other.v[i])
                    {
                        return false;
                    }
                }
                return false;
            }

            std::array<int, 2> v;
        };

        virtual std::array<Real, 3> GetGradient(std::array<Real, 3> const& pos) override
        {
            std::array<Real, 3> const zero{ (Real)0, (Real)0, (Real)0 };

            int x = static_cast<int>(pos[0]);
            if (x < 0 || x + 1 >= this->mXBound)
            {
                return zero;
            }

            int y = static_cast<int>(pos[1]);
            if (y < 0 || y + 1 >= this->mYBound)
            {
                return zero;
            }

            int z = static_cast<int>(pos[2]);
            if (z < 0 || z + 1 >= this->mZBound)
            {
                return zero;
            }

            // Get image values at corners of voxel.
            int i000 = x + this->mXBound * (y + this->mYBound * z);
            int i100 = i000 + 1;
            int i010 = i000 + this->mXBound;
            int i110 = i010 + 1;
            int i001 = i000 + this->mXYBound;
            int i101 = i001 + 1;
            int i011 = i001 + this->mXBound;
            int i111 = i011 + 1;
            Real f000 = static_cast<Real>(this->mVoxels[i000]);
            Real f100 = static_cast<Real>(this->mVoxels[i100]);
            Real f010 = static_cast<Real>(this->mVoxels[i010]);
            Real f110 = static_cast<Real>(this->mVoxels[i110]);
            Real f001 = static_cast<Real>(this->mVoxels[i001]);
            Real f101 = static_cast<Real>(this->mVoxels[i101]);
            Real f011 = static_cast<Real>(this->mVoxels[i011]);
            Real f111 = static_cast<Real>(this->mVoxels[i111]);

            Real dx = pos[0] - static_cast<Real>(x);
            Real dy = pos[1] - static_cast<Real>(y);
            Real dz = pos[2] - static_cast<Real>(z);

            std::array<Real, 3> grad;

            if ((x & 1) ^ (y & 1) ^ (z & 1))
            {
                if (dx - dy - dz >= (Real)0)
                {
                    // 1205
                    grad[0] = +f100 - f000;
                    grad[1] = -f100 + f110;
                    grad[2] = -f100 + f101;
                }
                else if (dx - dy + dz <= (Real)0)
                {
                    // 3027
                    grad[0] = -f010 + f110;
                    grad[1] = +f010 - f000;
                    grad[2] = -f010 + f011;
                }
                else if (dx + dy - dz <= (Real)0)
                {
                    // 4750
                    grad[0] = -f001 + f101;
                    grad[1] = -f001 + f011;
                    grad[2] = +f001 - f000;
                }
                else if (dx + dy + dz >= (Real)0)
                {
                    // 6572
                    grad[0] = +f111 - f011;
                    grad[1] = +f111 - f101;
                    grad[2] = +f111 - f110;
                }
                else
                {
                    // 0752
                    grad[0] = (Real)0.5 * (-f000 - f011 + f101 + f110);
                    grad[1] = (Real)0.5 * (-f000 + f011 - f101 + f110);
                    grad[2] = (Real)0.5 * (-f000 + f011 + f101 - f110);
                }
            }
            else
            {
                if (dx + dy + dz <= (Real)1)
                {
                    // 0134
                    grad[0] = -f000 + f100;
                    grad[1] = -f000 + f010;
                    grad[2] = -f000 + f001;
                }
                else if (dx + dy - dz >= (Real)1)
                {
                    // 2316
                    grad[0] = +f110 - f010;
                    grad[1] = +f110 - f100;
                    grad[2] = -f110 + f111;
                }
                else if (dx - dy + dz >= (Real)1)
                {
                    // 5461
                    grad[0] = +f101 - f001;
                    grad[1] = -f101 + f111;
                    grad[2] = +f101 - f100;
                }
                else if (-dx + dy + dz >= (Real)1)
                {
                    // 7643
                    grad[0] = -f011 + f111;
                    grad[1] = +f011 - f001;
                    grad[2] = +f011 - f010;
                }
                else
                {
                    // 6314
                    grad[0] = (Real)0.5 * (f111 - f010 + f100 - f001);
                    grad[1] = (Real)0.5 * (f111 + f010 - f100 - f001);
                    grad[2] = (Real)0.5 * (f111 - f010 - f100 + f001);
                }
            }

            return grad;
        }

        int AddVertex(Vertex const& v)
        {
            auto iter = mVMap.find(v);
            if (iter != mVMap.end())
            {
                // Vertex already in map, just return its unique index.
                return iter->second;
            }
            else
            {
                // Vertex not in map, insert it and assign it a unique index.
                int i = mNextVertex++;
                mVMap.insert(std::make_pair(v, i));
                return i;
            }
        }

        void AddEdge(Vertex const& v0, Vertex const& v1)
        {
            int i0 = AddVertex(v0);
            int i1 = AddVertex(v1);
            mESet.insert(Edge(i0, i1));
        }

        void AddTriangle(Vertex const& v0, Vertex const& v1, Vertex const& v2)
        {
            int i0 = AddVertex(v0);
            int i1 = AddVertex(v1);
            int i2 = AddVertex(v2);

            // Nothing to do if triangle already exists.
            Triangle triangle(i0, i1, i2);
            if (mTSet.find(triangle) != mTSet.end())
            {
                return;
            }

            // Prevent double-sided triangles.
            std::swap(triangle.v[1], triangle.v[2]);
            if (mTSet.find(triangle) != mTSet.end())
            {
                return;
            }

            mESet.insert(Edge(i0, i1));
            mESet.insert(Edge(i1, i2));
            mESet.insert(Edge(i2, i0));

            mTSet.insert(triangle);
        }

        // Support for extraction with linear interpolation.
        void ProcessTetrahedron(
            int64_t x0, int64_t y0, int64_t z0, int64_t f0,
            int64_t x1, int64_t y1, int64_t z1, int64_t f1,
            int64_t x2, int64_t y2, int64_t z2, int64_t f2,
            int64_t x3, int64_t y3, int64_t z3, int64_t f3)
        {
            int64_t xn0, yn0, zn0, d0;
            int64_t xn1, yn1, zn1, d1;
            int64_t xn2, yn2, zn2, d2;
            int64_t xn3, yn3, zn3, d3;
            if (f0 != 0)
            {
                // convert to case +***
                if (f0 < 0)
                {
                    f0 = -f0;
                    f1 = -f1;
                    f2 = -f2;
                    f3 = -f3;
                }

                if (f1 > 0)
                {
                    if (f2 > 0)
                    {
                        if (f3 > 0)
                        {
                            // ++++
                            return;
                        }
                        else if (f3 < 0)
                        {
                            // +++-
                            d0 = f0 - f3;
                            xn0 = f0 * x3 - f3 * x0;
                            yn0 = f0 * y3 - f3 * y0;
                            zn0 = f0 * z3 - f3 * z0;
                            d1 = f1 - f3;
                            xn1 = f1 * x3 - f3 * x1;
                            yn1 = f1 * y3 - f3 * y1;
                            zn1 = f1 * z3 - f3 * z1;
                            d2 = f2 - f3;
                            xn2 = f2 * x3 - f3 * x2;
                            yn2 = f2 * y3 - f3 * y2;
                            zn2 = f2 * z3 - f3 * z2;
                            AddTriangle(
                                Vertex(xn0, d0, yn0, d0, zn0, d0),
                                Vertex(xn1, d1, yn1, d1, zn1, d1),
                                Vertex(xn2, d2, yn2, d2, zn2, d2));
                        }
                        else
                        {
                            // +++0
                            AddVertex(
                                Vertex(x3, 1, y3, 1, z3, 1));
                        }
                    }
                    else if (f2 < 0)
                    {
                        d0 = f0 - f2;
                        xn0 = f0 * x2 - f2 * x0;
                        yn0 = f0 * y2 - f2 * y0;
                        zn0 = f0 * z2 - f2 * z0;
                        d1 = f1 - f2;
                        xn1 = f1 * x2 - f2 * x1;
                        yn1 = f1 * y2 - f2 * y1;
                        zn1 = f1 * z2 - f2 * z1;

                        if (f3 > 0)
                        {
                            // ++-+
                            d2 = f3 - f2;
                            xn2 = f3 * x2 - f2 * x3;
                            yn2 = f3 * y2 - f2 * y3;
                            zn2 = f3 * z2 - f2 * z3;
                            AddTriangle(
                                Vertex(xn0, d0, yn0, d0, zn0, d0),
                                Vertex(xn1, d1, yn1, d1, zn1, d1),
                                Vertex(xn2, d2, yn2, d2, zn2, d2));
                        }
                        else if (f3 < 0)
                        {
                            // ++--
                            d2 = f0 - f3;
                            xn2 = f0 * x3 - f3 * x0;
                            yn2 = f0 * y3 - f3 * y0;
                            zn2 = f0 * z3 - f3 * z0;
                            d3 = f1 - f3;
                            xn3 = f1 * x3 - f3 * x1;
                            yn3 = f1 * y3 - f3 * y1;
                            zn3 = f1 * z3 - f3 * z1;
                            AddTriangle(
                                Vertex(xn0, d0, yn0, d0, zn0, d0),
                                Vertex(xn1, d1, yn1, d1, zn1, d1),
                                Vertex(xn2, d2, yn2, d2, zn2, d2));
                            AddTriangle(
                                Vertex(xn1, d1, yn1, d1, zn1, d1),
                                Vertex(xn3, d3, yn3, d3, zn3, d3),
                                Vertex(xn2, d2, yn2, d2, zn2, d2));
                        }
                        else
                        {
                            // ++-0
                            AddTriangle(
                                Vertex(xn0, d0, yn0, d0, zn0, d0),
                                Vertex(xn1, d1, yn1, d1, zn1, d1),
                                Vertex(x3, 1, y3, 1, z3, 1));
                        }
                    }
                    else
                    {
                        if (f3 > 0)
                        {
                            // ++0+
                            AddVertex(
                                Vertex(x2, 1, y2, 1, z2, 1));
                        }
                        else if (f3 < 0)
                        {
                            // ++0-
                            d0 = f0 - f3;
                            xn0 = f0 * x3 - f3 * x0;
                            yn0 = f0 * y3 - f3 * y0;
                            zn0 = f0 * z3 - f3 * z0;
                            d1 = f1 - f3;
                            xn1 = f1 * x3 - f3 * x1;
                            yn1 = f1 * y3 - f3 * y1;
                            zn1 = f1 * z3 - f3 * z1;
                            AddTriangle(
                                Vertex(xn0, d0, yn0, d0, zn0, d0),
                                Vertex(xn1, d1, yn1, d1, zn1, d1),
                                Vertex(x2, 1, y2, 1, z2, 1));
                        }
                        else
                        {
                            // ++00
                            AddEdge(
                                Vertex(x2, 1, y2, 1, z2, 1),
                                Vertex(x3, 1, y3, 1, z3, 1));
                        }
                    }
                }
                else if (f1 < 0)
                {
                    if (f2 > 0)
                    {
                        d0 = f0 - f1;
                        xn0 = f0 * x1 - f1 * x0;
                        yn0 = f0 * y1 - f1 * y0;
                        zn0 = f0 * z1 - f1 * z0;
                        d1 = f2 - f1;
                        xn1 = f2 * x1 - f1 * x2;
                        yn1 = f2 * y1 - f1 * y2;
                        zn1 = f2 * z1 - f1 * z2;

                        if (f3 > 0)
                        {
                            // +-++
                            d2 = f3 - f1;
                            xn2 = f3 * x1 - f1 * x3;
                            yn2 = f3 * y1 - f1 * y3;
                            zn2 = f3 * z1 - f1 * z3;
                            AddTriangle(
                                Vertex(xn0, d0, yn0, d0, zn0, d0),
                                Vertex(xn1, d1, yn1, d1, zn1, d1),
                                Vertex(xn2, d2, yn2, d2, zn2, d2));
                        }
                        else if (f3 < 0)
                        {
                            // +-+-
                            d2 = f0 - f3;
                            xn2 = f0 * x3 - f3 * x0;
                            yn2 = f0 * y3 - f3 * y0;
                            zn2 = f0 * z3 - f3 * z0;
                            d3 = f2 - f3;
                            xn3 = f2 * x3 - f3 * x2;
                            yn3 = f2 * y3 - f3 * y2;
                            zn3 = f2 * z3 - f3 * z2;
                            AddTriangle(
                                Vertex(xn0, d0, yn0, d0, zn0, d0),
                                Vertex(xn1, d1, yn1, d1, zn1, d1),
                                Vertex(xn2, d2, yn2, d2, zn2, d2));
                            AddTriangle(
                                Vertex(xn1, d1, yn1, d1, zn1, d1),
                                Vertex(xn3, d3, yn3, d3, zn3, d3),
                                Vertex(xn2, d2, yn2, d2, zn2, d2));
                        }
                        else
                        {
                            // +-+0
                            AddTriangle(
                                Vertex(xn0, d0, yn0, d0, zn0, d0),
                                Vertex(xn1, d1, yn1, d1, zn1, d1),
                                Vertex(x3, 1, y3, 1, z3, 1));
                        }
                    }
                    else if (f2 < 0)
                    {
                        d0 = f1 - f0;
                        xn0 = f1 * x0 - f0 * x1;
                        yn0 = f1 * y0 - f0 * y1;
                        zn0 = f1 * z0 - f0 * z1;
                        d1 = f2 - f0;
                        xn1 = f2 * x0 - f0 * x2;
                        yn1 = f2 * y0 - f0 * y2;
                        zn1 = f2 * z0 - f0 * z2;

                        if (f3 > 0)
                        {
                            // +--+
                            d2 = f1 - f3;
                            xn2 = f1 * x3 - f3 * x1;
                            yn2 = f1 * y3 - f3 * y1;
                            zn2 = f1 * z3 - f3 * z1;
                            d3 = f2 - f3;
                            xn3 = f2 * x3 - f3 * x2;
                            yn3 = f2 * y3 - f3 * y2;
                            zn3 = f2 * z3 - f3 * z2;
                            AddTriangle(
                                Vertex(xn0, d0, yn0, d0, zn0, d0),
                                Vertex(xn1, d1, yn1, d1, zn1, d1),
                                Vertex(xn2, d2, yn2, d2, zn2, d2));
                            AddTriangle(
                                Vertex(xn1, d1, yn1, d1, zn1, d1),
                                Vertex(xn3, d3, yn3, d3, zn3, d3),
                                Vertex(xn2, d2, yn2, d2, zn2, d2));
                        }
                        else if (f3 < 0)
                        {
                            // +---
                            d2 = f3 - f0;
                            xn2 = f3 * x0 - f0 * x3;
                            yn2 = f3 * y0 - f0 * y3;
                            zn2 = f3 * z0 - f0 * z3;
                            AddTriangle(
                                Vertex(xn0, d0, yn0, d0, zn0, d0),
                                Vertex(xn1, d1, yn1, d1, zn1, d1),
                                Vertex(xn2, d2, yn2, d2, zn2, d2));
                        }
                        else
                        {
                            // +--0
                            AddTriangle(
                                Vertex(xn0, d0, yn0, d0, zn0, d0),
                                Vertex(xn1, d1, yn1, d1, zn1, d1),
                                Vertex(x3, 1, y3, 1, z3, 1));
                        }
                    }
                    else
                    {
                        d0 = f1 - f0;
                        xn0 = f1 * x0 - f0 * x1;
                        yn0 = f1 * y0 - f0 * y1;
                        zn0 = f1 * z0 - f0 * z1;

                        if (f3 > 0)
                        {
                            // +-0+
                            d1 = f1 - f3;
                            xn1 = f1 * x3 - f3 * x1;
                            yn1 = f1 * y3 - f3 * y1;
                            zn1 = f1 * z3 - f3 * z1;
                            AddTriangle(
                                Vertex(xn0, d0, yn0, d0, zn0, d0),
                                Vertex(xn1, d1, yn1, d1, zn1, d1),
                                Vertex(x2, 1, y2, 1, z2, 1));
                        }
                        else if (f3 < 0)
                        {
                            // +-0-
                            d1 = f3 - f0;
                            xn1 = f3 * x0 - f0 * x3;
                            yn1 = f3 * y0 - f0 * y3;
                            zn1 = f3 * z0 - f0 * z3;
                            AddTriangle(
                                Vertex(xn0, d0, yn0, d0, zn0, d0),
                                Vertex(xn1, d1, yn1, d1, zn1, d1),
                                Vertex(x2, 1, y2, 1, z2, 1));
                        }
                        else
                        {
                            // +-00
                            AddTriangle(
                                Vertex(xn0, d0, yn0, d0, zn0, d0),
                                Vertex(x2, 1, y2, 1, z2, 1),
                                Vertex(x3, 1, y3, 1, z3, 1));
                        }
                    }
                }
                else
                {
                    if (f2 > 0)
                    {
                        if (f3 > 0)
                        {
                            // +0++
                            AddVertex(
                                Vertex(x1, 1, y1, 1, z1, 1));
                        }
                        else if (f3 < 0)
                        {
                            // +0+-
                            d0 = f0 - f3;
                            xn0 = f0 * x3 - f3 * x0;
                            yn0 = f0 * y3 - f3 * y0;
                            zn0 = f0 * z3 - f3 * z0;
                            d1 = f2 - f3;
                            xn1 = f2 * x3 - f3 * x2;
                            yn1 = f2 * y3 - f3 * y2;
                            zn1 = f2 * z3 - f3 * z2;
                            AddTriangle(
                                Vertex(xn0, d0, yn0, d0, zn0, d0),
                                Vertex(xn1, d1, yn1, d1, zn1, d1),
                                Vertex(x1, 1, y1, 1, z1, 1));
                        }
                        else
                        {
                            // +0+0
                            AddEdge(
                                Vertex(x1, 1, y1, 1, z1, 1),
                                Vertex(x3, 1, y3, 1, z3, 1));
                        }
                    }
                    else if (f2 < 0)
                    {
                        d0 = f2 - f0;
                        xn0 = f2 * x0 - f0 * x2;
                        yn0 = f2 * y0 - f0 * y2;
                        zn0 = f2 * z0 - f0 * z2;

                        if (f3 > 0)
                        {
                            // +0-+
                            d1 = f2 - f3;
                            xn1 = f2 * x3 - f3 * x2;
                            yn1 = f2 * y3 - f3 * y2;
                            zn1 = f2 * z3 - f3 * z2;
                            AddTriangle(
                                Vertex(xn0, d0, yn0, d0, zn0, d0),
                                Vertex(xn1, d1, yn1, d1, zn1, d1),
                                Vertex(x1, 1, y1, 1, z1, 1));
                        }
                        else if (f3 < 0)
                        {
                            // +0--
                            d1 = f0 - f3;
                            xn1 = f0 * x3 - f3 * x0;
                            yn1 = f0 * y3 - f3 * y0;
                            zn1 = f0 * z3 - f3 * z0;
                            AddTriangle(
                                Vertex(xn0, d0, yn0, d0, zn0, d0),
                                Vertex(xn1, d1, yn1, d1, zn1, d1),
                                Vertex(x1, 1, y1, 1, z1, 1));
                        }
                        else
                        {
                            // +0-0
                            AddTriangle(
                                Vertex(xn0, d0, yn0, d0, zn0, d0),
                                Vertex(x1, 1, y1, 1, z1, 1),
                                Vertex(x3, 1, y3, 1, z3, 1));
                        }
                    }
                    else
                    {
                        if (f3 > 0)
                        {
                            // +00+
                            AddEdge(
                                Vertex(x1, 1, y1, 1, z1, 1),
                                Vertex(x2, 1, y2, 1, z2, 1));
                        }
                        else if (f3 < 0)
                        {
                            // +00-
                            d0 = f0 - f3;
                            xn0 = f0 * x3 - f3 * x0;
                            yn0 = f0 * y3 - f3 * y0;
                            zn0 = f0 * z3 - f3 * z0;
                            AddTriangle(
                                Vertex(xn0, d0, yn0, d0, zn0, d0),
                                Vertex(x1, 1, y1, 1, z1, 1),
                                Vertex(x2, 1, y2, 1, z2, 1));
                        }
                        else
                        {
                            // +000
                            AddTriangle(
                                Vertex(x1, 1, y1, 1, z1, 1),
                                Vertex(x2, 1, y2, 1, z2, 1),
                                Vertex(x3, 1, y3, 1, z3, 1));
                        }
                    }
                }
            }
            else if (f1 != 0)
            {
                // convert to case 0+**
                if (f1 < 0)
                {
                    f1 = -f1;
                    f2 = -f2;
                    f3 = -f3;
                }

                if (f2 > 0)
                {
                    if (f3 > 0)
                    {
                        // 0+++
                        AddVertex(
                            Vertex(x0, 1, y0, 1, z0, 1));
                    }
                    else if (f3 < 0)
                    {
                        // 0++-
                        d0 = f2 - f3;
                        xn0 = f2 * x3 - f3 * x2;
                        yn0 = f2 * y3 - f3 * y2;
                        zn0 = f2 * z3 - f3 * z2;
                        d1 = f1 - f3;
                        xn1 = f1 * x3 - f3 * x1;
                        yn1 = f1 * y3 - f3 * y1;
                        zn1 = f1 * z3 - f3 * z1;
                        AddTriangle(
                            Vertex(xn0, d0, yn0, d0, zn0, d0),
                            Vertex(xn1, d1, yn1, d1, zn1, d1),
                            Vertex(x0, 1, y0, 1, z0, 1));
                    }
                    else
                    {
                        // 0++0
                        AddEdge(
                            Vertex(x0, 1, y0, 1, z0, 1),
                            Vertex(x3, 1, y3, 1, z3, 1));
                    }
                }
                else if (f2 < 0)
                {
                    d0 = f2 - f1;
                    xn0 = f2 * x1 - f1 * x2;
                    yn0 = f2 * y1 - f1 * y2;
                    zn0 = f2 * z1 - f1 * z2;

                    if (f3 > 0)
                    {
                        // 0+-+
                        d1 = f2 - f3;
                        xn1 = f2 * x3 - f3 * x2;
                        yn1 = f2 * y3 - f3 * y2;
                        zn1 = f2 * z3 - f3 * z2;
                        AddTriangle(
                            Vertex(xn0, d0, yn0, d0, zn0, d0),
                            Vertex(xn1, d1, yn1, d1, zn1, d1),
                            Vertex(x0, 1, y0, 1, z0, 1));
                    }
                    else if (f3 < 0)
                    {
                        // 0+--
                        d1 = f1 - f3;
                        xn1 = f1 * x3 - f3 * x1;
                        yn1 = f1 * y3 - f3 * y1;
                        zn1 = f1 * z3 - f3 * z1;
                        AddTriangle(
                            Vertex(xn0, d0, yn0, d0, zn0, d0),
                            Vertex(xn1, d1, yn1, d1, zn1, d1),
                            Vertex(x0, 1, y0, 1, z0, 1));
                    }
                    else
                    {
                        // 0+-0
                        AddTriangle(
                            Vertex(xn0, d0, yn0, d0, zn0, d0),
                            Vertex(x0, 1, y0, 1, z0, 1),
                            Vertex(x3, 1, y3, 1, z3, 1));
                    }
                }
                else
                {
                    if (f3 > 0)
                    {
                        // 0+0+
                        AddEdge(
                            Vertex(x0, 1, y0, 1, z0, 1),
                            Vertex(x2, 1, y2, 1, z2, 1));
                    }
                    else if (f3 < 0)
                    {
                        // 0+0-
                        d0 = f1 - f3;
                        xn0 = f1 * x3 - f3 * x1;
                        yn0 = f1 * y3 - f3 * y1;
                        zn0 = f1 * z3 - f3 * z1;
                        AddTriangle(
                            Vertex(xn0, d0, yn0, d0, zn0, d0),
                            Vertex(x0, 1, y0, 1, z0, 1),
                            Vertex(x2, 1, y2, 1, z2, 1));
                    }
                    else
                    {
                        // 0+00
                        AddTriangle(
                            Vertex(x0, 1, y0, 1, z0, 1),
                            Vertex(x2, 1, y2, 1, z2, 1),
                            Vertex(x3, 1, y3, 1, z3, 1));
                    }
                }
            }
            else if (f2 != 0)
            {
                // convert to case 00+*
                if (f2 < 0)
                {
                    f2 = -f2;
                    f3 = -f3;
                }

                if (f3 > 0)
                {
                    // 00++
                    AddEdge(
                        Vertex(x0, 1, y0, 1, z0, 1),
                        Vertex(x1, 1, y1, 1, z1, 1));
                }
                else if (f3 < 0)
                {
                    // 00+-
                    d0 = f2 - f3;
                    xn0 = f2 * x3 - f3 * x2;
                    yn0 = f2 * y3 - f3 * y2;
                    zn0 = f2 * z3 - f3 * z2;
                    AddTriangle(
                        Vertex(xn0, d0, yn0, d0, zn0, d0),
                        Vertex(x0, 1, y0, 1, z0, 1),
                        Vertex(x1, 1, y1, 1, z1, 1));
                }
                else
                {
                    // 00+0
                    AddTriangle(
                        Vertex(x0, 1, y0, 1, z0, 1),
                        Vertex(x1, 1, y1, 1, z1, 1),
                        Vertex(x3, 1, y3, 1, z3, 1));
                }
            }
            else if (f3 != 0)
            {
                // cases 000+ or 000-
                AddTriangle(
                    Vertex(x0, 1, y0, 1, z0, 1),
                    Vertex(x1, 1, y1, 1, z1, 1),
                    Vertex(x2, 1, y2, 1, z2, 1));
            }
            else
            {
                // case 0000
                AddTriangle(
                    Vertex(x0, 1, y0, 1, z0, 1),
                    Vertex(x1, 1, y1, 1, z1, 1),
                    Vertex(x2, 1, y2, 1, z2, 1));
                AddTriangle(
                    Vertex(x0, 1, y0, 1, z0, 1),
                    Vertex(x1, 1, y1, 1, z1, 1),
                    Vertex(x3, 1, y3, 1, z3, 1));
                AddTriangle(
                    Vertex(x0, 1, y0, 1, z0, 1),
                    Vertex(x2, 1, y2, 1, z2, 1),
                    Vertex(x3, 1, y3, 1, z3, 1));
                AddTriangle(
                    Vertex(x1, 1, y1, 1, z1, 1),
                    Vertex(x2, 1, y2, 1, z2, 1),
                    Vertex(x3, 1, y3, 1, z3, 1));
            }
        }

        std::map<Vertex, int> mVMap;
        std::set<Edge> mESet;
        std::set<Triangle> mTSet;
        int mNextVertex;
    };
}