mgcmincircle.cpp

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

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

// indices of points that support current minimum area circle
class Support
{
public:
    int m_iQuantity;
    int m_aiIndex[3];

    bool Contains (int iIndex)
    {
        for (int i = 0; i < m_iQuantity; i++)
        {
            if ( iIndex == m_aiIndex[i] )
                return true;
        }
        return false;
    }
};

// error checking
static const MgcReal gs_fEpsilon = 1e-12;

// All internal minimal circle calculations store the squared radius in the
// radius member of MgcCircle2.  Only at the end is a sqrt computed.

//---------------------------------------------------------------------------
static bool PointInsideCircle (const MgcVector2& rkP, const MgcCircle2& rkC)
{
    MgcVector2 kDiff = rkP - rkC.Center();
    MgcReal fTest = kDiff.SquaredLength() - rkC.Radius();
    return fTest <= gs_fEpsilon;  // theoretical is: test <= 0
}
//----------------------------------------------------------------------------
static MgcCircle2 ExactCircle1 (const MgcVector2& rkP)
{
    MgcCircle2 kMinimal;
    kMinimal.Center() = rkP;
    kMinimal.Radius() = 0.0;
    return kMinimal;
}
//----------------------------------------------------------------------------
static MgcCircle2 ExactCircle2 (const MgcVector2& rkP0,
    const MgcVector2& rkP1)
{
    MgcCircle2 kMinimal;
    kMinimal.Center() = 0.5*(rkP0 + rkP1);
    MgcVector2 kDiff = rkP1 - rkP0;
    kMinimal.Radius() = 0.25*kDiff.SquaredLength();
    return kMinimal;
}
//----------------------------------------------------------------------------
static MgcCircle2 ExactCircle3 (const MgcVector2& rkP0,
    const MgcVector2& rkP1, const MgcVector2& rkP2)
{
    MgcVector2 kE10 = rkP1 - rkP0;
    MgcVector2 kE20 = rkP2 - rkP0;

    MgcReal aafA[2][2];
    aafA[0][0] = kE10.x;  aafA[0][1] = kE10.y;
    aafA[1][0] = kE20.x;  aafA[1][1] = kE20.y;

    MgcReal afB[2];
    afB[0] = 0.5*kE10.SquaredLength();
    afB[1] = 0.5*kE20.SquaredLength();

    MgcCircle2 kMinimal;
    MgcReal fDet = aafA[0][0]*aafA[1][1]-aafA[0][1]*aafA[1][0];
    assert( MgcMath::Abs(fDet) > gs_fEpsilon );

    MgcReal fInvDet = 1.0/fDet;
    MgcVector2 kQ;
    kQ.x = (aafA[1][1]*afB[0]-aafA[0][1]*afB[1])*fInvDet;
    kQ.y = (aafA[0][0]*afB[1]-aafA[1][0]*afB[0])*fInvDet;
    kMinimal.Center() = rkP0 + kQ;
    kMinimal.Radius() = kQ.SquaredLength();

    return kMinimal;
}
//----------------------------------------------------------------------------
static MgcCircle2 UpdateSupport1 (int i, MgcVector2** apkPerm,
    Support& rkSupp)
{
    const MgcVector2& rkP0 = *apkPerm[rkSupp.m_aiIndex[0]];
    const MgcVector2& rkP1 = *apkPerm[i];

    MgcCircle2 kMinimal = ExactCircle2(rkP0,rkP1);
    rkSupp.m_iQuantity = 2;
    rkSupp.m_aiIndex[1] = i;

    return kMinimal;
}
//----------------------------------------------------------------------------
static MgcCircle2 UpdateSupport2 (int i, MgcVector2** apkPerm,
    Support& rkSupp)
{
    const MgcVector2& rkP0 = *apkPerm[rkSupp.m_aiIndex[0]];
    const MgcVector2& rkP1 = *apkPerm[rkSupp.m_aiIndex[1]];
    const MgcVector2& rkP2 = *apkPerm[i];

    MgcCircle2 akC[3];
    MgcReal fMinRSqr = MgcMath::INFINITY;
    int iIndex = -1;

    akC[0] = ExactCircle2(rkP0,rkP2);
    if ( PointInsideCircle(rkP1,akC[0]) )
    {
        fMinRSqr = akC[0].Radius();
        iIndex = 0;
    }

    akC[1] = ExactCircle2(rkP1,rkP2);
    if ( PointInsideCircle(rkP0,akC[1]) )
    {
        if ( akC[1].Radius() < fMinRSqr )
        {
            fMinRSqr = akC[1].Radius();
            iIndex = 1;
        }
    }

    MgcCircle2 kMinimal;

    if ( iIndex != -1 )
    {
        kMinimal = akC[iIndex];
        rkSupp.m_aiIndex[1-iIndex] = i;
    }
    else
    {
        kMinimal = ExactCircle3(rkP0,rkP1,rkP2);
        assert( kMinimal.Radius() <= fMinRSqr );
        rkSupp.m_iQuantity = 3;
        rkSupp.m_aiIndex[2] = i;
    }

    return kMinimal;
}
//----------------------------------------------------------------------------
static MgcCircle2 UpdateSupport3 (int i, MgcVector2** apkPerm,
    Support& rkSupp)
{
    const MgcVector2& rkP0 = *apkPerm[rkSupp.m_aiIndex[0]];
    const MgcVector2& rkP1 = *apkPerm[rkSupp.m_aiIndex[1]];
    const MgcVector2& rkP2 = *apkPerm[rkSupp.m_aiIndex[2]];
    const MgcVector2& rkP3 = *apkPerm[i];

    MgcCircle2 akC[6];
    MgcReal fMinRSqr = MgcMath::INFINITY;
    int iIndex = -1;

    akC[0] = ExactCircle2(rkP0,rkP3);
    if ( PointInsideCircle(rkP1,akC[0]) && PointInsideCircle(rkP2,akC[0]) )
    {
        fMinRSqr = akC[0].Radius();
        iIndex = 0;
    }

    akC[1] = ExactCircle2(rkP1,rkP3);
    if ( PointInsideCircle(rkP0,akC[1]) && PointInsideCircle(rkP2,akC[1]) )
    {
        if ( akC[1].Radius() < fMinRSqr )
        {
            fMinRSqr = akC[1].Radius();
            iIndex = 1;
        }
    }

    akC[2] = ExactCircle2(rkP2,rkP3);
    if ( PointInsideCircle(rkP0,akC[2]) && PointInsideCircle(rkP1,akC[2]) )
    {
        if ( akC[2].Radius() < fMinRSqr )
        {
            fMinRSqr = akC[2].Radius();
            iIndex = 2;
        }
    }

    akC[3] = ExactCircle3(rkP0,rkP1,rkP3);
    if ( PointInsideCircle(rkP2,akC[3]) )
    {
        if ( akC[3].Radius() < fMinRSqr )
        {
            fMinRSqr = akC[3].Radius();
            iIndex = 3;
        }
    }

    akC[4] = ExactCircle3(rkP0,rkP2,rkP3);
    if ( PointInsideCircle(rkP1,akC[4]) )
    {
        if ( akC[4].Radius() < fMinRSqr )
        {
            fMinRSqr = akC[4].Radius();
            iIndex = 4;
        }
    }

    akC[5] = ExactCircle3(rkP1,rkP2,rkP3);
    if ( PointInsideCircle(rkP0,akC[5]) )
    {
        if ( akC[5].Radius() < fMinRSqr )
        {
            fMinRSqr = akC[5].Radius();
            iIndex = 5;
        }
    }

    assert( iIndex != -1 );

    MgcCircle2 kMinimal = akC[iIndex];

    switch ( iIndex )
    {
    case 0:
        rkSupp.m_iQuantity = 2;
        rkSupp.m_aiIndex[1] = i;
        break;
    case 1:
        rkSupp.m_iQuantity = 2;
        rkSupp.m_aiIndex[0] = i;
        break;
    case 2:
        rkSupp.m_iQuantity = 2;
        rkSupp.m_aiIndex[0] = rkSupp.m_aiIndex[2];
        rkSupp.m_aiIndex[1] = i;
        break;
    case 3:
        rkSupp.m_aiIndex[2] = i;
        break;
    case 4:
        rkSupp.m_aiIndex[1] = i;
        break;
    case 5:
        rkSupp.m_aiIndex[0] = i;
        break;
    }

    return kMinimal;
}
//----------------------------------------------------------------------------
typedef MgcCircle2 (*UpdateFunction)(int,MgcVector2**,Support&);
static UpdateFunction Update[4] =
{
    0,
    UpdateSupport1,
    UpdateSupport2,
    UpdateSupport3
};
//----------------------------------------------------------------------------
MgcCircle2 MgcMinCircle (int iQuantity, const MgcVector2* akPoint)
{
    // initialize random number generator
    static bool s_bFirstTime = true;
    if ( s_bFirstTime )
    {
        srand(367);
        s_bFirstTime = false;
    }

    MgcCircle2 kMinimal;
    Support kSupp;

    if ( iQuantity >= 1 )
    {
        // create identity permutation (0,1,...,iQuantity-1)
        MgcVector2** apkPerm = new MgcVector2*[iQuantity];
        int i;
        for (i = 0; i < iQuantity; i++)
            apkPerm[i] = (MgcVector2*) &akPoint[i];
        
        // generate random permutation
        for (i = iQuantity-1; i > 0; i--)
        {
            int j = rand() % (i+1);
            if ( j != i )
            {
                MgcVector2* pSave = apkPerm[i];
                apkPerm[i] = apkPerm[j];
                apkPerm[j] = pSave;
            }
        }
        
        kMinimal = ExactCircle1(*apkPerm[0]);
        kSupp.m_iQuantity = 1;
        kSupp.m_aiIndex[0] = 0;
        i = 1;
        while ( i < iQuantity )
        {
            if ( !kSupp.Contains(i) )
            {
                if ( !PointInsideCircle(*apkPerm[i],kMinimal) )
                {
                    kMinimal = Update[kSupp.m_iQuantity](i,apkPerm,kSupp);
                    i = 0;
                    continue;
                }
            }
            i++;
        }
        
        delete[] apkPerm;
    }
    else
    {
        assert( false );
    }

    kMinimal.Radius() = MgcMath::Sqrt(kMinimal.Radius());
    return kMinimal;
}
//----------------------------------------------------------------------------