mgcimplicitsurface.cpp

3D Game Engine Design Source Code非常棒

// 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 "MgcImplicitSurface.h"
#include "MgcRTLib.h"

//----------------------------------------------------------------------------
MgcImplicitSurface::MgcImplicitSurface (Function oF, Function aoDF[3],
    Function aoD2F[6])
{
    m_oF = oF;
    memcpy(m_aoDF,aoDF,3*sizeof(Function));
    memcpy(m_aoD2F,aoD2F,6*sizeof(Function));
}
//----------------------------------------------------------------------------
bool MgcImplicitSurface::IsOnSurface (MgcReal fX, MgcReal fY,
    MgcReal fZ, MgcReal fTolerance) const
{
    return MgcMath::Abs(m_oF(fX,fY,fZ)) <= fTolerance;
}
//----------------------------------------------------------------------------
MgcVector3 MgcImplicitSurface::GetGradient (MgcReal fX, MgcReal fY,
    MgcReal fZ) const
{
    MgcReal fDx = m_aoDF[0](fX,fY,fZ);
    MgcReal fDy = m_aoDF[1](fX,fY,fZ);
    MgcReal fDz = m_aoDF[2](fX,fY,fZ);
    return MgcVector3(fDx,fDy,fDz);
}
//----------------------------------------------------------------------------
MgcMatrix3 MgcImplicitSurface::GetHessian (MgcReal fX, MgcReal fY,
    MgcReal fZ) const
{
    MgcReal fDxx = m_aoD2F[0](fX,fY,fZ);
    MgcReal fDxy = m_aoD2F[1](fX,fY,fZ);
    MgcReal fDxz = m_aoD2F[2](fX,fY,fZ);
    MgcReal fDyy = m_aoD2F[3](fX,fY,fZ);
    MgcReal fDyz = m_aoD2F[4](fX,fY,fZ);
    MgcReal fDzz = m_aoD2F[5](fX,fY,fZ);
    return MgcMatrix3(fDxx,fDxy,fDxz,fDxy,fDyy,fDyz,fDxz,fDyz,fDzz);
}
//----------------------------------------------------------------------------
void MgcImplicitSurface::GetFrame (MgcReal fX, MgcReal fY, MgcReal fZ,
    MgcVector3& kTangent0, MgcVector3& kTangent1, MgcVector3& kNormal) const
{
    kNormal = GetGradient(fX,fY,fZ);

    if ( MgcMath::Abs(kNormal.x) >= MgcMath::Abs(kNormal.y)
    &&   MgcMath::Abs(kNormal.x) >= MgcMath::Abs(kNormal.z) )
    {
        kTangent0.x = -kNormal.y;
        kTangent0.y = kNormal.x;
        kTangent0.z = 0.0;
    }
    else
    {
        kTangent0.x = 0.0;
        kTangent0.y = kNormal.z;
        kTangent0.z = -kNormal.y;
    }

    kTangent0.Unitize();
    kTangent1 = kNormal.Cross(kTangent0);
}
//----------------------------------------------------------------------------
bool MgcImplicitSurface::ComputePrincipalCurvatureInfo (MgcReal fX,
    MgcReal fY, MgcReal fZ, MgcReal& rfCurv0, MgcReal& rfCurv1,
    MgcVector3& rkDir0, MgcVector3& rkDir1)
{
    // Principal curvatures and directions for implicitly defined surfaces
    // F(x,y,z) = 0.
    //
    // DF = (Fx,Fy,Fz), L = Length(DF)
    //
    // D^2 F = +-           -+
    //         | Fxx Fxy Fxz |
    //         | Fxy Fyy Fyz |
    //         | Fxz Fyz Fzz |
    //         +-           -+
    //
    // adj(D^2 F) = +-                                                 -+
    //              | Fyy*Fzz-Fyz*Fyz  Fyz*Fxz-Fxy*Fzz  Fxy*Fyz-Fxz*Fyy |
    //              | Fyz*Fxz-Fxy*Fzz  Fxx*Fzz-Fxz*Fxz  Fxy*Fxz-Fxx*Fyz |
    //              | Fxy*Fyz-Fxz*Fyy  Fxy*Fxz-Fxx*Fyz  Fxx*Fyy-Fxy*Fxy |
    //              +-                                                 -+
    //
    // Gaussian curvature = [DF^t adj(D^2 F) DF]/L^4
    // 
    // Mean curvature = 0.5*[trace(D^2 F)/L - (DF^t D^2 F DF)/L^3]

    // first derivatives
    MgcVector3 kGradient = GetGradient(fX,fY,fZ);
    MgcReal fL = kGradient.Length();
    const MgcReal fTolerance = 1e-08;
    if ( fL <= fTolerance )
        return false;

    MgcReal fDxDx = kGradient.x*kGradient.x;
    MgcReal fDxDy = kGradient.x*kGradient.y;
    MgcReal fDxDz = kGradient.x*kGradient.z;
    MgcReal fDyDy = kGradient.y*kGradient.y;
    MgcReal fDyDz = kGradient.y*kGradient.z;
    MgcReal fDzDz = kGradient.z*kGradient.z;

    MgcReal fInvL = 1.0/fL;
    MgcReal fInvL2 = fInvL*fInvL;
    MgcReal fInvL3 = fInvL*fInvL2;
    MgcReal fInvL4 = fInvL2*fInvL2;

    // second derivatives
    MgcReal fDxx = m_aoD2F[0](fX,fY,fZ);
    MgcReal fDxy = m_aoD2F[1](fX,fY,fZ);
    MgcReal fDxz = m_aoD2F[2](fX,fY,fZ);
    MgcReal fDyy = m_aoD2F[3](fX,fY,fZ);
    MgcReal fDyz = m_aoD2F[4](fX,fY,fZ);
    MgcReal fDzz = m_aoD2F[5](fX,fY,fZ);

    // mean curvature
    MgcReal fMCurv = 0.5*fInvL3*(fDxx*(fDyDy+fDzDz) + fDyy*(fDxDx+fDzDz) +
        fDzz*(fDxDx+fDyDy) - 2.0*(fDxy*fDxDy+fDxz*fDxDz+fDyz*fDyDz));

    // Gaussian curvature
    MgcReal fGCurv = fInvL4*(fDxDx*(fDyy*fDzz-fDyz*fDyz) +
        fDyDy*(fDxx*fDzz-fDxz*fDxz) + fDzDz*(fDxx*fDyy-fDxy*fDxy) +
        2.0*(fDxDy*(fDxz*fDyz-fDxy*fDzz) + fDxDz*(fDxy*fDyz-fDxz*fDyy) +
        fDyDz*(fDxy*fDxz-fDxx*fDyz)));

    // solve for principal curvatures
    MgcReal fDiscr = MgcMath::Sqrt(fabs(fMCurv*fMCurv-fGCurv));
    rfCurv0 = fMCurv - fDiscr;
    rfCurv1 = fMCurv + fDiscr;

    MgcReal fM00 = ((-1.0 + fDxDx*fInvL2)*fDxx)*fInvL + (fDxDy*fDxy)*fInvL3 +
        (fDxDz*fDxz)*fInvL3;
    MgcReal fM01 = ((-1.0 + fDxDx*fInvL2)*fDxy)*fInvL + (fDxDy*fDyy)*fInvL3 +
        (fDxDz*fDyz)*fInvL3;
    MgcReal fM02 = ((-1.0 + fDxDx*fInvL2)*fDxz)*fInvL + (fDxDy*fDyz)*fInvL3 +
        (fDxDz*fDzz)*fInvL3;
    MgcReal fM10 = (fDxDy*fDxx)*fInvL3 + ((-1.0 + fDyDy*fInvL2)*fDxy)*fInvL +
        (fDyDz*fDxz)*fInvL3;
    MgcReal fM11 = (fDxDy*fDxy)*fInvL3 + ((-1.0 + fDyDy*fInvL2)*fDyy)*fInvL +
        (fDyDz*fDyz)*fInvL3;
    MgcReal fM12 = (fDxDy*fDxz)*fInvL3 + ((-1.0 + fDyDy*fInvL2)*fDyz)*fInvL +
        (fDyDz*fDzz)*fInvL3;
    MgcReal fM20 = (fDxDz*fDxx)*fInvL3 + (fDyDz*fDxy)*fInvL3 + ((-1.0 +
        fDzDz*fInvL2)*fDxz)*fInvL;
    MgcReal fM21 = (fDxDz*fDxy)*fInvL3 + (fDyDz*fDyy)*fInvL3 + ((-1.0 +
        fDzDz*fInvL2)*fDyz)*fInvL;
    MgcReal fM22 = (fDxDz*fDxz)*fInvL3 + (fDyDz*fDyz)*fInvL3 + ((-1.0 +
        fDzDz*fInvL2)*fDzz)*fInvL;

    // solve for principal directions
    MgcReal fTmp1 = fM00 + rfCurv0;
    MgcReal fTmp2 = fM11 + rfCurv0;
    MgcReal fTmp3 = fM22 + rfCurv0;

    MgcVector3 akU[3];
    MgcReal afLength[3];

    akU[0].x = fM01*fM12-fM02*fTmp2;
    akU[0].y = fM02*fM10-fM12*fTmp1;
    akU[0].z = fTmp1*fTmp2-fM01*fM10;
    afLength[0] = akU[0].Length();

    akU[1].x = fM01*fTmp3-fM02*fM21;
    akU[1].y = fM02*fM20-fTmp1*fTmp3;
    akU[1].z = fTmp1*fM21-fM01*fM20;
    afLength[1] = akU[1].Length();

    akU[2].x = fTmp2*fTmp3-fM12*fM21;
    akU[2].y = fM12*fM20-fM10*fTmp3;
    akU[2].z = fM10*fM21-fM20*fTmp2;
    afLength[2] = akU[2].Length();

    int iMaxIndex = 0;
    MgcReal fMax = afLength[0];
    if ( afLength[1] > fMax )
    {
        iMaxIndex = 1;
        fMax = afLength[1];
    }
    if ( afLength[2] > fMax )
        iMaxIndex = 2;

    MgcReal fInvLength = 1.0/afLength[iMaxIndex];
    akU[iMaxIndex] *= fInvLength;

    rkDir1 = akU[iMaxIndex];
    rkDir0 = rkDir1.UnitCross(kGradient);

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