wmlintpqdrnonuniform2.cpp

3D Game Engine Design Source Code非常棒

// 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 "WmlIntpQdrNonuniform2.h"
using namespace Wml;

//----------------------------------------------------------------------------
template <class Real>
IntpQdrNonuniform2<Real>::IntpQdrNonuniform2 (int iVertexQuantity,
    Vector2<Real>* akVertex, Real* afF, Real* afFx, Real* afFy)
    :
    Delaunay2<Real>(iVertexQuantity,akVertex)
{
    assert( afF && afFx && afFy );
    m_afF = afF;
    m_afFx = afFx;
    m_afFy = afFy;
    ProcessTriangles();
}
//----------------------------------------------------------------------------
template <class Real>
IntpQdrNonuniform2<Real>::IntpQdrNonuniform2 (int iVertexQuantity,
    Vector2<Real>* akVertex, Real* afF)
    :
    Delaunay2<Real>(iVertexQuantity,akVertex)
{
    assert( afF );
    m_afF = afF;
    EstimateDerivatives();
    ProcessTriangles();
}
//----------------------------------------------------------------------------
template <class Real>
IntpQdrNonuniform2<Real>::IntpQdrNonuniform2 (Delaunay2<Real>& rkNet,
    Real* afF, Real* afFx, Real* afFy)
    :
    Delaunay2<Real>(rkNet)
{
    assert( afF && afFx && afFy );
    m_afF = afF;
    m_afFx = afFx;
    m_afFy = afFy;
    ProcessTriangles();
}
//----------------------------------------------------------------------------
template <class Real>
IntpQdrNonuniform2<Real>::IntpQdrNonuniform2 (Delaunay2<Real>& rkNet,
    Real* afF)
    :
    Delaunay2<Real>(rkNet)
{
    assert( afF );
    m_afF = afF;
    EstimateDerivatives();
    ProcessTriangles();
}
//----------------------------------------------------------------------------
template <class Real>
IntpQdrNonuniform2<Real>::~IntpQdrNonuniform2 ()
{
    delete[] m_afF;
    delete[] m_afFx;
    delete[] m_afFy;
    delete[] m_akTData;
    delete[] m_akECenter;
}
//----------------------------------------------------------------------------
template <class Real>
void IntpQdrNonuniform2<Real>::EstimateDerivatives ()
{
    m_afFx = new Real[m_iVertexQuantity];
    m_afFy = new Real[m_iVertexQuantity];
    Real* afFz = new Real[m_iVertexQuantity];
    memset(m_afFx,0,m_iVertexQuantity*sizeof(Real));
    memset(m_afFy,0,m_iVertexQuantity*sizeof(Real));
    memset(afFz,0,m_iVertexQuantity*sizeof(Real));

    // accumulate normals at spatial locations (averaging process)
    int i;
    for (i = 0; i < m_iTriangleQuantity; i++)
    {
        typename Delaunay2<Real>::Triangle& rkTri = m_akTriangle[i];

        // get three vertices of triangle
        int j0 = rkTri.m_aiVertex[0];
        int j1 = rkTri.m_aiVertex[1];
        int j2 = rkTri.m_aiVertex[2];

        // compute normal vector of triangle (with positive z-component)
        Real fDx1 = m_akVertex[j1].X() - m_akVertex[j0].X();
        Real fDy1 = m_akVertex[j1].Y() - m_akVertex[j0].Y();
        Real fDz1 = m_afF[j1] - m_afF[j0];
        Real fDx2 = m_akVertex[j2].X() - m_akVertex[j0].X();
        Real fDy2 = m_akVertex[j2].Y() - m_akVertex[j0].Y();
        Real fDz2 = m_afF[j2] - m_afF[j0];
        Real fNx = fDy1*fDz2 - fDy2*fDz1;
        Real fNy = fDz1*fDx2 - fDz2*fDx1;
        Real fNz = fDx1*fDy2 - fDx2*fDy1;
        if ( fNz < (Real)0.0 )
        {
            fNx = -fNx;
            fNy = -fNy;
            fNz = -fNz;
        }

        m_afFx[j0] += fNx;  m_afFy[j0] += fNy;  afFz[j0] += fNz;
        m_afFx[j1] += fNx;  m_afFy[j1] += fNy;  afFz[j1] += fNz;
        m_afFx[j2] += fNx;  m_afFy[j2] += fNy;  afFz[j2] += fNz;
    }

    // scale the normals to form (x,y,-1)
    for (i = 0; i < m_iVertexQuantity; i++)
    {
        if ( Math<Real>::FAbs(afFz[i]) > Math<Real>::EPSILON )
        {
            Real fInv = -((Real)1.0)/afFz[i];
            m_afFx[i] *= fInv;
            m_afFy[i] *= fInv;
        }
        else
        {
            m_afFx[i] = (Real)0.0;
            m_afFy[i] = (Real)0.0;
        }
    }

    delete[] afFz;
}
//----------------------------------------------------------------------------
template <class Real>
void IntpQdrNonuniform2<Real>::ProcessTriangles ()
{
    // Add degenerate triangles to boundary triangles so that interpolation
    // at the boundary can be treated in the same way as interpolation in
    // the interior.

    // add quadratic data to triangle network
    m_akTData = new TriangleData[m_iTriangleQuantity];
    m_akECenter = new Vector2<Real>[m_iExtraTriangleQuantity];

    // compute centers of inscribed circles for triangles
    int iT;
    for (iT = 0; iT < m_iTriangleQuantity; iT++)
    {
        typename Delaunay2<Real>::Triangle& rkTri = m_akTriangle[iT];
        Delaunay2<Real>::ComputeInscribedCenter(
            m_akVertex[rkTri.m_aiVertex[0]],
            m_akVertex[rkTri.m_aiVertex[1]],
            m_akVertex[rkTri.m_aiVertex[2]],
            m_akTData[iT].m_kCenter);
    }

    // compute centers of edges on convex hull
    int iE = 0;
    for (iT = 0; iT < m_iTriangleQuantity; iT++)
    {
        typename Delaunay2<Real>::Triangle& rkTri = m_akTriangle[iT];
        for (int j = 0; j < 3; j++)
        {
            if ( rkTri.m_aiAdjacent[j] >= m_iTriangleQuantity )
            {
                m_akECenter[iE] = ((Real)0.5)*(
                    m_akVertex[rkTri.m_aiVertex[j]] +
                    m_akVertex[rkTri.m_aiVertex[(j+1) % 3]]);
                iE++;
            }
        }
    }

    // compute cross-edge intersections
    for (iT = 0; iT < m_iTriangleQuantity; iT++)
        ComputeCrossEdgeIntersections(iT);

    // compute Bezier coefficients
    for (iT = 0; iT < m_iTriangleQuantity; iT++)
        ComputeCoefficients(iT);
}
//----------------------------------------------------------------------------
template <class Real>
void IntpQdrNonuniform2<Real>::ComputeCrossEdgeIntersections (int iT)
{
    typename Delaunay2<Real>::Triangle& rkTri = m_akTriangle[iT];
    const Vector2<Real>& rkV0 = m_akVertex[rkTri.m_aiVertex[0]];
    const Vector2<Real>& rkV1 = m_akVertex[rkTri.m_aiVertex[1]];
    const Vector2<Real>& rkV2 = m_akVertex[rkTri.m_aiVertex[2]];

    Vector2<Real> akU[3];
    for (int i = 0; i < 3; i++)
    {
        int iA = rkTri.m_aiAdjacent[i];
        if ( iA < m_iTriangleQuantity )
            akU[i] = m_akTData[iA].m_kCenter;
        else
            akU[i] = m_akECenter[iA - m_iTriangleQuantity];
    }

    Real fM00, fM01, fM10, fM11, fR0, fR1, fInvDet;

    // intersection on edge <V0,V1>
    fM00 = rkV0.Y() - rkV1.Y();
    fM01 = rkV1.X() - rkV0.X();
    fM10 = m_akTData[iT].m_kCenter.Y() - akU[0].Y();
    fM11 = akU[0].X() - m_akTData[iT].m_kCenter.X();
    fR0  = fM00*rkV0.X() + fM01*rkV0.Y();
    fR1  = fM10*m_akTData[iT].m_kCenter.X() +
        fM11*m_akTData[iT].m_kCenter.Y();
    fInvDet = ((Real)1.0)/(fM00*fM11 - fM01*fM10);
    m_akTData[iT].m_akIntersect[0].X() = (fM11*fR0-fM01*fR1)*fInvDet;
    m_akTData[iT].m_akIntersect[0].Y() = (fM00*fR1-fM10*fR0)*fInvDet;

    // intersection on edge <V1,V2>
    fM00 = rkV1.Y() - rkV2.Y();
    fM01 = rkV2.X() - rkV1.X();
    fM10 = m_akTData[iT].m_kCenter.Y() - akU[1].Y();
    fM11 = akU[1].X() - m_akTData[iT].m_kCenter.X();
    fR0  = fM00*rkV1.X() + fM01*rkV1.Y();
    fR1  = fM10*m_akTData[iT].m_kCenter.X() +
        fM11*m_akTData[iT].m_kCenter.Y();
    fInvDet = ((Real)1.0)/(fM00*fM11 - fM01*fM10);
    m_akTData[iT].m_akIntersect[1].X() = (fM11*fR0-fM01*fR1)*fInvDet;
    m_akTData[iT].m_akIntersect[1].Y() = (fM00*fR1-fM10*fR0)*fInvDet;

    // intersection on edge <V0,V2>
    fM00 = rkV0.Y() - rkV2.Y();
    fM01 = rkV2.X() - rkV0.X();
    fM10 = m_akTData[iT].m_kCenter.Y() - akU[2].Y();
    fM11 = akU[2].X() - m_akTData[iT].m_kCenter.X();
    fR0  = fM00*rkV0.X() + fM01*rkV0.Y();
    fR1  = fM10*m_akTData[iT].m_kCenter.X() +
        fM11*m_akTData[iT].m_kCenter.Y();
    fInvDet = ((Real)1.0)/(fM00*fM11 - fM01*fM10);