mgccirclescribe.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 "MgcCircleScribe.h"
#include "MgcLinearSystem.h"

static const MgcReal gs_fTolerance = 1e-06;

//----------------------------------------------------------------------------
bool MgcCircumscribe (const MgcVector2& rkV0, const MgcVector2& rkV1,
    const MgcVector2& rkV2, MgcCircle2& rkCircle)
{
    MgcVector2 kE10 = rkV1 - rkV0;
    MgcVector2 kE20 = rkV2 - rkV0;

    MgcReal aafA[2][2] =
    {
        kE10.x, kE10.y,
        kE20.x, kE20.y
    };

    MgcReal afB[2] =
    {
        0.5*kE10.SquaredLength(),
        0.5*kE20.SquaredLength()
    };

    MgcLinearSystem kSystem;
    MgcVector2 kSol;
    if ( kSystem.Solve2(aafA,afB,(MgcReal*)&kSol) )
    {
        rkCircle.Center() = rkV0 + kSol;
        rkCircle.Radius() = kSol.Length();
        return true;
    }
    else
    {
        return false;
    }
}
//----------------------------------------------------------------------------
bool MgcCircumscribe (const MgcVector3& rkV0, const MgcVector3& rkV1,
    const MgcVector3& rkV2, MgcCircle3& rkCircle)
{
    MgcVector3 kE02 = rkV0 - rkV2;
    MgcVector3 kE12 = rkV1 - rkV2;
    MgcReal fE02E02 = kE02.Dot(kE02);
    MgcReal fE02E12 = kE02.Dot(kE12);
    MgcReal fE12E12 = kE12.Dot(kE12);
    MgcReal fDet = fE02E02*fE12E12 - fE02E12*fE02E12;
    if ( MgcMath::Abs(fDet) < gs_fTolerance )
        return false;

    MgcReal fHalfInvDet = 0.5/fDet;
    MgcReal fU0 = fHalfInvDet*fE12E12*(fE02E02 - fE02E12);
    MgcReal fU1 = fHalfInvDet*fE02E02*(fE12E12 - fE02E12);
    MgcReal fU2 = 1.0-fU0-fU1;
    MgcVector3 kTmp = fU0*kE02 + fU1*kE12;

    rkCircle.Center() = rkV2 + kTmp;
    rkCircle.Radius() = kTmp.Length();

    rkCircle.N() = kE02.UnitCross(kE12);
    MgcVector3& rkN = rkCircle.N();
    MgcVector3& rkU = rkCircle.U();
    if ( MgcMath::Abs(rkN.x) >= MgcMath::Abs(rkN.y)
    &&   MgcMath::Abs(rkN.x) >= MgcMath::Abs(rkN.z) )
    {
        rkU.x = -rkN.y;
        rkU.y = rkN.x;
        rkU.z = 0.0;
    }
    else
    {
        rkU.x = 0.0;
        rkU.y = rkN.z;
        rkU.z = -rkN.y;
    }

    rkU.Unitize();
    rkCircle.V() = rkCircle.N().Cross(rkCircle.U());

    return true;
}
//----------------------------------------------------------------------------
bool MgcCircumscribe (const MgcVector3& rkV0, const MgcVector3& rkV1,
    const MgcVector3& rkV2, const MgcVector3& rkV3, MgcSphere& rkSphere)
{
    MgcVector3 kE10 = rkV1 - rkV0;
    MgcVector3 kE20 = rkV2 - rkV0;
    MgcVector3 kE30 = rkV3 - rkV0;

    MgcReal aafA[3][3] =
    {
        kE10.x, kE10.y, kE10.z,
        kE20.x, kE20.y, kE20.z,
        kE30.x, kE30.y, kE30.z
    };

    MgcReal afB[3] =
    {
        0.5*kE10.SquaredLength(),
        0.5*kE20.SquaredLength(),
        0.5*kE30.SquaredLength()
    };

    MgcLinearSystem kSystem;
    MgcVector3 kSol;
    if ( kSystem.Solve3(aafA,afB,(MgcReal*)&kSol) )
    {
        rkSphere.Center() = rkV0 + kSol;
        rkSphere.Radius() = kSol.Length();
        return true;
    }
    else
    {
        return false;
    }
}
//----------------------------------------------------------------------------
bool MgcInscribe (const MgcVector2& rkV0, const MgcVector2& rkV1,
    const MgcVector2& rkV2, MgcCircle2& rkCircle)
{
    // edges
    MgcVector2 kE0 = rkV1 - rkV0;
    MgcVector2 kE1 = rkV2 - rkV1;
    MgcVector2 kE2 = rkV0 - rkV2;

    // normals
    MgcVector2 kN0 = kE0.UnitCross();
    MgcVector2 kN1 = kE1.UnitCross();
    MgcVector2 kN2 = kE2.UnitCross();

    MgcReal fA0 = kN1.Dot(kE0);
    if ( MgcMath::Abs(fA0) < gs_fTolerance )
        return false;

    MgcReal fA1 = kN2.Dot(kE1);
    if ( MgcMath::Abs(fA1) < gs_fTolerance )
        return false;

    MgcReal fA2 = kN0.Dot(kE2);
    if ( MgcMath::Abs(fA2) < gs_fTolerance )
        return false;

    MgcReal fInvA0 = 1.0/fA0;
    MgcReal fInvA1 = 1.0/fA1;
    MgcReal fInvA2 = 1.0/fA2;

    rkCircle.Radius() = 1.0/(fInvA0 + fInvA1 + fInvA2);
    rkCircle.Center() = rkCircle.Radius()*(fInvA0*rkV0 + fInvA1*rkV1 +
        fInvA2*rkV2);

    return true;
}
//----------------------------------------------------------------------------
bool MgcInscribe (const MgcVector3& rkV0, const MgcVector3& rkV1,
    const MgcVector3& rkV2, MgcCircle3& rkCircle)
{
    // edges
    MgcVector3 kE0 = rkV1 - rkV0;
    MgcVector3 kE1 = rkV2 - rkV1;
    MgcVector3 kE2 = rkV0 - rkV2;

    // plane normal
    rkCircle.N() = kE1.Cross(kE0);

    // edge normals within the plane
    MgcVector3 kN0 = rkCircle.N().UnitCross(kE0);
    MgcVector3 kN1 = rkCircle.N().UnitCross(kE1);
    MgcVector3 kN2 = rkCircle.N().UnitCross(kE2);

    MgcReal fA0 = kN1.Dot(kE0);
    if ( MgcMath::Abs(fA0) < gs_fTolerance )
        return false;

    MgcReal fA1 = kN2.Dot(kE1);
    if ( MgcMath::Abs(fA1) < gs_fTolerance )
        return false;

    MgcReal fA2 = kN0.Dot(kE2);
    if ( MgcMath::Abs(fA2) < gs_fTolerance )
        return false;

    MgcReal fInvA0 = 1.0/fA0;
    MgcReal fInvA1 = 1.0/fA1;
    MgcReal fInvA2 = 1.0/fA2;

    rkCircle.Radius() = 1.0/(fInvA0 + fInvA1 + fInvA2);
    rkCircle.Center() = rkCircle.Radius()*(fInvA0*rkV0 + fInvA1*rkV1 +
        fInvA2*rkV2);

    rkCircle.U() = kN0;
    rkCircle.V() = rkCircle.N().Cross(rkCircle.U());

    return true;
}
//----------------------------------------------------------------------------
bool MgcInscribe (const MgcVector3& rkV0, const MgcVector3& rkV1,
    const MgcVector3& rkV2, const MgcVector3& rkV3, MgcSphere& rkSphere)
{
    // edges
    MgcVector3 kE10 = rkV1 - rkV0;
    MgcVector3 kE20 = rkV2 - rkV0;
    MgcVector3 kE30 = rkV3 - rkV0;
    MgcVector3 kE21 = rkV2 - rkV1;
    MgcVector3 kE31 = rkV3 - rkV1;

    // normals
    MgcVector3 kN0 = kE31.Cross(kE21);
    MgcVector3 kN1 = kE20.Cross(kE30);
    MgcVector3 kN2 = kE30.Cross(kE10);
    MgcVector3 kN3 = kE10.Cross(kE20);

    // unitize normals
    if ( MgcMath::Abs(kN0.Unitize()) < gs_fTolerance )
        return false;
    if ( MgcMath::Abs(kN1.Unitize()) < gs_fTolerance )
        return false;
    if ( MgcMath::Abs(kN2.Unitize()) < gs_fTolerance )
        return false;
    if ( MgcMath::Abs(kN3.Unitize()) < gs_fTolerance )
        return false;

    MgcReal aafA[3][3] =
    {
        kN1.x-kN0.x, kN1.y-kN0.y, kN1.z-kN0.z,
        kN2.x-kN0.x, kN2.y-kN0.y, kN2.z-kN0.z,
        kN3.x-kN0.x, kN3.y-kN0.y, kN3.z-kN0.z
    };

    MgcReal afB[3] =
    {
        0.0,
        0.0,
        -kN3.Dot(kE30)
    };

    MgcLinearSystem kSystem;
    MgcVector3 kSol;
    if ( kSystem.Solve3(aafA,afB,(MgcReal*)&kSol) )
    {
        rkSphere.Center() = rkV3 + kSol;
        rkSphere.Radius() = MgcMath::Abs(kN0.Dot(kSol));
        return true;
    }
    else
    {
        return false;
    }
}
//----------------------------------------------------------------------------