// Geometric Tools, Inc.
// http://www.geometrictools.com
// Copyright (c) 1998-2005.  All Rights Reserved
//
// The Wild Magic Library (WM3) source code is supplied under the terms of
// the license agreement
//     http://www.geometrictools.com/License/WildMagic3License.pdf
// and may not be copied or disclosed except in accordance with the terms
// of that agreement.

#include "Wm3Delaunay3.h"
#include "Wm3DelPolyhedronFace.h"
#include "Wm3Mapper3.h"
#include "Wm3ETManifoldMesh.h"
#include "Wm3Delaunay2.h"
#include "Wm3Query3.h"
#include "Wm3Query3Int64.h"
#include "Wm3Query3TInteger.h"
#include "Wm3Query3TRational.h"

#include <algorithm>
#include <map>
#include <stack>
#include <vector>

// Indexing for the vertices of the triangle opposite a vertex.  The triangle
// opposite vertex j is
//   <gs_aaiIndex[j][0],gs_aaiIndex[j][1],gs_aaiIndex[j][2]>
// and is listed in counterclockwise order when viewed from outside the
// tetrahedron.
static const int gs_aaiIndex[4][3] = { {1,2,3}, {0,3,2}, {0,1,3}, {0,2,1} };

namespace Wm3
{
//----------------------------------------------------------------------------
template <class Real>
Delaunay3<Real>::Delaunay3 (int iVertexQuantity, Vector3<Real>* akVertex,
    Real fEpsilon, bool bOwner, int eComputationType)
    :
    Delaunay<Real>(iVertexQuantity,fEpsilon,bOwner,eComputationType),
    m_kLineOrigin(Vector3<Real>::ZERO),
    m_kLineDirection(Vector3<Real>::ZERO),
    m_kPlaneOrigin(Vector3<Real>::ZERO)
{
    assert( akVertex );
    m_akVertex = akVertex;
    m_iUniqueVertexQuantity = 0;
    m_akPlaneDirection[0] = Vector3<Real>::ZERO;
    m_akPlaneDirection[1] = Vector3<Real>::ZERO;
    m_akSVertex = 0;
    m_pkQuery = 0;
    m_iPathLast = -1;
    m_aiPath = 0;
    m_iLastFaceV0 = -1;
    m_iLastFaceV1 = -1;
    m_iLastFaceV2 = -1;
    m_iLastFaceOpposite = -1;
    m_iLastFaceOppositeIndex = -1;

    Mapper3<Real> kMapper(m_iVertexQuantity,m_akVertex,fEpsilon);
    if ( kMapper.GetDimension() == 0 )
    {
        // The values of m_iDimension, m_aiIndex, and m_aiAdjacent were
        // already initialized by the Delaunay base class.
        return;
    }

    int i;
    if ( kMapper.GetDimension() == 1 )
    {
        // The set is (nearly) collinear.  The caller is responsible for
        // creating a Delaunay1 object.
        m_iDimension = 1;
        m_kLineOrigin = kMapper.GetOrigin();
        m_kLineDirection = kMapper.GetDirection(0);
        return;
    }

    if ( kMapper.GetDimension() == 2 )
    {
        // The set is (nearly) coplanar.  The caller is responsible for
        // creating a Delaunay2 object.
        m_iDimension = 2;
        m_kPlaneOrigin = kMapper.GetOrigin();
        m_akPlaneDirection[0] = kMapper.GetDirection(0);
        m_akPlaneDirection[1] = kMapper.GetDirection(1);
        return;
    }

    m_iDimension = 3;

    // Allocate storage for the input vertices and the supertetrahedron
    // vertices.
    m_akSVertex = new Vector3<Real>[m_iVertexQuantity+4];

    if (eComputationType != CT_RATIONAL)
    {
        // Transform the vertices to the cube [0,1]^3.
        m_kMin = kMapper.GetMin();
        m_fScale = ((Real)1.0)/kMapper.GetMaxRange();
        for (i = 0; i < m_iVertexQuantity; i++)
        {
            m_akSVertex[i] = (m_akVertex[i] - m_kMin)*m_fScale;
        }

        // Construct the supertetrahedron to contain [0,1]^3.
        m_aiSV[0] = m_iVertexQuantity++;
        m_aiSV[1] = m_iVertexQuantity++;
        m_aiSV[2] = m_iVertexQuantity++;
        m_aiSV[3] = m_iVertexQuantity++;
        m_akSVertex[m_aiSV[0]] = Vector3<Real>((Real)-1.0,(Real)-1.0,
            (Real)-1.0);
        m_akSVertex[m_aiSV[1]] = Vector3<Real>((Real)+6.0,(Real)-1.0,
            (Real)-1.0);
        m_akSVertex[m_aiSV[2]] = Vector3<Real>((Real)-1.0,(Real)+6.0,
            (Real)-1.0);
        m_akSVertex[m_aiSV[3]] = Vector3<Real>((Real)-1.0,(Real)-1.0,
            (Real)+6.0);

        Real fExpand;
        if (eComputationType == CT_INT64)
        {
            // Scale the vertices to the cube [0,2^{10}]^3 to allow use of
            // 64-bit integers for tetrahedralization.
            fExpand = (Real)(1 << 10);
            m_pkQuery = new Query3Int64<Real>(m_iVertexQuantity,m_akSVertex);
        }
        else if (eComputationType == CT_INTEGER)
        {
            // Scale the vertices to the cube [0,2^{20}]^3 to get more
            // precision for TInteger than for 64-bit integers for
            // tetrahedralization.
            fExpand = (Real)(1 << 20);
            m_pkQuery = new Query3TInteger<Real>(m_iVertexQuantity,
                m_akSVertex);
        }
        else // eComputationType == CT_REAL
        {
            // No scaling for floating point.
            fExpand = (Real)1.0;
            m_pkQuery = new Query3<Real>(m_iVertexQuantity,m_akSVertex);
        }

        m_fScale *= fExpand;
        for (i = 0; i < m_iVertexQuantity; i++)
        {
            m_akSVertex[i] *= fExpand;
        }
    }
    else
    {
        // No transformation needed for exact rational arithmetic.
        m_kMin = Vector3<Real>::ZERO;
        m_fScale = (Real)1.0;
        size_t uiSize = m_iVertexQuantity*sizeof(Vector3<Real>);
        memcpy(m_akSVertex,m_akVertex,uiSize);

        // Construct the supertriangle to contain [min,max].
        Vector3<Real> kMin = kMapper.GetMin();
        Vector3<Real> kMax = kMapper.GetMax();
        Vector3<Real> kDelta = kMax - kMin;
        Vector3<Real> kSMin = kMin - kDelta;
        Vector3<Real> kSMax = kMax + ((Real)5.0)*kDelta;
        m_aiSV[0] = m_iVertexQuantity++;
        m_aiSV[1] = m_iVertexQuantity++;
        m_aiSV[2] = m_iVertexQuantity++;
        m_aiSV[3] = m_iVertexQuantity++;
        m_akSVertex[m_aiSV[0]] = kSMin;
        m_akSVertex[m_aiSV[1]] = Vector3<Real>(kSMax[0],kSMin[1],kSMin[2]);
        m_akSVertex[m_aiSV[2]] = Vector3<Real>(kSMin[0],kSMax[1],kSMin[2]);
        m_akSVertex[m_aiSV[3]] = Vector3<Real>(kSMin[0],kSMin[1],kSMax[2]);

        m_pkQuery = new Query3TRational<Real>(m_iVertexQuantity,m_akSVertex);
    }

    DelTetrahedron<Real>* pkTetra = new DelTetrahedron<Real>(m_aiSV[0],
        m_aiSV[1],m_aiSV[2],m_aiSV[3]);
    m_kTetrahedron.insert(pkTetra);

    // Incrementally update the tetrahedralization.  The set of processed
    // points is maintained to eliminate duplicates, either in the original
    // input points or in the points obtained by snap rounding.
    std::set<Vector3<Real> > kProcessed;
    for (i = 0; i < m_iVertexQuantity-4; i++)
    {
        if ( kProcessed.find(m_akSVertex[i]) == kProcessed.end() )
        {
            Update(i);
            kProcessed.insert(m_akSVertex[i]);
        }
    }
    m_iUniqueVertexQuantity = (int)kProcessed.size();

    // Remove tetrahedra sharing a vertex of the supertetrahedron.
    RemoveTetrahedra();

    // Assign integer values to the tetrahedra for use by the caller.
    std::map<DelTetrahedron<Real>*,int> kPermute;
    typename std::set<DelTetrahedron<Real>*>::iterator pkTIter =
        m_kTetrahedron.begin();
    for (i = 0; pkTIter != m_kTetrahedron.end(); pkTIter++)
    {
        pkTetra = *pkTIter;
        kPermute[pkTetra] = i++;
    }
    kPermute[0] = -1;

    // Put Delaunay tetrahedra into an array (vertices and adjacency info).
    m_iSimplexQuantity = (int)m_kTetrahedron.size();
    if ( m_iSimplexQuantity > 0 )
    {
        m_aiIndex = new int[4*m_iSimplexQuantity];
        m_aiAdjacent = new int[4*m_iSimplexQuantity];
        i = 0;
        pkTIter = m_kTetrahedron.begin();
        for (/**/; pkTIter != m_kTetrahedron.end(); pkTIter++)
        {
            pkTetra = *pkTIter;
            m_aiIndex[i] = pkTetra->V[0];
            m_aiAdjacent[i++] = kPermute[pkTetra->A[0]];
            m_aiIndex[i] = pkTetra->V[1];
            m_aiAdjacent[i++] = kPermute[pkTetra->A[1]];
            m_aiIndex[i] = pkTetra->V[2];
            m_aiAdjacent[i++] = kPermute[pkTetra->A[2]];
            m_aiIndex[i] = pkTetra->V[3];
            m_aiAdjacent[i++] = kPermute[pkTetra->A[3]];
        }
        assert( i == 4*m_iSimplexQuantity );

        m_iPathLast = -1;
        m_aiPath = new int[m_iSimplexQuantity+1];
    }

    // Restore the vertex count to the original (discards the vertices of the
    // supertetrahedron).
    m_iVertexQuantity -= 4;
}
//----------------------------------------------------------------------------
template <class Real>
Delaunay3<Real>::~Delaunay3 ()
{
    delete m_pkQuery;
    delete[] m_akSVertex;
    delete[] m_aiPath;
    if ( m_bOwner )
        delete[] m_akVertex;

    typename std::set<DelTetrahedron<Real>*>::iterator pkTIter =
        m_kTetrahedron.begin();
    for (/**/; pkTIter != m_kTetrahedron.end(); ++pkTIter)
        delete *pkTIter;
}
//----------------------------------------------------------------------------
template <class Real>
const Vector3<Real>* Delaunay3<Real>::GetVertices () const
{
    return m_akVertex;
}
//----------------------------------------------------------------------------
template <class Real>
int Delaunay3<Real>::GetUniqueVertexQuantity () const
{
    return m_iUniqueVertexQuantity;
}
//----------------------------------------------------------------------------
template <class Real>
const Vector3<Real>& Delaunay3<Real>::GetLineOrigin () const
{
    return m_kLineOrigin;
}
//----------------------------------------------------------------------------
template <class Real>
const Vector3<Real>& Delaunay3<Real>::GetLineDirection () const
{
    return m_kLineDirection;
}
//----------------------------------------------------------------------------
template <class Real>
Delaunay1<Real>* Delaunay3<Real>::GetDelaunay1 () const
{
    assert( m_iDimension == 1 );
    if ( m_iDimension != 1 )
        return 0;

    Real* afProjection = new Real[m_iVertexQuantity];
    for (int i = 0; i < m_iVertexQuantity; i++)
    {
        Vector3<Real> kDiff = m_akVertex[i] - m_kLineOrigin;
        afProjection[i] = m_kLineDirection.Dot(kDiff);
    }

    return new Delaunay1<Real>(m_iVertexQuantity,afProjection,m_fEpsilon,
        true,m_eComputationType);
}
//----------------------------------------------------------------------------
template <class Real>
const Vector3<Real>& Delaunay3<Real>::GetPlaneOrigin () const
{
    return m_kPlaneOrigin;
}
//----------------------------------------------------------------------------
template <class Real>
const Vector3<Real>& Delaunay3<Real>::GetPlaneDirection (int i) const
{
    assert( 0 <= i && i < 2 );
    return m_akPlaneDirection[i];
}
//----------------------------------------------------------------------------
template <class Real>
Delaunay2<Real>* Delaunay3<Real>::GetDelaunay2 () const
{
    assert( m_iDimension == 2 );
    if ( m_iDimension != 2 )
        return 0;

    Vector2<Real>* akProjection = new Vector2<Real>[m_iVertexQuantity];
    for (int i = 0; i < m_iVertexQuantity; i++)
    {
        Vector3<Real> kDiff = m_akVertex[i] - m_kPlaneOrigin;
        akProjection[i][0] = m_akPlaneDirection[0].Dot(kDiff);
        akProjection[i][1] = m_akPlaneDirection[1].Dot(kDiff);
    }

    return new Delaunay2<Real>(m_iVertexQuantity,akProjection,m_fEpsilon,
        true,m_eComputationType);
}
//----------------------------------------------------------------------------
template <class Real>
bool Delaunay3<Real>::GetHull (int& riTQuantity, int*& raiIndex) const
{
    assert( m_iDimension == 3 );
    if ( m_iDimension != 3 )
        return false;

    riTQuantity = 0;
    raiIndex = 0;

    // Count the number of triangles that are not shared by two tetrahedra.
    int i, iAdjQuantity = 4*m_iSimplexQuantity;
    for (i = 0; i < iAdjQuantity; i++)
    {
        if ( m_aiAdjacent[i] == -1 )
            riTQuantity++;
    }
    assert( riTQuantity > 0 );
    if ( riTQuantity == 0 )
        return false;

    // Enumerate the triangles.
    raiIndex = new int[3*riTQuantity];
    int* piIndex = raiIndex;
    for (i = 0; i < iAdjQuantity; i++)
    {
        if ( m_aiAdjacent[i] == -1 )
        {
            int iTetra = i/4, iFace = i%4;
            for (int j = 0; j < 4; j++)
            {
                if ( j != iFace )
                    *piIndex++ = m_aiIndex[4*iTetra+j];
            }
            if ( (iFace % 2) == 0 )
            {
                int iSave = *(piIndex-1);
                *(piIndex-1) = *(piIndex-2);
                *(piIndex-2) = iSave;
            }
        }
    }

    return true;
}
//----------------------------------------------------------------------------
template <class Real>
int Delaunay3<Real>::GetContainingTetrahedron (const Vector3<Real>& rkP)
    const
{
    assert(m_iDimension == 3);
    if (m_iDimension != 3)
    {
        return -1;
    }

    // convert to scaled coordinates
    Vector3<Real> kXFrmP = (rkP - m_kMin)*m_fScale;

    // start at first tetrahedron in mesh
    int iIndex = (m_iPathLast >= 0 ? m_aiPath[m_iPathLast] : 0);
    m_iPathLast = -1;
    m_iLastFaceV0 = -1;
    m_iLastFaceV1 = -1;
    m_iLastFaceV2 = -1;
    m_iLastFaceOpposite = -1;
    m_iLastFaceOppositeIndex = -1;

    // use tetrahedron faces as binary separating planes
    for (int i = 0; i < m_iSimplexQuantity; i++)
    {
        m_aiPath[++m_iPathLast] = iIndex;

        int* aiV = &m_aiIndex[4*iIndex];

        // <V1,V2,V3> counterclockwise when viewed outside tetrahedron
        if (m_pkQuery->ToPlane(kXFrmP,aiV[1],aiV[2],aiV[3]) > 0)
        {
            iIndex = m_aiAdjacent[4*iIndex];
            if (iIndex == -1)
            {
                m_iLastFaceV0 = aiV[1];
                m_iLastFaceV1 = aiV[2];
                m_iLastFaceV2 = aiV[3];
                m_iLastFaceOpposite = aiV[0];
                m_iLastFaceOppositeIndex = 0;
                return -1;
            }
            continue;
        }

        // <V0,V3,V2> counterclockwise when viewed outside tetrahedron
        if (m_pkQuery->ToPlane(kXFrmP,aiV[0],aiV[2],aiV[3]) < 0)
        {
            iIndex = m_aiAdjacent[4*iIndex+1];
            if (iIndex == -1)
            {
                m_iLastFaceV0 = aiV[0];
                m_iLastFaceV1 = aiV[2];
                m_iLastFaceV2 = aiV[3];
                m_iLastFaceOpposite = aiV[1];
                m_iLastFaceOppositeIndex = 1;
                return -1;
            }
            continue;
        }

        // <V0,V1,V3> counterclockwise when viewed outside tetrahedron
        if (m_pkQuery->ToPlane(kXFrmP,aiV[0],aiV[1],aiV[3]) > 0)
        {
            iIndex = m_aiAdjacent[4*iIndex+2];
            if (iIndex == -1)
            {
                m_iLastFaceV0 = aiV[0];
                m_iLastFaceV1 = aiV[1];
                m_iLastFaceV2 = aiV[3];
                m_iLastFaceOpposite = aiV[2];
                m_iLastFaceOppositeIndex = 2;
                return -1;
            }
            continue;
        }

        // <V0,V2,V1> counterclockwise when viewed outside tetrahedron
        if (m_pkQuery->ToPlane(kXFrmP,aiV[0],aiV[1],aiV[2]) < 0)
        {
            iIndex = m_aiAdjacent[4*iIndex+3];
            if (iIndex == -1)
            {
                m_iLastFaceV0 = aiV[0];
                m_iLastFaceV1 = aiV[1];
                m_iLastFaceV2 = aiV[2];
                m_iLastFaceOpposite = aiV[3];
                m_iLastFaceOppositeIndex = 3;
                return -1;
            }
            continue;
        }

        m_iLastFaceV0 = -1;
        m_iLastFaceV1 = -1;
        m_iLastFaceV2 = -1;
        m_iLastFaceOppositeIndex = -1;
        return iIndex;
    }

    return -1;
}
//----------------------------------------------------------------------------
template <class Real>
int Delaunay3<Real>::GetPathLast () const
{
    return m_iPathLast;
}
//----------------------------------------------------------------------------
template <class Real>
const int* Delaunay3<Real>::GetPath () const
{
    return m_aiPath;
}
//----------------------------------------------------------------------------
template <class Real>
int Delaunay3<Real>::GetLastFace (int& riV0, int& riV1, int& riV2,
    int& riV3) const
{
    riV0 = m_iLastFaceV0;
    riV1 = m_iLastFaceV1;
    riV2 = m_iLastFaceV2;
    riV3 = m_iLastFaceOpposite;
    return m_iLastFaceOppositeIndex;
}
//----------------------------------------------------------------------------
template <class Real>
bool Delaunay3<Real>::GetVertexSet (int i, Vector3<Real> akV[4]) const
{
    assert( m_iDimension == 3 );
    if ( m_iDimension != 3 )
        return false;

    if ( 0 <= i && i < m_iSimplexQuantity )
    {
        akV[0] = m_akVertex[m_aiIndex[4*i  ]];
        akV[1] = m_akVertex[m_aiIndex[4*i+1]];
        akV[2] = m_akVertex[m_aiIndex[4*i+2]];
        akV[3] = m_akVertex[m_aiIndex[4*i+3]];
        return true;
    }

    return false;
}
//----------------------------------------------------------------------------
template <class Real>
bool Delaunay3<Real>::GetIndexSet (int i, int aiIndex[4]) const
{
    assert( m_iDimension == 3 );
    if ( m_iDimension != 3 )
        return false;

    if ( 0 <= i && i < m_iSimplexQuantity )
    {
        aiIndex[0] = m_aiIndex[4*i  ];
        aiIndex[1] = m_aiIndex[4*i+1];
        aiIndex[2] = m_aiIndex[4*i+2];
        aiIndex[3] = m_aiIndex[4*i+3];
        return true;
    }

    return false;
}
//----------------------------------------------------------------------------
template <class Real>
bool Delaunay3<Real>::GetAdjacentSet (int i, int aiAdjacent[4]) const
{
    assert( m_iDimension == 3 );
    if ( m_iDimension != 3 )
        return false;

    if ( 0 <= i && i < m_iSimplexQuantity )
    {
        aiAdjacent[0] = m_aiAdjacent[4*i  ];
        aiAdjacent[1] = m_aiAdjacent[4*i+1];
        aiAdjacent[2] = m_aiAdjacent[4*i+2];
        aiAdjacent[3] = m_aiAdjacent[4*i+3];
        return true;
    }

    return false;
}
//----------------------------------------------------------------------------
template <class Real>
bool Delaunay3<Real>::GetBarycentricSet (int i, const Vector3<Real>& rkP,
    Real afBary[4]) const
{
    assert( m_iDimension == 3 );
    if ( m_iDimension != 3 )
        return false;

    if ( 0 <= i && i < m_iSimplexQuantity )
    {
        Vector3<Real> kV0 = m_akVertex[m_aiIndex[4*i  ]];
        Vector3<Real> kV1 = m_akVertex[m_aiIndex[4*i+1]];
        Vector3<Real> kV2 = m_akVertex[m_aiIndex[4*i+2]];
        Vector3<Real> kV3 = m_akVertex[m_aiIndex[4*i+3]];
        rkP.GetBarycentrics(kV0,kV1,kV2,kV3,afBary);
        return true;
    }

    return false;
}
//----------------------------------------------------------------------------
template <class Real>
void Delaunay3<Real>::Update (int i)
{
    // Locate the tetrahedron containing vertex i.
    DelTetrahedron<Real>* pkTetra = GetContainingTetrahedron(i);

    // Locate and remove the tetrahedra forming the insertion polyhedron.
    std::stack<DelTetrahedron<Real>*> kStack;
    ETManifoldMesh kPolyhedron(0,DelPolyhedronFace<Real>::TCreator);
    kStack.push(pkTetra);
    pkTetra->OnStack = true;
    int j, iV0, iV1, iV2;
    DelPolyhedronFace<Real>* pkFace;
    while ( !kStack.empty() )
    {
        pkTetra = kStack.top();
        kStack.pop();
        pkTetra->OnStack = false;
        for (j = 0; j < 4; j++)
        {
            DelTetrahedron<Real>* pkAdj = pkTetra->A[j];
            if ( pkAdj )
            {
                // Detach tetrahedron and adjacent tetrahedron from each
                // other.
                int iNullIndex = pkTetra->DetachFrom(j,pkAdj);

                if ( pkAdj->IsInsertionComponent(i,pkTetra,m_pkQuery,m_aiSV) )
                {
                    if ( !pkAdj->OnStack )
                    {
                        // Adjacent triangle inside insertion polyhedron.
                        kStack.push(pkAdj);
                        pkAdj->OnStack = true;
                    }
                }
                else
                {
                    // Adjacent tetrahedron outside insertion polyhedron.
                    iV0 = pkTetra->V[gs_aaiIndex[j][0]];
                    iV1 = pkTetra->V[gs_aaiIndex[j][1]];
                    iV2 = pkTetra->V[gs_aaiIndex[j][2]];
                    pkFace = (DelPolyhedronFace<Real>*)
                        kPolyhedron.InsertTriangle(iV0,iV1,iV2);
                    pkFace->NullIndex = iNullIndex;
                    pkFace->Tetra = pkAdj;
                }
            }
            else
            {
                // The tetrahedron is in the insertion polyhedron, but the
                // adjacent one does not exist.  This means one of two things:
                // (1) We are at a face of the supertetrahedron, and that
                //     face is part of the insertion polyhedron.
                // (2) We are at a face that was recently shared by the
                //     tetrahedron and the adjacent, but we detached those
                //     tetrahedra from each other.  These faces should be
                //     ignored.

                iV0 = pkTetra->V[gs_aaiIndex[j][0]];
                if ( IsSupervertex(iV0) )
                {
                    iV1 = pkTetra->V[gs_aaiIndex[j][1]];
                    if ( IsSupervertex(iV1) )
                    {
                        iV2 = pkTetra->V[gs_aaiIndex[j][2]];
                        if ( IsSupervertex(iV2) )
                        {
                            pkFace = (DelPolyhedronFace<Real>*)
                                kPolyhedron.InsertTriangle(iV0,iV1,iV2);
                            pkFace->NullIndex = -1;
                            pkFace->Tetra = 0;
                        }
                    }
                }
            }
        }
        m_kTetrahedron.erase(pkTetra);
        delete pkTetra;
    }

    // Insert the new tetrahedra formed by the input point and the faces of
    // the insertion polyhedron.
    const ETManifoldMesh::TMap& rkTMap = kPolyhedron.GetTriangles();
    assert( rkTMap.size() >= 4 && kPolyhedron.IsClosed() );
    typename ETManifoldMesh::TMapCIterator pkTIter;
    for (pkTIter = rkTMap.begin(); pkTIter != rkTMap.end(); pkTIter++)
    {
        pkFace = (DelPolyhedronFace<Real>*)pkTIter->second;

        // Create and insert the new tetrahedron.
        pkTetra = new DelTetrahedron<Real>(i,pkFace->V[0],pkFace->V[1],
            pkFace->V[2]);
        m_kTetrahedron.insert(pkTetra);

        // Establish the adjacency links across the polyhedron face.
        pkTetra->A[0] = pkFace->Tetra;
        if ( pkFace->Tetra )
            pkFace->Tetra->A[pkFace->NullIndex] = pkTetra;

        // Update the faces's tetrahedron pointer to point to the newly
        // created tetrahedron.  This information is used later to establish
        // the links between the new tetrahedra.
        pkFace->Tetra = pkTetra;
    }

    // Establish the adjacency links between the new tetrahedra.
    DelPolyhedronFace<Real>* pkAdjFace;
    for (pkTIter = rkTMap.begin(); pkTIter != rkTMap.end(); pkTIter++)
    {
        pkFace = (DelPolyhedronFace<Real>*)pkTIter->second;

        pkAdjFace = (DelPolyhedronFace<Real>*)pkFace->T[0];
        pkFace->Tetra->A[3] = pkAdjFace->Tetra;
        assert( SharesFace(3,pkFace->Tetra,pkAdjFace->Tetra) );

        pkAdjFace = (DelPolyhedronFace<Real>*)pkFace->T[1];
        pkFace->Tetra->A[1] = pkAdjFace->Tetra;
        assert( SharesFace(1,pkFace->Tetra,pkAdjFace->Tetra) );

        pkAdjFace = (DelPolyhedronFace<Real>*)pkFace->T[2];
        pkFace->Tetra->A[2] = pkAdjFace->Tetra;
        assert( SharesFace(2,pkFace->Tetra,pkAdjFace->Tetra) );
    }
}
//----------------------------------------------------------------------------
template <class Real>
DelTetrahedron<Real>* Delaunay3<Real>::GetContainingTetrahedron (int i) const
{
    // Locate which tetrahedron in the current mesh contains vertex i.  By
    // construction, there must be such a tetrahedron (the vertex cannot be
    // outside the supertetrahedron).

    DelTetrahedron<Real>* pkTetra = *m_kTetrahedron.begin();
    int iTQuantity = (int)m_kTetrahedron.size();
    for (int iT = 0; iT < iTQuantity; iT++)
    {
        int* aiV = pkTetra->V;

        // <V1,V2,V3> counterclockwise when viewed outside tetrahedron
        if ( m_pkQuery->ToPlane(i,aiV[1],aiV[2],aiV[3]) > 0 )
        {
            pkTetra = pkTetra->A[0];
            if ( !pkTetra )
                break;
            continue;
        }

        // <V0,V3,V2> counterclockwise when viewed outside tetrahedron
        if ( m_pkQuery->ToPlane(i,aiV[0],aiV[2],aiV[3]) < 0 )
        {
            pkTetra = pkTetra->A[1];
            if ( !pkTetra )
                break;
            continue;
        }

        // <V0,V1,V3> counterclockwise when viewed outside tetrahedron
        if ( m_pkQuery->ToPlane(i,aiV[0],aiV[1],aiV[3]) > 0 )
        {
            pkTetra = pkTetra->A[2];
            if ( !pkTetra )
                break;
            continue;
        }

        // <V0,V2,V1> counterclockwise when viewed outside tetrahedron
        if ( m_pkQuery->ToPlane(i,aiV[0],aiV[1],aiV[2]) < 0 )
        {
            pkTetra = pkTetra->A[3];
            if ( !pkTetra )
                break;
            continue;
        }

        return pkTetra;
    }

    assert( false );
    return 0;
}
//----------------------------------------------------------------------------
template <class Real>
void Delaunay3<Real>::RemoveTetrahedra ()
{
    // Identify those triangles sharing a vertex of the supertetrahedron.
    std::set<DelTetrahedron<Real>*> kRemoveTetra;
    DelTetrahedron<Real>* pkTetra;
    typename std::set<DelTetrahedron<Real>*>::iterator pkTIter =
        m_kTetrahedron.begin();
    for (/**/; pkTIter != m_kTetrahedron.end(); pkTIter++)
    {
        pkTetra = *pkTIter;
        for (int j = 0; j < 4; j++)
        {
            if ( IsSupervertex(pkTetra->V[j]) )
            {
                kRemoveTetra.insert(pkTetra);
                break;
            }
        }
    }

    // Remove the tetrahedra from the mesh.
    pkTIter = kRemoveTetra.begin();
    for (/**/; pkTIter != kRemoveTetra.end(); pkTIter++)
    {
        pkTetra = *pkTIter;
        for (int j = 0; j < 4; j++)
        {
            // Break the links with adjacent tetrahedra.
            DelTetrahedron<Real>* pkAdj = pkTetra->A[j];
            if ( pkAdj )
            {
                for (int k = 0; k < 4; k++)
                {
                    if ( pkAdj->A[k] == pkTetra )
                    {
                        pkAdj->A[k] = 0;
                        break;
                    }
                }
            }
        }
        m_kTetrahedron.erase(pkTetra);
        delete pkTetra;
    }
}
//----------------------------------------------------------------------------
template <class Real>
bool Delaunay3<Real>::IsSupervertex (int i) const
{
    for (int j = 0; j < 4; j++)
    {
        if ( i == m_aiSV[j] )
            return true;
    }
    return false;
}
//----------------------------------------------------------------------------
template <class Real>
bool Delaunay3<Real>::SharesFace (int i, DelTetrahedron<Real>* pkFace,
    DelTetrahedron<Real>* pkAdj)
{
    int aiF[3], iCount = 0, j;
    for (j = 0; j < 4; j++)
    {
        if ( j != i )
            aiF[iCount++] = pkFace->V[j];
    }

    for (i = 0; i < 4; i++)
    {
        if ( pkAdj->V[i] != aiF[0] &&
             pkAdj->V[i] != aiF[1] &&
             pkAdj->V[i] != aiF[2] )
        {
            break;
        }
    }
    if ( i == 4 )
        return false;

    int aiA[3];
    for (j = 0, iCount = 0; j < 4; j++)
    {
        if ( j != i )
            aiA[iCount++] = pkAdj->V[j];
    }

    if ( aiF[0] > aiF[1] )
    {
        j = aiF[0];
        aiF[0] = aiF[1];
        aiF[1] = j;
    }
    if ( aiF[1] > aiF[2] )
    {
        j = aiF[1];
        aiF[1] = aiF[2];
        aiF[2] = j;
    }
    if ( aiF[0] > aiF[1] )
    {
        j = aiF[0];
        aiF[0] = aiF[1];
        aiF[1] = j;
    }

    if ( aiA[0] > aiA[1] )
    {
        j = aiA[0];
        aiA[0] = aiA[1];
        aiA[1] = j;
    }
    if ( aiA[1] > aiA[2] )
    {
        j = aiA[1];
        aiA[1] = aiA[2];
        aiA[2] = j;
    }
    if ( aiA[0] > aiA[1] )
    {
        j = aiA[0];
        aiA[0] = aiA[1];
        aiA[1] = j;
    }

    if ( aiA[0] != aiF[0] || aiA[1] != aiF[1] || aiA[2] != aiF[2] )
        return false;

    return true;
}
//----------------------------------------------------------------------------

//----------------------------------------------------------------------------
// explicit instantiation
//----------------------------------------------------------------------------
template WM3_ITEM class Delaunay3<float>;
template WM3_ITEM class Delaunay3<double>;
//----------------------------------------------------------------------------
}


