wmlintpqdrnonuniform2.cpp

3D Game Engine Design Source Code非常棒

    m_akTData[iT].m_akIntersect[2].X() = (fM11*fR0-fM01*fR1)*fInvDet;
    m_akTData[iT].m_akIntersect[2].Y() = (fM00*fR1-fM10*fR0)*fInvDet;
}
//----------------------------------------------------------------------------
template <class Real>
void IntpQdrNonuniform2<Real>::ComputeCoefficients (int iT)
{
    typename Delaunay2<Real>::Triangle& rkTri = m_akTriangle[iT];
    TriangleData& rkTData = m_akTData[iT];

    // get sample data at main triangle vertices
    Vector2<Real> akV[3];
    Jet afJet[3];
    int i;
    for (i = 0; i < 3; i++)
    {
        int j = rkTri.m_aiVertex[i];
        akV[i] = m_akVertex[j];
        afJet[i].m_fF = m_afF[j];
        afJet[i].m_fFx = m_afFx[j];
        afJet[i].m_fFy = m_afFy[j];
    }

    Vector2<Real> akU[3];
    for (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];
    }

    // compute intermediate terms
    Real afCenT[3], afCen0[3], afCen1[3], afCen2[3];

    ComputeBarycenter(akV[0],akV[1],akV[2],rkTData.m_kCenter,afCenT);
    ComputeBarycenter(akV[0],akV[1],akV[2],akU[0],afCen0);
    ComputeBarycenter(akV[0],akV[1],akV[2],akU[1],afCen1);
    ComputeBarycenter(akV[0],akV[1],akV[2],akU[2],afCen2);

    Real fAlpha = (afCenT[1]*afCen1[0]-afCenT[0]*afCen1[1]) /
        (afCen1[0]-afCenT[0]);
    Real fBeta = (afCenT[2]*afCen2[1]-afCenT[1]*afCen2[2]) /
        (afCen2[1]-afCenT[1]);
    Real fGamma = (afCenT[0]*afCen0[2]-afCenT[2]*afCen0[0]) /
        (afCen0[2]-afCenT[2]);
    Real fOmAlpha = (Real)1.0 - fAlpha;
    Real fOmBeta  = (Real)1.0 - fBeta;
    Real fOmGamma = (Real)1.0 - fGamma;

    Real fTmp, afA[9], afB[9];

    fTmp = afCenT[0]*akV[0].X()+afCenT[1]*akV[1].X()+afCenT[2]*akV[2].X();
    afA[0] = ((Real)0.5)*(fTmp-akV[0].X());
    afA[1] = ((Real)0.5)*(fTmp-akV[1].X());
    afA[2] = ((Real)0.5)*(fTmp-akV[2].X());
    afA[3] = ((Real)0.5)*fBeta*(akV[2].X()-akV[0].X());
    afA[4] = ((Real)0.5)*fOmGamma*(akV[1].X()-akV[0].X());
    afA[5] = ((Real)0.5)*fGamma*(akV[0].X()-akV[1].X());
    afA[6] = ((Real)0.5)*fOmAlpha*(akV[2].X()-akV[1].X());
    afA[7] = ((Real)0.5)*fAlpha*(akV[1].X()-akV[2].X());
    afA[8] = ((Real)0.5)*fOmBeta*(akV[0].X()-akV[2].X());

    fTmp = afCenT[0]*akV[0].Y()+afCenT[1]*akV[1].Y()+afCenT[2]*akV[2].Y();
    afB[0] = ((Real)0.5)*(fTmp-akV[0].Y());
    afB[1] = ((Real)0.5)*(fTmp-akV[1].Y());
    afB[2] = ((Real)0.5)*(fTmp-akV[2].Y());
    afB[3] = ((Real)0.5)*fBeta*(akV[2].Y()-akV[0].Y());
    afB[4] = ((Real)0.5)*fOmGamma*(akV[1].Y()-akV[0].Y());
    afB[5] = ((Real)0.5)*fGamma*(akV[0].Y()-akV[1].Y());
    afB[6] = ((Real)0.5)*fOmAlpha*(akV[2].Y()-akV[1].Y());
    afB[7] = ((Real)0.5)*fAlpha*(akV[1].Y()-akV[2].Y());
    afB[8] = ((Real)0.5)*fOmBeta*(akV[0].Y()-akV[2].Y());

    // compute Bezier coefficients
    rkTData.m_afCoeff[ 2] = afJet[0].m_fF;
    rkTData.m_afCoeff[ 4] = afJet[1].m_fF;
    rkTData.m_afCoeff[ 6] = afJet[2].m_fF;

    rkTData.m_afCoeff[14] = afJet[0].m_fF + afA[0]*afJet[0].m_fFx +
        afB[0]*afJet[0].m_fFy;

    rkTData.m_afCoeff[ 7] = afJet[0].m_fF + afA[3]*afJet[0].m_fFx +
        afB[3]*afJet[0].m_fFy;

    rkTData.m_afCoeff[ 8] = afJet[0].m_fF + afA[4]*afJet[0].m_fFx +
        afB[4]*afJet[0].m_fFy;

    rkTData.m_afCoeff[16] = afJet[1].m_fF + afA[1]*afJet[1].m_fFx +
        afB[1]*afJet[1].m_fFy;

    rkTData.m_afCoeff[ 9] = afJet[1].m_fF + afA[5]*afJet[1].m_fFx +
        afB[5]*afJet[1].m_fFy;

    rkTData.m_afCoeff[10] = afJet[1].m_fF + afA[6]*afJet[1].m_fFx +
        afB[6]*afJet[1].m_fFy;

    rkTData.m_afCoeff[18] = afJet[2].m_fF + afA[2]*afJet[2].m_fFx +
        afB[2]*afJet[2].m_fFy;

    rkTData.m_afCoeff[11] = afJet[2].m_fF + afA[7]*afJet[2].m_fFx +
        afB[7]*afJet[2].m_fFy;

    rkTData.m_afCoeff[12] = afJet[2].m_fF + afA[8]*afJet[2].m_fFx +
        afB[8]*afJet[2].m_fFy;

    rkTData.m_afCoeff[ 5] = fAlpha*rkTData.m_afCoeff[10] +
        fOmAlpha*rkTData.m_afCoeff[11];

    rkTData.m_afCoeff[17] = fAlpha*rkTData.m_afCoeff[16] +
        fOmAlpha*rkTData.m_afCoeff[18];

    rkTData.m_afCoeff[ 1] = fBeta*rkTData.m_afCoeff[12] +
        fOmBeta*rkTData.m_afCoeff[ 7];

    rkTData.m_afCoeff[13] = fBeta*rkTData.m_afCoeff[18] +
        fOmBeta*rkTData.m_afCoeff[14];

    rkTData.m_afCoeff[ 3] = fGamma*rkTData.m_afCoeff[ 8] +
        fOmGamma*rkTData.m_afCoeff[ 9];

    rkTData.m_afCoeff[15] = fGamma*rkTData.m_afCoeff[14] +
        fOmGamma*rkTData.m_afCoeff[16]; 

    rkTData.m_afCoeff[ 0] = afCenT[0]*rkTData.m_afCoeff[14] + 
        afCenT[1]*rkTData.m_afCoeff[16] + afCenT[2]*rkTData.m_afCoeff[18];
}
//----------------------------------------------------------------------------
template <class Real>
bool IntpQdrNonuniform2<Real>::Evaluate (const Vector2<Real>& rkPoint,
    Real& rfF, Real& rfFx, Real& rfFy)
{
    // determine which triangle contains the target point
    Vector2<Real> kV0, kV1, kV2;
    int i;
    for (i = 0; i < m_iTriangleQuantity; i++)
    {
        typename Delaunay2<Real>::Triangle& rkTri = m_akTriangle[i];
        kV0 = m_akVertex[rkTri.m_aiVertex[0]];
        kV1 = m_akVertex[rkTri.m_aiVertex[1]];
        kV2 = m_akVertex[rkTri.m_aiVertex[2]];

        if ( InTriangle(kV0,kV1,kV2,rkPoint) )
            break;
    }
    if ( i == m_iTriangleQuantity )
    {
        // point is outside interpolation region
        return false;
    }

    // the input point is in this triangle
    typename Delaunay2<Real>::Triangle& rkTri = m_akTriangle[i];
    TriangleData& rkTData = m_akTData[i];

    // determine which of the six subtriangles contains the target point
    Vector2<Real> kSub0 = rkTData.m_kCenter;
    Vector2<Real> kSub1;
    Vector2<Real> kSub2 = rkTData.m_akIntersect[2];
    int iIndex;
    for (iIndex = 1; iIndex <= 6; iIndex++)
    {
        kSub1 = kSub2;
        if ( iIndex % 2 )
            kSub2 = m_akVertex[rkTri.m_aiVertex[iIndex/2]];
        else
            kSub2 = rkTData.m_akIntersect[iIndex/2-1];

        if ( InTriangle(kSub0,kSub1,kSub2,rkPoint) )
            break;
    }

    // This should not happen theoretically, but it can happen due to
    // numerical round-off errors.  Just in case, select an index and go
    // with it.  Probably better is to keep track of the dot products in
    // InTriangle and find the one closest to zero and use a triangle that
    // contains the edge as the one that contains the input point.
    assert( iIndex <= 6 );
    if ( iIndex > 6 )
        iIndex = 1;

    // compute barycentric coordinates with respect to subtriangle
    Real afBary[3];
    ComputeBarycenter(kSub0,kSub1,kSub2,rkPoint,afBary);

    // fetch Bezier control points
    Real afBez[6] =
    {
        rkTData.m_afCoeff[0],
        rkTData.m_afCoeff[12 + iIndex],
        rkTData.m_afCoeff[13 + (iIndex % 6)],
        rkTData.m_afCoeff[iIndex],
        rkTData.m_afCoeff[6 + iIndex],
        rkTData.m_afCoeff[1 + (iIndex % 6)]
    };

    // evaluate Bezier quadratic
    rfF = afBary[0]*(afBez[0]*afBary[0] + afBez[1]*afBary[1] +
        afBez[2]*afBary[2]) + afBary[1]*(afBez[1]*afBary[0] +
        afBez[3]*afBary[1] + afBez[4]*afBary[2]) + afBary[2]*(
        afBez[2]*afBary[0] + afBez[4]*afBary[1] + afBez[5]*afBary[2]);

    // evaluate barycentric derivatives of F
    Real fFu = ((Real)2.0)*(afBez[0]*afBary[0] + afBez[1]*afBary[1] +
        afBez[2]*afBary[2]);
    Real fFv = ((Real)2.0)*(afBez[1]*afBary[0] + afBez[3]*afBary[1] +
        afBez[4]*afBary[2]);
    Real fFw = ((Real)2.0)*(afBez[2]*afBary[0] + afBez[4]*afBary[1] +
        afBez[5]*afBary[2]);
    Real fDuw = fFu - fFw;
    Real fDvw = fFv - fFw;

    // convert back to (x,y) coordinates
    Real fM00 = kSub0.X() - kSub2.X();
    Real fM10 = kSub0.Y() - kSub2.Y();
    Real fM01 = kSub1.X() - kSub2.X();
    Real fM11 = kSub1.Y() - kSub2.Y();
    Real fInvDet = ((Real)1.0)/(fM00*fM11 - fM10*fM01);

    rfFx = fInvDet*(fM11*fDuw - fM10*fDvw);
    rfFy = fInvDet*(fM00*fDvw - fM01*fDuw);

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

//----------------------------------------------------------------------------
// explicit instantiation
//----------------------------------------------------------------------------
namespace Wml
{
template class WML_ITEM IntpQdrNonuniform2<float>;
template class WML_ITEM IntpQdrNonuniform2<double>;
}
//----------------------------------------------------------------------------