wmlintpakimauniform2.cpp

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

// 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 "WmlIntpAkimaUniform2.h"
#include "WmlMath.h"
using namespace Wml;

//----------------------------------------------------------------------------
template <class Real>
IntpAkimaUniform2<Real>::IntpAkimaUniform2 (int iXBound, int iYBound,
    Real fXMin, Real fXSpacing, Real fYMin, Real fYSpacing, Real** aafF)
{
    // At least a 3x3 block of data points are needed to construct the
    // estimates of the boundary derivatives.
    assert( iXBound >= 3 && iYBound >= 3 && aafF );
    assert( fXSpacing > (Real)0.0 && fYSpacing > (Real)0.0 );

    m_iXBound = iXBound;
    m_iYBound = iYBound;
    m_fXMin = fXMin;
    m_fXSpacing = fXSpacing;
    m_fYMin = fYMin;
    m_fYSpacing = fYSpacing;
    m_aafF = aafF;

    int iXBoundM1 = iXBound - 1, iYBoundM1 = iYBound - 1;
    m_iQuantity = iXBound*iYBound;
    m_fXMax = fXMin + fXSpacing*iXBoundM1;
    m_fYMax = fYMin + fYSpacing*iYBoundM1;

    // compute slopes
    Real fInvDX = ((Real)1.0)/fXSpacing, fInvDY = ((Real)1.0)/fYSpacing;
    Real fInvDXDY = fInvDX*fInvDY;
    Real** aafXSlope;
    Allocate2D(iXBound+3,iYBound,aafXSlope);  // xslope[y][x]
    Real** aafYSlope;
    Allocate2D(iYBound+3,iXBound,aafYSlope);  // yslope[x][y]

    int iX, iY;
    for (iY = 0; iY < iYBound; iY++)
    {
        for (iX = 0; iX < iXBoundM1; iX++)
            aafXSlope[iY][iX+2] = (aafF[iY][iX+1]-aafF[iY][iX])*fInvDX;

        aafXSlope[iY][1] = ((Real)2.0)*aafXSlope[iY][2] - aafXSlope[iY][3];
        aafXSlope[iY][0] = ((Real)2.0)*aafXSlope[iY][1] - aafXSlope[iY][2];
        aafXSlope[iY][iXBound+1] = ((Real)2.0)*aafXSlope[iY][iXBound] -
            aafXSlope[iY][iXBound-1];
        aafXSlope[iY][iXBound+2] = ((Real)2.0)*aafXSlope[iY][iXBound+1] -
            aafXSlope[iY][iXBound];
    }

    for (iX = 0; iX < iXBound; iX++)
    {
        for (iY = 0; iY < iYBoundM1; iY++)
            aafYSlope[iX][iY+2] = (aafF[iY+1][iX]-aafF[iY][iX])*fInvDY;

        aafYSlope[iX][1] = ((Real)2.0)*aafYSlope[iX][2] - aafYSlope[iX][3];
        aafYSlope[iX][0] = ((Real)2.0)*aafYSlope[iX][1] - aafYSlope[iX][2];
        aafYSlope[iX][iYBound+1] = ((Real)2.0)*aafYSlope[iX][iYBound] -
            aafYSlope[iX][iYBound-1];
        aafYSlope[iX][iYBound+2] = ((Real)2.0)*aafYSlope[iX][iYBound+1] -
            aafYSlope[iX][iYBound];
    }

    // construct first-order derivatives
    Real** aafFX;
    Allocate2D(iXBound,iYBound,aafFX);
    Real** aafFY;
    Allocate2D(iXBound,iYBound,aafFY);

    for (iY = 0; iY < iYBound; iY++)
    {
        for (iX = 0; iX < iXBound; iX++)
            aafFX[iY][iX] = ComputeDerivative(aafXSlope[iY]+iX);
    }

    for (iX = 0; iX < iXBound; iX++)
    {
        for (iY = 0; iY < iYBound; iY++)
            aafFY[iY][iX] = ComputeDerivative(aafYSlope[iX]+iY);
    }

    // construct second-order derivatives
    Real** aafFXY;
    Allocate2D(iXBound,iYBound,aafFXY);

    unsigned int iX0 = iXBoundM1, iX1 = iX0-1,  iX2 = iX1-1;
    unsigned int iY0 = iYBoundM1, iY1 = iY0-1,  iY2 = iY1-1;

    // corners
    aafFXY[0][0] = ((Real)0.25)*fInvDXDY*(
         ((Real) 9.0)*aafF[0][0]
        -((Real)12.0)*aafF[0][1]
        +((Real) 3.0)*aafF[0][2]
        -((Real)12.0)*aafF[1][0]
        +((Real)16.0)*aafF[1][1]
        -((Real) 4.0)*aafF[1][2]
        +((Real) 3.0)*aafF[2][0]
        -((Real) 4.0)*aafF[2][1]
        +             aafF[2][2]);

    aafFXY[0][iXBoundM1] = ((Real)0.25)*fInvDXDY*(
         ((Real) 9.0)*aafF[0][iX0]
        -((Real)12.0)*aafF[0][iX1]
        +((Real) 3.0)*aafF[0][iX2]
        -((Real)12.0)*aafF[1][iX0]
        +((Real)16.0)*aafF[1][iX1]
        -((Real) 4.0)*aafF[1][iX2]
        +((Real) 3.0)*aafF[2][iX0]
        -((Real) 4.0)*aafF[2][iX1]
        +             aafF[2][iX2]);

    aafFXY[iYBoundM1][0] = ((Real)0.25)*fInvDXDY*(
         ((Real)9.0)*aafF[iY0][0]
        -((Real)12.0)*aafF[iY0][1]
        +((Real) 3.0)*aafF[iY0][2]
        -((Real)12.0)*aafF[iY1][0]
        +((Real)16.0)*aafF[iY1][1]
        -((Real) 4.0)*aafF[iY1][2]
        +((Real) 3.0)*aafF[iY2][0]
        -((Real) 4.0)*aafF[iY2][1]
        +             aafF[iY2][2]);

    aafFXY[iYBoundM1][iXBoundM1] = ((Real)0.25)*fInvDXDY*(
         ((Real)9.0)*aafF[iY0][iX0]
        -((Real)12.0)*aafF[iY0][iX1]
        +((Real) 3.0)*aafF[iY0][iX2]
        -((Real)12.0)*aafF[iY1][iX0]
        +((Real)16.0)*aafF[iY1][iX1]
        -((Real) 4.0)*aafF[iY1][iX2]
        +((Real) 3.0)*aafF[iY2][iX0]
        -((Real) 4.0)*aafF[iY2][iX1]
        +             aafF[iY2][iX2]);

    // x-edges
    for (iX = 1; iX < iXBoundM1; iX++)
    {
        aafFXY[0][iX] = ((Real)0.25)*fInvDXDY*(
            ((Real)3.0)*(aafF[0][iX-1] - aafF[0][iX+1])
           -((Real)4.0)*(aafF[1][iX-1] - aafF[1][iX+1])
           +            (aafF[2][iX-1] - aafF[2][iX+1]));

        aafFXY[iYBoundM1][iX] = ((Real)0.25)*fInvDXDY*(
            3.0f*(aafF[iY0][iX-1] - aafF[iY0][iX+1])
           -4.0f*(aafF[iY1][iX-1] - aafF[iY1][iX+1])
           +     (aafF[iY2][iX-1] - aafF[iY2][iX+1]));
    }

    // y-edges
    for (iY = 1; iY < iYBoundM1; iY++)
    {
        aafFXY[iY][0] = ((Real)0.25)*fInvDXDY*(
            ((Real)3.0)*(aafF[iY-1][0] - aafF[iY+1][0])
           -((Real)4.0)*(aafF[iY-1][1] - aafF[iY+1][1])
           +            (aafF[iY-1][2] - aafF[iY+1][2]));

        aafFXY[iY][iXBoundM1] = ((Real)0.25)*fInvDXDY*(
            ((Real)3.0)*(aafF[iY-1][iX0] - aafF[iY+1][iX0])
           -((Real)4.0)*(aafF[iY-1][iX1] - aafF[iY+1][iX1])
           +            (aafF[iY-1][iX2] - aafF[iY+1][iX2]));
    }

    // interior
    for (iY = 1; iY < iYBoundM1; iY++)
    {
        for (iX = 1; iX < iXBoundM1; iX++)
        {
            aafFXY[iY][iX] = ((Real)0.25)*fInvDXDY*(aafF[iY-1][iX-1] -
                aafF[iY-1][iX+1] - aafF[iY+1][iX-1] + aafF[iY+1][iX+1]);
        }
    }

    // construct polynomials
    Allocate2D(iXBoundM1,iYBoundM1,m_aakPoly);
    for (iY = 0; iY < iYBoundM1; iY++)
    {
        for (iX = 0; iX < iXBoundM1; iX++)
        {
            // Note the 'transposing' of the 2x2 blocks (to match notation
            // used in the polynomial definition).
            Real aafG[2][2] =
            {
                {aafF[iY][iX], aafF[iY+1][iX]},
                {aafF[iY][iX+1], aafF[iY+1][iX+1]}
            };

            Real aafGX[2][2] =
            {
                {aafFX[iY][iX], aafFX[iY+1][iX]},
                {aafFX[iY][iX+1], aafFX[iY+1][iX+1]}
            };

            Real aafGY[2][2] =
            {
                {aafFY[iY][iX], aafFY[iY+1][iX]},
                {aafFY[iY][iX+1], aafFY[iY+1][iX+1]}
            };

            Real aafGXY[2][2] =
            {
                {aafFXY[iY][iX], aafFXY[iY+1][iX]},
                {aafFXY[iY][iX+1], aafFXY[iY+1][iX+1]}
            };

            Construct(m_aakPoly[iY][iX],aafG,aafGX,aafGY,aafGXY);
        }
    }

    Deallocate2D(aafXSlope);
    Deallocate2D(aafYSlope);
    Deallocate2D(aafFX);
    Deallocate2D(aafFY);
    Deallocate2D(aafFXY);
}
//----------------------------------------------------------------------------
template <class Real>
IntpAkimaUniform2<Real>::~IntpAkimaUniform2 ()
{
    Deallocate2D(m_aakPoly);
}
//----------------------------------------------------------------------------
template <class Real>
int IntpAkimaUniform2<Real>::GetXBound () const
{
    return m_iXBound;
}
//----------------------------------------------------------------------------
template <class Real>