mgcquadraticnetwork.cpp

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

    fR1  = fM10*m_akTData[iT].m_kCenter.x + fM11*m_akTData[iT].m_kCenter.y;
    fInvDet = 1.0/(fM00*fM11 - fM01*fM10);
    m_akTData[iT].m_akIntersect[2].x = (fM11*fR0-fM01*fR1)*fInvDet;
    m_akTData[iT].m_akIntersect[2].y = (fM00*fR1-fM10*fR0)*fInvDet;
}
//----------------------------------------------------------------------------
void MgcQuadraticNetwork::ComputeCoefficients (int iT)
{
    Triangle& rkTri = m_akTriangle[iT];
    TriangleData& rkTData = m_akTData[iT];

    // get sample data at main triangle vertices
    MgcVector2 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];
    }

    MgcVector2 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
    MgcReal 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);

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

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

    fTmp = afCenT[0]*akV[0].x+afCenT[1]*akV[1].x+afCenT[2]*akV[2].x;
    afA[0] = 0.5f*(fTmp-akV[0].x);
    afA[1] = 0.5f*(fTmp-akV[1].x);
    afA[2] = 0.5f*(fTmp-akV[2].x);
    afA[3] = 0.5f*fBeta*(akV[2].x-akV[0].x);
    afA[4] = 0.5f*fOmGamma*(akV[1].x-akV[0].x);
    afA[5] = 0.5f*fGamma*(akV[0].x-akV[1].x);
    afA[6] = 0.5f*fOmAlpha*(akV[2].x-akV[1].x);
    afA[7] = 0.5f*fAlpha*(akV[1].x-akV[2].x);
    afA[8] = 0.5f*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] = 0.5f*(fTmp-akV[0].y);
    afB[1] = 0.5f*(fTmp-akV[1].y);
    afB[2] = 0.5f*(fTmp-akV[2].y);
    afB[3] = 0.5f*fBeta*(akV[2].y-akV[0].y);
    afB[4] = 0.5f*fOmGamma*(akV[1].y-akV[0].y);
    afB[5] = 0.5f*fGamma*(akV[0].y-akV[1].y);
    afB[6] = 0.5f*fOmAlpha*(akV[2].y-akV[1].y);
    afB[7] = 0.5f*fAlpha*(akV[1].y-akV[2].y);
    afB[8] = 0.5f*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];
}
//----------------------------------------------------------------------------
bool MgcQuadraticNetwork::Evaluate (const MgcVector2& rkPoint, MgcReal& rfF,
    MgcReal& rfFx, MgcReal& rfFy)
{
    // determine which triangle contains the target point
    MgcVector2 kV0, kV1, kV2;
    int i;
    for (i = 0; i < m_iTriangleQuantity; i++)
    {
        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
    Triangle& rkTri = m_akTriangle[i];
    TriangleData& rkTData = m_akTData[i];

    // determine which of the six subtriangles contains the target point
    MgcVector2 kSub0 = rkTData.m_kCenter;
    MgcVector2 kSub1;
    MgcVector2 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
    MgcReal afBary[3];
    ComputeBarycenter(kSub0,kSub1,kSub2,rkPoint,afBary);

    // fetch Bezier control points
    MgcReal 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
    MgcReal fFu = 2.0*(afBez[0]*afBary[0] + afBez[1]*afBary[1] +
        afBez[2]*afBary[2]);
    MgcReal fFv = 2.0*(afBez[1]*afBary[0] + afBez[3]*afBary[1] +
        afBez[4]*afBary[2]);
    MgcReal fFw = 2.0*(afBez[2]*afBary[0] + afBez[4]*afBary[1] +
        afBez[5]*afBary[2]);
    MgcReal fDuw = fFu - fFw;
    MgcReal fDvw = fFv - fFw;

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

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

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