wmlextractsurfacetetra.cpp

来自「3D Game Engine Design Source Code非常棒」· C++ 代码 · 共 1,297 行 · 第 1/3 页

// Magic Software, Inc.
// http://www.magic-software.com
// http://www.wild-magic.com
// Copyright (c) 2003.  All Rights Reserved
//
// The Wild Magic Library (WML) source code is supplied under the terms of
// the license agreement http://www.magic-software.com/License/WildMagic.pdf
// and may not be copied or disclosed except in accordance with the terms of
// that agreement.

#include "WmlExtractSurfaceTetra.h"
#include "WmlVETMesh.h"
using namespace Wml;

#include <algorithm>
using namespace std;

typedef map<Vector3f,int> VMap;
typedef VMap::iterator VMapIterator;
typedef map<TriangleKey,int> TMap;
typedef TMap::iterator TMapIterator;

//----------------------------------------------------------------------------
ExtractSurfaceTetra::ExtractSurfaceTetra (int iXBound, int iYBound,
    int iZBound, int* aiData)
{
    assert( iXBound > 0 && iYBound > 0 && iZBound > 0 && aiData );
    m_iXBound = iXBound;
    m_iYBound = iYBound;
    m_iZBound = iZBound;
    m_iXYBound = iXBound*iYBound;
    m_iXYZBound = m_iXYBound*iZBound;
    m_aiData = aiData;
}
//----------------------------------------------------------------------------
float ExtractSurfaceTetra::GetFunction (const Vector3f& rkP) const
{
    int iX = (int) rkP.X();
    if ( iX < 0 || iX >= m_iXBound-1 )
        return 0.0f;

    int iY = (int) rkP.Y();
    if ( iY < 0 || iY >= m_iYBound-1 )
        return 0.0f;

    int iZ = (int) rkP.Z();
    if ( iZ < 0 || iZ >= m_iZBound-1 )
        return 0.0f;

    float fDX = rkP.X() - iX, fDY = rkP.Y() - iY, fDZ = rkP.Z() - iZ;

    int i000 = iX + m_iXBound*(iY + m_iYBound*iZ);
    int i100 = i000 + 1;
    int i010 = i000 + m_iXBound;
    int i110 = i100 + m_iXBound;
    int i001 = i000 + m_iXYBound;
    int i101 = i100 + m_iXYBound;
    int i011 = i010 + m_iXYBound;
    int i111 = i110 + m_iXYBound;
    float fF000 = (float) m_aiData[i000];
    float fF100 = (float) m_aiData[i100];
    float fF010 = (float) m_aiData[i010];
    float fF110 = (float) m_aiData[i110];
    float fF001 = (float) m_aiData[i001];
    float fF101 = (float) m_aiData[i101];
    float fF011 = (float) m_aiData[i011];
    float fF111 = (float) m_aiData[i111];
    float fC0, fC1, fC2, fInterp;

    if ( (iX & 1) ^ (iY & 1) ^ (iZ & 1) )
    {
        if ( fDX - fDY - fDZ >= 0.0f )
        {
            // 1205
            fInterp =
                (1.0f-(1.0f-fDX)-fDY-fDZ)*fF100 +
                (1.0f-fDX)*fF000 +
                fDY*fF110 +
                fDZ*fF101;
        }
        else if ( fDX - fDY + fDZ <= 0.0f )
        {
            // 3027
            fInterp =
                (1.0f-fDX-(1.0f-fDY)-fDZ)*fF010 +
                fDX*fF110 +
                (1.0f-fDY)*fF000 +
                fDZ*fF011;
        }
        else if ( fDX + fDY - fDZ <= 0.0f )
        {
            // 4750
            fInterp =
                (1.0f-fDX-fDY-(1-fDZ))*fF001 +
                fDX*fF101 +
                fDY*fF011 +
                (1.0f-fDZ)*fF000;
        }
        else if ( fDX + fDY + fDZ >= 2.0f )
        {
            // 6572
            fInterp =
                (1.0f-(1.0f-fDX)-(1.0f-fDY)-(1.0f-fDZ))*fF111 +
                (1.0f-fDX)*fF011 +
                (1.0f-fDY)*fF101 +
                (1.0f-fDZ)*fF110;
        }
        else
        {
            // 0752
            fC0 = 0.5f*(-fDX+fDY+fDZ);
            fC1 = 0.5f*(fDX-fDY+fDZ);
            fC2 = 0.5f*(fDX+fDY-fDZ);
            fInterp =
                (1.0f-fC0-fC1-fC2)*fF000 +
                fC0*fF011 +
                fC1*fF101 +
                fC2*fF110;
        }
    }
    else
    {
        if ( fDX + fDY + fDZ <= 1.0f )
        {
            // 0134
            fInterp =
                (1.0f-fDX-fDY-fDZ)*fF000 +
                fDX*fF100 +
                fDY*fF010 +
                fDZ*fF001;
        }
        else if ( fDX + fDY - fDZ >= 1.0f )
        {
            // 2316
            fInterp =
                (1.0f-(1.0f-fDX)-(1.0f-fDY)-fDZ)*fF110 +
                (1.0f-fDX)*fF010 +
                (1.0f-fDY)*fF100 +
                fDZ*fF111;
        }
        else if ( fDX - fDY + fDZ >= 1.0f )
        {
            // 5461
            fInterp =
                (1.0f-(1.0f-fDX)-fDY-(1.0f-fDZ))*fF101 +
                (1.0f-fDX)*fF001 +
                fDY*fF111 +
                (1.0f-fDZ)*fF100;
        }
        else if ( -fDX + fDY + fDZ >= 1.0f )
        {
            // 7643
            fInterp =
                (1.0f-fDX-(1.0f-fDY)-(1.0f-fDZ))*fF011 +
                fDX*fF111 +
                (1.0f-fDY)*fF001 +
                (1.0f-fDZ)*fF010;
        }
        else
        {
            // 6314
            fC0 = 0.5f*((1.0f-fDX)-(1.0f-fDY)+(1.0f-fDZ));
            fC1 = 0.5f*(-(1.0f-fDX)+(1.0f-fDY)+(1.0f-fDZ));
            fC2 = 0.5f*((1.0f-fDX)+(1.0f-fDY)-(1.0f-fDZ));
            fInterp =
                (1.0f-fC0-fC1-fC2)*fF111 +
                fC0*fF010 +
                fC1*fF100 +
                fC2*fF001;
        }
    }

    return fInterp;
}
//----------------------------------------------------------------------------
Vector3f ExtractSurfaceTetra::GetGradient (const Vector3f& rkP) const
{
    int iX = (int) rkP.X();
    if ( iX < 0 || iX >= m_iXBound-1 )
        return Vector3f::ZERO;

    int iY = (int) rkP.Y();
    if ( iY < 0 || iY >= m_iYBound-1 )
        return Vector3f::ZERO;

    int iZ = (int) rkP.Z();
    if ( iZ < 0 || iZ >= m_iZBound-1 )
        return Vector3f::ZERO;

    float fDX = rkP.X() - iX, fDY = rkP.Y() - iY, fDZ = rkP.Z() - iZ;

    int i000 = iX + m_iXBound*(iY + m_iYBound*iZ);
    int i100 = i000 + 1;
    int i010 = i000 + m_iXBound;
    int i110 = i100 + m_iXBound;
    int i001 = i000 + m_iXYBound;
    int i101 = i100 + m_iXYBound;
    int i011 = i010 + m_iXYBound;
    int i111 = i110 + m_iXYBound;
    float fF000 = (float) m_aiData[i000];
    float fF100 = (float) m_aiData[i100];
    float fF010 = (float) m_aiData[i010];
    float fF110 = (float) m_aiData[i110];
    float fF001 = (float) m_aiData[i001];
    float fF101 = (float) m_aiData[i101];
    float fF011 = (float) m_aiData[i011];
    float fF111 = (float) m_aiData[i111];
    Vector3f kInterp;

    if ( (iX & 1) ^ (iY & 1) ^ (iZ & 1) )
    {
        if ( fDX - fDY - fDZ >= 0.0f )
        {
            // 1205
            kInterp.X() = + fF100 - fF000;
            kInterp.Y() = - fF100 + fF110;
            kInterp.Z() = - fF100 + fF101; 
        }
        else if ( fDX - fDY + fDZ <= 0.0f )
        {
            // 3027
            kInterp.X() = - fF010 + fF110;
            kInterp.Y() = + fF010 - fF000;
            kInterp.Z() = - fF010 + fF011;
        }
        else if ( fDX + fDY - fDZ <= 0.0f )
        {
            // 4750
            kInterp.X() = - fF001 + fF101;
            kInterp.Y() = - fF001 + fF011;
            kInterp.Z() = + fF001 - fF000;
        }
        else if ( fDX + fDY + fDZ >= 2.0f )
        {
            // 6572
            kInterp.X() = + fF111 - fF011;
            kInterp.Y() = + fF111 - fF101;
            kInterp.Z() = + fF111 - fF110;
        }
        else
        {
            // 0752
            kInterp.X() = 0.5f*(-fF000-fF011+fF101+fF110);
            kInterp.Y() = 0.5f*(-fF000+fF011-fF101+fF110);
            kInterp.Z() = 0.5f*(-fF000+fF011+fF101-fF110);
        }
    }
    else
    {
        if ( fDX + fDY + fDZ <= 1.0f )
        {
            // 0134
            kInterp.X() = - fF000 + fF100;
            kInterp.Y() = - fF000 + fF010;
            kInterp.Z() = - fF000 + fF001;
        }
        else if ( fDX + fDY - fDZ >= 1.0f )
        {
            // 2316
            kInterp.X() = + fF110 - fF010;
            kInterp.Y() = + fF110 - fF100;
            kInterp.Z() = - fF110 + fF111;
        }
        else if ( fDX - fDY + fDZ >= 1.0f )
        {
            // 5461
            kInterp.X() = + fF101 - fF001;
            kInterp.Y() = - fF101 + fF111;
            kInterp.Z() = + fF101 - fF100;
        }
        else if ( -fDX + fDY + fDZ >= 1.0f )
        {
            // 7643
            kInterp.X() = - fF011 + fF111;
            kInterp.Y() = + fF011 - fF001;
            kInterp.Z() = + fF011 - fF010;
        }
        else
        {
            // 6314
            kInterp.X() = 0.5f*(fF111-fF010+fF100-fF001);
            kInterp.Y() = 0.5f*(fF111+fF010-fF100-fF001);
            kInterp.Z() = 0.5f*(fF111-fF010-fF100+fF001);
        }
    }

    return kInterp;
}
//----------------------------------------------------------------------------
void ExtractSurfaceTetra::ExtractContour (int iLevel, vector<Vector3f>& rkVA,
    vector<TriangleKey>& rkTA)
{
    // NOTE:  Isolated edges are computed, but not reported to the caller.
    // You can modify the code to return these if so desired.
    VtxMap kVMap;
    ESet kESet;
    TSet kTSet;
    m_iNextIndex = 0;

    // adjust image so level set is F(x,y,z) = 0
    int i;
    for (i = 0; i < m_iXYZBound; i++)
        m_aiData[i] = m_aiData[i] - iLevel;

    int iXBoundM1 = m_iXBound - 1;
    int iYBoundM1 = m_iYBound - 1;
    int iZBoundM1 = m_iZBound - 1;
    for (int iZ = 0, iZP = 1; iZ < iZBoundM1; iZ++, iZP++)
    {
        int iZParity = (iZ & 1);

        for (int iY = 0, iYP = 1; iY < iYBoundM1; iY++, iYP++)
        {
            int iYParity = (iY & 1);

            for (int iX = 0, iXP = 1; iX < iXBoundM1; iX++, iXP++)
            {
                int iXParity = (iX & 1);

                int i000 = iX + m_iXBound*(iY + m_iYBound*iZ);
                int i100 = i000 + 1;
                int i010 = i000 + m_iXBound;
                int i110 = i100 + m_iXBound;
                int i001 = i000 + m_iXYBound;
                int i101 = i100 + m_iXYBound;
                int i011 = i010 + m_iXYBound;
                int i111 = i110 + m_iXYBound;
                int iF000 = m_aiData[i000];
                int iF100 = m_aiData[i100];
                int iF010 = m_aiData[i010];
                int iF110 = m_aiData[i110];
                int iF001 = m_aiData[i001];
                int iF101 = m_aiData[i101];
                int iF011 = m_aiData[i011];
                int iF111 = m_aiData[i111];

                if ( iXParity ^ iYParity ^ iZParity )
                {
                    // 1205
                    ProcessTetrahedron(kVMap,kESet,kTSet,
                        iXP,iY,iZ,iF100,iXP,iYP,iZ,iF110,iX,iY,iZ,iF000,iXP,
                        iY,iZP,iF101);

                    // 3027
                    ProcessTetrahedron(kVMap,kESet,kTSet,
                        iX,iYP,iZ,iF010,iX,iY,iZ,iF000,iXP,iYP,iZ,iF110,iX,
                        iYP,iZP,iF011);

                    // 4750
                    ProcessTetrahedron(kVMap,kESet,kTSet,
                        iX,iY,iZP,iF001,iX,iYP,iZP,iF011,iXP,iY,iZP,iF101,iX,
                        iY,iZ,iF000);

                    // 6572
                    ProcessTetrahedron(kVMap,kESet,kTSet,
                        iXP,iYP,iZP,iF111,iXP,iY,iZP,iF101,iX,iYP,iZP,iF011,
                        iXP,iYP,iZ,iF110);

                    // 0752
                    ProcessTetrahedron(kVMap,kESet,kTSet,
                        iX,iY,iZ,iF000,iX,iYP,iZP,iF011,iXP,iY,iZP,iF101,iXP,
                        iYP,iZ,iF110);
                }
                else
                {
                    // 0134
                    ProcessTetrahedron(kVMap,kESet,kTSet,
                        iX,iY,iZ,iF000,iXP,iY,iZ,iF100,iX,iYP,iZ,iF010,iX,iY,
                        iZP,iF001);

                    // 2316
                    ProcessTetrahedron(kVMap,kESet,kTSet,
                        iXP,iYP,iZ,iF110,iX,iYP,iZ,iF010,iXP,iY,iZ,iF100,iXP,
                        iYP,iZP,iF111);

                    // 5461
                    ProcessTetrahedron(kVMap,kESet,kTSet,
                        iXP,iY,iZP,iF101,iX,iY,iZP,iF001,iXP,iYP,iZP,iF111,
                        iXP,iY,iZ,iF100);

                    // 7643
                    ProcessTetrahedron(kVMap,kESet,kTSet,
                        iX,iYP,iZP,iF011,iXP,iYP,iZP,iF111,iX,iY,iZP,iF001,iX,
                        iYP,iZ,iF010);

                    // 6314
                    ProcessTetrahedron(kVMap,kESet,kTSet,
                        iXP,iYP,iZP,iF111,iX,iYP,iZ,iF010,iXP,iY,iZ,iF100,iX,
                        iY,iZP,iF001);
                }
            }
        }
    }

    // readjust image so level set is F(x,y,z) = L
    for (i = 0; i < m_iXYZBound; i++)
        m_aiData[i] = m_aiData[i] + iLevel;

    // pack vertices into an array
    rkVA.resize(kVMap.size());
    if ( rkVA.size() > 0 )
    {
        VtxMapIterator pkVIter;
        for (pkVIter = kVMap.begin(); pkVIter != kVMap.end(); pkVIter++)
        {
            const Vertex& rkV = pkVIter->first;
            rkV.GetTriple(rkVA[pkVIter->second]);
        }
    }

    // pack edges into an array (computed, but not reported to caller)
    vector<EdgeKey> kEA;
    kEA.resize(kESet.size());
    if ( kEA.size() > 0 )
        copy(kESet.begin(),kESet.end(),&kEA.front());

    // pack triangles into an array
    rkTA.resize(kTSet.size());
    if ( rkTA.size() > 0 )
        copy(kTSet.begin(),kTSet.end(),&rkTA.front());
}
//----------------------------------------------------------------------------
void ExtractSurfaceTetra::MakeUnique (vector<Vector3f>& rkVA,
    vector<TriangleKey>& rkTA)
{
    int iVQuantity = (int)rkVA.size();
    int iTQuantity = (int)rkTA.size();
    if ( iVQuantity == 0 || iTQuantity == 0 )
        return;

    // use a hash table to generate unique storage
    VMap kVMap;
    VMapIterator pkVIter;