// 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.2020.01.10

#pragma once

#include <Mathematics/Vector.h>

namespace WwiseGTE
{
    // Template alias for convenience.
    template <typename Real>
    using Vector3 = Vector<3, Real>;

    // Cross, UnitCross, and DotCross have a template parameter N that should
    // be 3 or 4.  The latter case supports affine vectors in 4D (last
    // component w = 0) when you want to use 4-tuples and 4x4 matrices for
    // affine algebra.

    // Compute the cross product using the formal determinant:
    //   cross = det{{e0,e1,e2},{x0,x1,x2},{y0,y1,y2}}
    //         = (x1*y2-x2*y1, x2*y0-x0*y2, x0*y1-x1*y0)
    // where e0 = (1,0,0), e1 = (0,1,0), e2 = (0,0,1), v0 = (x0,x1,x2), and
    // v1 = (y0,y1,y2).
    template <int N, typename Real>
    Vector<N, Real> Cross(Vector<N, Real> const& v0, Vector<N, Real> const& v1)
    {
        static_assert(N == 3 || N == 4, "Dimension must be 3 or 4.");

        Vector<N, Real> result;
        result.MakeZero();
        result[0] = v0[1] * v1[2] - v0[2] * v1[1];
        result[1] = v0[2] * v1[0] - v0[0] * v1[2];
        result[2] = v0[0] * v1[1] - v0[1] * v1[0];
        return result;
    }

    // Compute the normalized cross product.
    template <int N, typename Real>
    Vector<N, Real> UnitCross(Vector<N, Real> const& v0, Vector<N, Real> const& v1, bool robust = false)
    {
        static_assert(N == 3 || N == 4, "Dimension must be 3 or 4.");
        Vector<N, Real> unitCross = Cross(v0, v1);
        Normalize(unitCross, robust);
        return unitCross;
    }

    // Compute Dot((x0,x1,x2),Cross((y0,y1,y2),(z0,z1,z2)), the triple scalar
    // product of three vectors, where v0 = (x0,x1,x2), v1 = (y0,y1,y2), and
    // v2 is (z0,z1,z2).
    template <int N, typename Real>
    Real DotCross(Vector<N, Real> const& v0, Vector<N, Real> const& v1,
        Vector<N, Real> const& v2)
    {
        static_assert(N == 3 || N == 4, "Dimension must be 3 or 4.");
        return Dot(v0, Cross(v1, v2));
    }

    // Compute a right-handed orthonormal basis for the orthogonal complement
    // of the input vectors.  The function returns the smallest length of the
    // unnormalized vectors computed during the process.  If this value is
    // nearly zero, it is possible that the inputs are linearly dependent
    // (within numerical round-off errors).  On input, numInputs must be 1 or
    // 2 and v[0] through v[numInputs-1] must be initialized.  On output, the
    // vectors v[0] through v[2] form an orthonormal set.
    template <typename Real>
    Real ComputeOrthogonalComplement(int numInputs, Vector3<Real>* v, bool robust = false)
    {
        if (numInputs == 1)
        {
            if (std::fabs(v[0][0]) > std::fabs(v[0][1]))
            {
                v[1] = { -v[0][2], (Real)0, +v[0][0] };
            }
            else
            {
                v[1] = { (Real)0, +v[0][2], -v[0][1] };
            }
            numInputs = 2;
        }

        if (numInputs == 2)
        {
            v[2] = Cross(v[0], v[1]);
            return Orthonormalize<3, Real>(3, v, robust);
        }

        return (Real)0;
    }

    // Compute the barycentric coordinates of the point P with respect to the
    // tetrahedron <V0,V1,V2,V3>, P = b0*V0 + b1*V1 + b2*V2 + b3*V3, where
    // b0 + b1 + b2 + b3 = 1.  The return value is 'true' iff {V0,V1,V2,V3} is
    // a linearly independent set.  Numerically, this is measured by
    // |det[V0 V1 V2 V3]| <= epsilon.  The values bary[] are valid only when
    // the return value is 'true' but set to zero when the return value is
    // 'false'.
    template <typename Real>
    bool ComputeBarycentrics(Vector3<Real> const& p, Vector3<Real> const& v0,
        Vector3<Real> const& v1, Vector3<Real> const& v2, Vector3<Real> const& v3,
        Real bary[4], Real epsilon = (Real)0)
    {
        // Compute the vectors relative to V3 of the tetrahedron.
        Vector3<Real> diff[4] = { v0 - v3, v1 - v3, v2 - v3, p - v3 };

        Real det = DotCross(diff[0], diff[1], diff[2]);
        if (det < -epsilon || det > epsilon)
        {
            Real invDet = ((Real)1) / det;
            bary[0] = DotCross(diff[3], diff[1], diff[2]) * invDet;
            bary[1] = DotCross(diff[3], diff[2], diff[0]) * invDet;
            bary[2] = DotCross(diff[3], diff[0], diff[1]) * invDet;
            bary[3] = (Real)1 - bary[0] - bary[1] - bary[2];
            return true;
        }

        for (int i = 0; i < 4; ++i)
        {
            bary[i] = (Real)0;
        }
        return false;
    }

    // Get intrinsic information about the input array of vectors.  The return
    // value is 'true' iff the inputs are valid (numVectors > 0, v is not
    // null, and epsilon >= 0), in which case the class members are valid.
    template <typename Real>
    class IntrinsicsVector3
    {
    public:
        // The constructor sets the class members based on the input set.
        IntrinsicsVector3(int numVectors, Vector3<Real> const* v, Real inEpsilon)
            :
            epsilon(inEpsilon),
            dimension(0),
            maxRange((Real)0),
            origin{ (Real)0, (Real)0, (Real)0 },
            extremeCCW(false)
        {
            min[0] = (Real)0;
            min[1] = (Real)0;
            min[2] = (Real)0;
            direction[0] = { (Real)0, (Real)0, (Real)0 };
            direction[1] = { (Real)0, (Real)0, (Real)0 };
            direction[2] = { (Real)0, (Real)0, (Real)0 };
            extreme[0] = 0;
            extreme[1] = 0;
            extreme[2] = 0;
            extreme[3] = 0;

            if (numVectors > 0 && v && epsilon >= (Real)0)
            {
                // Compute the axis-aligned bounding box for the input vectors.
                // Keep track of the indices into 'vectors' for the current
                // min and max.
                int j, indexMin[3], indexMax[3];
                for (j = 0; j < 3; ++j)
                {
                    min[j] = v[0][j];
                    max[j] = min[j];
                    indexMin[j] = 0;
                    indexMax[j] = 0;
                }

                int i;
                for (i = 1; i < numVectors; ++i)
                {
                    for (j = 0; j < 3; ++j)
                    {
                        if (v[i][j] < min[j])
                        {
                            min[j] = v[i][j];
                            indexMin[j] = i;
                        }
                        else if (v[i][j] > max[j])
                        {
                            max[j] = v[i][j];
                            indexMax[j] = i;
                        }
                    }
                }

                // Determine the maximum range for the bounding box.
                maxRange = max[0] - min[0];
                extreme[0] = indexMin[0];
                extreme[1] = indexMax[0];
                Real range = max[1] - min[1];
                if (range > maxRange)
                {
                    maxRange = range;
                    extreme[0] = indexMin[1];
                    extreme[1] = indexMax[1];
                }
                range = max[2] - min[2];
                if (range > maxRange)
                {
                    maxRange = range;
                    extreme[0] = indexMin[2];
                    extreme[1] = indexMax[2];
                }

                // The origin is either the vector of minimum x0-value, vector
                // of minimum x1-value, or vector of minimum x2-value.
                origin = v[extreme[0]];

                // Test whether the vector set is (nearly) a vector.
                if (maxRange <= epsilon)
                {
                    dimension = 0;
                    for (j = 0; j < 3; ++j)
                    {
                        extreme[j + 1] = extreme[0];
                    }
                    return;
                }

                // Test whether the vector set is (nearly) a line segment.  We
                // need {direction[2],direction[3]} to span the orthogonal
                // complement of direction[0].
                direction[0] = v[extreme[1]] - origin;
                Normalize(direction[0], false);
                if (std::fabs(direction[0][0]) > std::fabs(direction[0][1]))
                {
                    direction[1][0] = -direction[0][2];
                    direction[1][1] = (Real)0;
                    direction[1][2] = +direction[0][0];
                }
                else
                {
                    direction[1][0] = (Real)0;
                    direction[1][1] = +direction[0][2];
                    direction[1][2] = -direction[0][1];
                }
                Normalize(direction[1], false);
                direction[2] = Cross(direction[0], direction[1]);

                // Compute the maximum distance of the points from the line
                // origin + t * direction[0].
                Real maxDistance = (Real)0;
                Real distance, dot;
                extreme[2] = extreme[0];
                for (i = 0; i < numVectors; ++i)
                {
                    Vector3<Real> diff = v[i] - origin;
                    dot = Dot(direction[0], diff);
                    Vector3<Real> proj = diff - dot * direction[0];
                    distance = Length(proj, false);
                    if (distance > maxDistance)
                    {
                        maxDistance = distance;
                        extreme[2] = i;
                    }
                }

                if (maxDistance <= epsilon * maxRange)
                {
                    // The points are (nearly) on the line
                    // origin + t * direction[0].
                    dimension = 1;
                    extreme[2] = extreme[1];
                    extreme[3] = extreme[1];
                    return;
                }

                // Test whether the vector set is (nearly) a planar polygon.
                // The point v[extreme[2]] is farthest from the line:
                // origin + t * direction[0].  The vector
                // v[extreme[2]] - origin is not necessarily perpendicular to
                // direction[0], so project out the direction[0] component so
                // that the result is perpendicular to direction[0].
                direction[1] = v[extreme[2]] - origin;
                dot = Dot(direction[0], direction[1]);
                direction[1] -= dot * direction[0];
                Normalize(direction[1], false);

                // We need direction[2] to span the orthogonal complement of
                // {direction[0],direction[1]}.
                direction[2] = Cross(direction[0], direction[1]);

                // Compute the maximum distance of the points from the plane
                // origin+t0 * direction[0] + t1 * direction[1].
                maxDistance = (Real)0;
                Real maxSign = (Real)0;
                extreme[3] = extreme[0];
                for (i = 0; i < numVectors; ++i)
                {
                    Vector3<Real> diff = v[i] - origin;
                    distance = Dot(direction[2], diff);
                    Real sign = (distance > (Real)0 ? (Real)1 :
                        (distance < (Real)0 ? (Real)-1 : (Real)0));
                    distance = std::fabs(distance);
                    if (distance > maxDistance)
                    {
                        maxDistance = distance;
                        maxSign = sign;
                        extreme[3] = i;
                    }
                }

                if (maxDistance <= epsilon * maxRange)
                {
                    // The points are (nearly) on the plane
                    // origin + t0 * direction[0] + t1 * direction[1].
                    dimension = 2;
                    extreme[3] = extreme[2];
                    return;
                }

                dimension = 3;
                extremeCCW = (maxSign > (Real)0);
                return;
            }
        }

        // A nonnegative tolerance that is used to determine the intrinsic
        // dimension of the set.
        Real epsilon;

        // The intrinsic dimension of the input set, computed based on the
        // nonnegative tolerance mEpsilon.
        int dimension;

        // Axis-aligned bounding box of the input set.  The maximum range is
        // the larger of max[0]-min[0], max[1]-min[1], and max[2]-min[2].
        Real min[3], max[3];
        Real maxRange;

        // Coordinate system.  The origin is valid for any dimension d.  The
        // unit-length direction vector is valid only for 0 <= i < d.  The
        // extreme index is relative to the array of input points, and is also
        // valid only for 0 <= i < d.  If d = 0, all points are effectively
        // the same, but the use of an epsilon may lead to an extreme index
        // that is not zero.  If d = 1, all points effectively lie on a line
        // segment.  If d = 2, all points effectively line on a plane.  If
        // d = 3, the points are not coplanar.
        Vector3<Real> origin;
        Vector3<Real> direction[3];

        // The indices that define the maximum dimensional extents.  The
        // values extreme[0] and extreme[1] are the indices for the points
        // that define the largest extent in one of the coordinate axis
        // directions.  If the dimension is 2, then extreme[2] is the index
        // for the point that generates the largest extent in the direction
        // perpendicular to the line through the points corresponding to
        // extreme[0] and extreme[1].  Furthermore, if the dimension is 3,
        // then extreme[3] is the index for the point that generates the
        // largest extent in the direction perpendicular to the triangle
        // defined by the other extreme points.  The tetrahedron formed by the
        // points V[extreme[0]], V[extreme[1]], V[extreme[2]], and
        // V[extreme[3]] is clockwise or counterclockwise, the condition
        // stored in extremeCCW.
        int extreme[4];
        bool extremeCCW;
    };
}