mgcquadricsurface.cpp

《3D游戏引擎设计》的源码

// Magic Software, Inc.
// http://www.magic-software.com
// Copyright (c) 2000, All Rights Reserved
//
// Source code from Magic Software is supplied under the terms of a license
// agreement and may not be copied or disclosed except in accordance with the
// terms of that agreement.  The various license agreements may be found at
// the Magic Software web site.  This file is subject to the license
//
// FREE SOURCE CODE
// http://www.magic-software.com/License/free.pdf

#include "MgcEigen.h"
#include "MgcQuadricSurface.h"
#include "MgcRTLib.h"

//----------------------------------------------------------------------------
MgcQuadricSurface::MgcQuadricSurface (const MgcMatrix3& rkA,
    const MgcVector3& rkB, MgcReal fC)
    :
    m_kA(rkA), m_kB(rkB)
{
    m_fC = fC;
}
//----------------------------------------------------------------------------
void MgcQuadricSurface::GetCharacterization (Type& eType,
    MgcReal afD[3]) const
{
    eType = QST_MAX_TYPE;

    MgcEigen kES(3);
    int iRow, iCol;
    for (iRow = 0; iRow < 3; iRow++)
    {
        for (iCol = 0; iCol < 3; iCol++)
            kES.Matrix(iRow,iCol) = m_kA[iRow][iCol];
    }

    kES.IncrSortEigenStuff3();
    for (iRow = 0; iRow < 3; iRow++)
        afD[iRow] = kES.GetEigenvalue(iRow);

    const MgcReal fEpsilon = 1e-06;
    if ( MgcMath::Abs(afD[2]) > fEpsilon )
    {
        // matrix A has at least one nonzero eigenvalue, change coordinates
        MgcVector3 kE = MgcVector3::ZERO;
        for (iRow = 0; iRow < 3; iRow++)
        {
            for (iCol = 0; iCol < 3; iCol++)
                kE[iRow] += kES.Matrix(iRow,iCol)*m_kB[iCol];
        }

        MgcReal fF;
        if ( MgcMath::Abs(afD[1]) > fEpsilon )
        {
            if ( MgcMath::Abs(afD[0]) > fEpsilon )
            {
                // matrix A has three nonzero eigenvalues
                fF = -m_fC;
                for (iRow = 0; iRow <= 2; iRow++)
                    fF += 0.25*kE[iRow]*kE[iRow]/afD[iRow];

                if ( fF > fEpsilon )
                {
                    if ( afD[2] < 0.0 )
                        eType = QST_NONE;
                    else if ( afD[0] > 0.0 )
                        eType = QST_ELLIPSOID;
                    else if ( afD[1] > 0.0 )
                        eType = QST_HYPERBOLOID_ONE_SHEET;
                    else
                        eType = QST_HYPERBOLOID_TWO_SHEETS;
                }
                else if ( fF < -fEpsilon )
                {
                    if ( afD[0] > 0.0 )
                        eType = QST_NONE;
                    else if ( afD[2] < 0.0 )
                        eType = QST_ELLIPSOID;
                    else if ( afD[1] > 0.0 )
                        eType = QST_HYPERBOLOID_ONE_SHEET;
                    else
                        eType = QST_HYPERBOLOID_TWO_SHEETS;
                }
                else
                {
                    if ( afD[0] > 0.0 || afD[2] < 0.0 )
                        eType = QST_POINT;
                    else
                        eType = QST_ELLIPTIC_CONE;
                }
            }
            else
            {
                // matrix A has two nonzero eigenvalues (d1,d2)
                if ( MgcMath::Abs(kE[0]) > fEpsilon )
                {
                    fF = -m_fC;
                    for (iRow = 1; iRow <= 2; iRow++)
                        fF += 0.25*kE[iRow]*kE[iRow]/afD[iRow];

                    if ( fF > fEpsilon )
                    {
                        if ( afD[1] > 0.0 )
                        {
                            if ( afD[2] > 0.0 )
                                eType = QST_ELLIPTIC_CYLINDER;
                            else
                                eType = QST_HYPERBOLIC_CYLINDER;
                        }
                        else
                        {
                            if ( afD[2] > 0.0 )
                                eType = QST_HYPERBOLIC_CYLINDER;
                            else
                                eType = QST_NONE;
                        }
                    }
                    else if ( fF < -fEpsilon )
                    {
                        if ( afD[1] > 0.0 )
                        {
                            if ( afD[2] > 0.0 )
                                eType = QST_NONE;
                            else
                                eType = QST_HYPERBOLIC_CYLINDER;
                        }
                        else
                        {
                            if ( afD[2] > 0.0 )
                                eType = QST_HYPERBOLIC_CYLINDER;
                            else
                                eType = QST_ELLIPTIC_CYLINDER;
                        }
                    }
                    else
                    {
                        if ( afD[1]*afD[2] > 0.0 )
                            eType = QST_LINE;
                        else
                            eType = QST_TWO_PLANES;
                    }
                }
                else
                {
                    if ( afD[1]*afD[2] > 0.0 )
                        eType = QST_ELLIPTIC_PARABOLOID;
                    else
                        eType = QST_HYPERBOLIC_PARABOLOID;
                }
            }
        }
        else
        {
            // matrix A has one nonzero eigenvalues (d2)
            if ( MgcMath::Abs(kE[0]) > fEpsilon
            ||   MgcMath::Abs(kE[1]) > fEpsilon )
            {
                eType = QST_PARABOLIC_CYLINDER;
            }
            else
            {
                fF = 0.25*kE[2]*kE[2]/afD[2] - m_fC;
                if ( fF*afD[2] > 0.0 )
                    eType = QST_PLANE;
                else
                    eType = QST_NONE;
            }
        }
    }
    else
    {
        // matrix A is (numerically) the zero matrix
        if ( m_kB.Length() > fEpsilon )
            eType = QST_PLANE;
        else
            eType = QST_NONE;
    }

    assert( eType != QST_MAX_TYPE );
}
//----------------------------------------------------------------------------
int MgcQuadricSurface::VertexIndex (const ConvexPolyhedron& rkPoly,
    const Vertex* pkV)
{
    return int(pkV - rkPoly.m_apkVertex);
}
//---------------------------------------------------------------------------
int MgcQuadricSurface::EdgeIndex (const ConvexPolyhedron& rkPoly,
    const Edge* pkE)
{
    return int(pkE - rkPoly.m_apkEdge);
}
//---------------------------------------------------------------------------
int MgcQuadricSurface::TriangleIndex (const ConvexPolyhedron& rkPoly,
    const Triangle* pkT)
{
    return int(pkT - rkPoly.m_apkTriangle);
}
//---------------------------------------------------------------------------
int MgcQuadricSurface::AdjacentOrient (const Triangle* pkT,
    const Triangle* pkA)
{
    if ( pkA->m_apkAdjacent[0] == pkT )
        return 0;
    if ( pkA->m_apkAdjacent[1] == pkT )
        return 1;
    return 2;
}
//---------------------------------------------------------------------------
void MgcQuadricSurface::ComputeCentroid (ConvexPolyhedron& rkPoly,
    int iNumVertices)
{
    rkPoly.m_kCentroid = MgcVector3::ZERO;
    for (int i = 0; i < iNumVertices; i++)
        rkPoly.m_kCentroid += *rkPoly.m_apkVertex[i].m_pkPoint;

    float fInvNumVertices = 1.0/iNumVertices;
    rkPoly.m_kCentroid *= fInvNumVertices;
}
//---------------------------------------------------------------------------
void MgcQuadricSurface::RayIntersectSphere (const MgcVector3& rkCen,
    const MgcVector3& rkMid, MgcVector3* pkIntersect)
{
    // ray is X(t) = cen+t*(mid-cen) = cen+t*dir
    MgcVector3 kDir = rkMid - rkCen;

    // Find t > 0 for which Dot(X(t),X(t)) = 1.  This yields a quadratic
    // equation a2*t^2+a1*t+a0 = 0 for which we want the positive root.
    MgcReal fA2 = kDir.SquaredLength();
    MgcReal fA1 = 2.0*kDir.Dot(rkCen);
    MgcReal fA0 = rkCen.SquaredLength() - 1.0;
    MgcReal fDiscr = MgcMath::Sqrt(MgcMath::Abs(fA1*fA1-4.0*fA0*fA2));
    MgcReal fT = (-fA1+fDiscr)/(2.0*fA2);

    *pkIntersect = rkCen + fT*kDir;
}
//---------------------------------------------------------------------------
void MgcQuadricSurface::Expand (unsigned int uiSteps,
    ConvexPolyhedron& rkPoly)
{
    // compute centroid of seed polyhedron
    ComputeCentroid(rkPoly,rkPoly.m_iNumVertices);

    // save seed polyhedron (shallow copy)
    ConvexPolyhedron kSeed = rkPoly;

    // allocate space for fully subdivided polyhedron
    for (unsigned int uiS = 1; uiS <= uiSteps; uiS++)
    {
        rkPoly.m_iNumVertices = rkPoly.m_iNumVertices + rkPoly.m_iNumEdges;
        rkPoly.m_iNumEdges = 2*rkPoly.m_iNumEdges + 3*rkPoly.m_iNumTriangles;
        rkPoly.m_iNumTriangles = 4*rkPoly.m_iNumTriangles;
    }

    rkPoly.m_apkVertex = new Vertex[rkPoly.m_iNumVertices];
    rkPoly.m_apkEdge = new Edge[rkPoly.m_iNumEdges];
    rkPoly.m_apkTriangle = new Triangle[rkPoly.m_iNumTriangles];

    // duplicate data structure

    // duplicate vertices, compute centroid
    int iV, iE;
    for (iV = 0; iV < kSeed.m_iNumVertices; iV++)
    {
        // transfer points from seed to rkPoly
        Vertex& rkVS = kSeed.m_apkVertex[iV];
        Vertex& rkVP = rkPoly.m_apkVertex[iV];
        rkVP.m_pkPoint = rkVS.m_pkPoint;

        // duplicate edge information
        rkVP.m_iNumEdges = rkVS.m_iNumEdges;
        int iQuantity = rkVS.m_iNumEdges*sizeof(Edge*);
        rkVP.m_apkEdge = new Edge*[iQuantity];

        for (iE = 0; iE < rkVS.m_iNumEdges; iE++)
        {
            rkVP.m_apkEdge[iE] = &rkPoly.m_apkEdge[EdgeIndex(kSeed,
                rkVS.m_apkEdge[iE])];
        }

        delete[] rkVS.m_apkEdge;
    }

    delete[] kSeed.m_apkVertex;

    // duplicate edges
    for (iE = 0; iE < kSeed.m_iNumEdges; iE++)
    {
        Edge& rkES = kSeed.m_apkEdge[iE];
        Edge& rkEP = rkPoly.m_apkEdge[iE];

        for (int i = 0; i < 2; i++)
        {
            rkEP.m_apkVertex[i] = &rkPoly.m_apkVertex[
                VertexIndex(kSeed,rkES.m_apkVertex[i])];
            rkEP.m_apkTriangle[i] = &rkPoly.m_apkTriangle[
                TriangleIndex(kSeed,rkES.m_apkTriangle[iE])];
        }

        // For testing purposes, but not necessary for the algorithm.
        // This allows the display program to show the subdivision
        // structure.
        rkEP.m_uiStep = 0;
    }

    delete[] kSeed.m_apkEdge;

    // duplicate triangles
    for (int iT = 0; iT < kSeed.m_iNumTriangles; iT++)
    {
        Triangle& rkTS = kSeed.m_apkTriangle[iT];
        Triangle& rkTP = rkPoly.m_apkTriangle[iT];

        for (int i = 0; i < 3; i++)
        {
            rkTP.m_apkVertex[i] = &rkPoly.m_apkVertex[
                VertexIndex(kSeed,rkTS.m_apkVertex[i])];
            rkTP.m_apkEdge[i] = &rkPoly.m_apkEdge[
                EdgeIndex(kSeed,rkTS.m_apkEdge[i])];
            rkTP.m_apkAdjacent[i] = &rkPoly.m_apkTriangle[
                TriangleIndex(kSeed,rkTS.m_apkAdjacent[i])];
        }
    }

    delete[] kSeed.m_apkTriangle;
}
//---------------------------------------------------------------------------
void MgcQuadricSurface::TessellateSphere (unsigned int uiSteps,
    ConvexPolyhedron& rkPoly)
{
    if ( uiSteps == 0 )
        return;

    // indices to beginning of new slots
    int iVMax = rkPoly.m_iNumVertices;
    int iEMax = rkPoly.m_iNumEdges;
    int iTMax = rkPoly.m_iNumTriangles;

    // increase memory allocation for polyhedron to support subdivision
    Expand(uiSteps,rkPoly);

    // subdivide polyhedron
    for (unsigned int uiS = 1; uiS <= uiSteps; uiS++)
    {
        // Add new point locations to vertices.  Each vertex added in the
        // subdivision always has 6 edges sharing it.  The edge pointers
        // will be set later to point to the subdivided edges.
        //
        // Split old edges in half.  The vertex pointers are set to the new
        // vertices, but the triangle pointers are the old ones and must be
        // updated later to point to the subdivide triangles.
        int iE;
        for (iE = 0; iE < iEMax; iE++)
        {
            Edge* pkE0 = &rkPoly.m_apkEdge[iE];
            Edge* pkE1 = &rkPoly.m_apkEdge[iEMax+iE];

            // get end points of edge
            MgcVector3& rkP0 = *pkE0->m_apkVertex[0]->m_pkPoint;
            MgcVector3& rkP1 = *pkE0->m_apkVertex[1]->m_pkPoint;

            // New vertex is intersection of sphere and ray from polyhedron
            // centroid through midpoint of edge.
            MgcVector3 kMid = 0.5*(rkP0 + rkP1);
            Vertex* pkM = &rkPoly.m_apkVertex[iVMax+iE];
            pkM->m_pkPoint = new MgcVector3;
            RayIntersectSphere(rkPoly.m_kCentroid,kMid,pkM->m_pkPoint);

            // numEdges will be incremented dynamically to 6 later on
            pkM->m_iNumEdges = 0;
            pkM->m_apkEdge = new Edge*[6];

            // edge[iE] = <V0,V1>, M = (V0+V1)/2, E0 = <V0,M>, E1 = <V1,M>
            pkE1->m_apkVertex[0] = pkE0->m_apkVertex[1];
            pkE0->m_apkVertex[1] = pkM;
            pkE1->m_apkVertex[1] = pkM;