testapproximation.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.html/free.pdf

#include <MgcCore.pkg>
#include <MgcApproximation.pkg>
#include <MgcNumerics.pkg>

#define DRAND (2.0*rand()/(MgcReal)RAND_MAX-1.0)

//----------------------------------------------------------------------------
void TestLineFit2 ()
{
    const int n = 256;

    MgcVector2* pt = new MgcVector2[n];
    MgcVector2 dir(1.0,2.0);
    MgcVector2 off(0.1,0.1);

    for (int i = 0; i < n; i++) 
    {
        // generate random point on true line
        MgcReal tval = DRAND;
        pt[i].x = off.x + tval * dir.x;
        pt[i].y = off.y + tval * dir.y;

        // generate random perturbation of point
        const MgcReal amp = 0.1;
        pt[i].x += amp*DRAND;
        pt[i].y += amp*DRAND;
    }

    // unitize the true direction
    MgcReal len = sqrt(dir.x*dir.x+dir.y*dir.y);
    if ( len > 0 ) 
    {
        dir.x /= len;
        dir.y /= len;
    }

    // find true line offset which is closest to origin
    MgcReal offDotDir = off.x*dir.x+off.y*dir.y;
    off.x -= offDotDir * dir.x;
    off.y -= offDotDir * dir.y;

    // print the true parameters
    cout << "TRUE:" << endl;
    cout << "offset = (" 
         << off.x << ',' 
         << off.y << ')' << endl;
    cout << "direction = ("
         << dir.x << ','
         << dir.y << ')' << endl;
    cout << endl;

    MgcOrthogonalLineFit(n,pt,off,dir);
    offDotDir = off.x*dir.x+off.y*dir.y;
    off.x -= offDotDir * dir.x;
    off.y -= offDotDir * dir.y;
    cout << "APPROXIMATE:" << endl;
    cout << "offset = (" 
         << off.x << ',' 
         << off.y << ')' << endl;
    cout << "direction = ("
         << dir.x << ','
         << dir.y << ')' << endl;

    delete[] pt;
}
//----------------------------------------------------------------------------
void TestLineFit3 ()
{
    const int n = 256;

    MgcVector3* pt = new MgcVector3[n];
    MgcVector3 dir(1.0,2.0,3.0);
    MgcVector3 off(0.1,0.1,0.1);

    for (int i = 0; i < n; i++) 
    {
        // generate random point on true line
        MgcReal tval = DRAND;
        pt[i].x = off.x + tval * dir.x;
        pt[i].y = off.y + tval * dir.y;
        pt[i].z = off.z + tval * dir.z;

        // generate random perturbation of point
        const MgcReal amp = 0.1;
        pt[i].x += amp*DRAND;
        pt[i].y += amp*DRAND;
        pt[i].z += amp*DRAND;
    }

    // unitize the true direction
    MgcReal len = sqrt(dir.x*dir.x+dir.y*dir.y+dir.z*dir.z);
    if ( len > 0 ) 
    {
        dir.x /= len;
        dir.y /= len;
        dir.z /= len;
    }

    // find true line offset which is closest to origin
    MgcReal offDotDir = off.x*dir.x+off.y*dir.y+off.z*dir.z;
    off.x -= offDotDir * dir.x;
    off.y -= offDotDir * dir.y;
    off.z -= offDotDir * dir.z;

    // print the true parameters
    cout << "TRUE:" << endl;
    cout << "offset = (" 
         << off.x << ',' 
         << off.y << ',' 
         << off.z << ')' << endl;
    cout << "direction = ("
         << dir.x << ','
         << dir.y << ','
         << dir.z << ')' << endl;
    cout << endl;

    MgcOrthogonalLineFit(n,pt,off,dir);
    offDotDir = off.x*dir.x+off.y*dir.y+off.z*dir.z;
    off.x -= offDotDir * dir.x;
    off.y -= offDotDir * dir.y;
    off.z -= offDotDir * dir.z;
    cout << "APPROXIMATE:" << endl;
    cout << "offset = (" 
         << off.x << ',' 
         << off.y << ',' 
         << off.z << ')' << endl;
    cout << "direction = ("
         << dir.x << ','
         << dir.y << ','
         << dir.z << ')' << endl;

    delete[] pt;
}
//----------------------------------------------------------------------------
void TestPlaneFit ()
{
    const int n = 256;

    MgcVector3* pt = new MgcVector3[n];
    MgcVector3 nor( 1.0, 2.0, 3.0 );
    MgcVector3 off( 0.1, 0.1, 0.1 );

    MgcVector3 tan1( -2.0, 1.0, 0.0 );
    MgcVector3 tan2( -3.0, -6.0, 5.0 );
    for (int i = 0; i < n; i++) 
    {
        // generate random point on true line
        MgcReal sval = DRAND, tval = DRAND;
        pt[i].x = off.x + sval * tan1.x + tval * tan2.x;
        pt[i].y = off.y + sval * tan1.y + tval * tan2.y;
        pt[i].z = off.z + sval * tan1.z + tval * tan2.z;

        // generate random perturbation of point
        const MgcReal amp = 0.1;
        pt[i].x += amp*DRAND;
        pt[i].y += amp*DRAND;
        pt[i].z += amp*DRAND;
    }

    // unitize the true normal
    MgcReal len = sqrt(nor.x*nor.x+nor.y*nor.y+nor.z*nor.z);
    if ( len > 0 ) 
    {
        nor.x /= len;
        nor.y /= len;
        nor.z /= len;
    }

    // find true plane offset which is closest to origin
    MgcReal offDotNor = off.x*nor.x+off.y*nor.y+off.z*nor.z;
    off.x = offDotNor * nor.x;
    off.y = offDotNor * nor.y;
    off.z = offDotNor * nor.z;

    // print the true parameters
    cout << "TRUE:" << endl;
    cout << "offset = (" 
         << off.x << ',' 
         << off.y << ',' 
         << off.z << ')' << endl;
    cout << "normal = ("
         << nor.x << ','
         << nor.y << ','
         << nor.z << ')' << endl;
    cout << endl;

    MgcReal minEnergy = MgcOrthogonalPlaneFit(n,pt,off,nor);
    offDotNor = off.x*nor.x+off.y*nor.y+off.z*nor.z;
    off.x = offDotNor * nor.x;
    off.y = offDotNor * nor.y;
    off.z = offDotNor * nor.z;
    cout << "APPROXIMATE:" << endl;
    cout << "offset = (" 
         << off.x << ',' 
         << off.y << ',' 
         << off.z << ')' << endl;
    cout << "normal = ("
         << nor.x << ','
         << nor.y << ','
         << nor.z << ')' << endl;
    cout << "energy = " << minEnergy << endl;

    delete[] pt;
}
//----------------------------------------------------------------------------
void TestCircleFit0()
{
    // points perturbed slightly from a circle

    // exact circle
    MgcVector2 kCenter(1.0,2.0);
    MgcReal fRadius = 1.0;

    // perturbed points
    const int N = 128;
    MgcVector2 point[N];

    for (int i = 0; i < N; i++)
    {
        MgcReal pa = kCenter.x + 0.1*DRAND;
        MgcReal pb = kCenter.y + 0.1*DRAND;
        MgcReal pr = fRadius + 0.1*DRAND;
        MgcReal angle = DRAND*MgcMath::PI;
        point[i].x = pa + pr*MgcMath::Cos(angle);
        point[i].y = pb + pr*MgcMath::Sin(angle);
    }

    bool bSuccess = MgcCircleFit(N,point,kCenter,fRadius);
}
//----------------------------------------------------------------------------
void TestCircleFit1()
{
    // randomly generated points

    const int N = 128;
    MgcVector2 point[N];

    for (int i = 0; i < N; i++)
    {
        point[i].x = DRAND;
        point[i].y = DRAND;
    }

    MgcVector2 kCenter;
    MgcReal fRadius;
    bool bSuccess = MgcCircleFit(N,point,kCenter,fRadius);
}
//----------------------------------------------------------------------------
void TestSphereFit0 ()
{
    // exact sphere
    MgcVector3 kCenter(1.0,2.0,3.0);
    MgcReal fRadius = 1.0;

    // perturbed points
    const int N = 128;
    MgcVector3 point[N];

    for (int i = 0; i < N; i++)
    {
        MgcReal pa = kCenter.x + 0.1*DRAND;
        MgcReal pb = kCenter.y + 0.1*DRAND;
        MgcReal pc = kCenter.z + 0.1*DRAND;
        MgcReal pr = fRadius + 0.1*DRAND;
        MgcReal theta = DRAND*MgcMath::PI;
        MgcReal phi = 0.5*(DRAND+1.0)*MgcMath::PI;
        MgcReal cst = MgcMath::Cos(theta);
        MgcReal snt = MgcMath::Sin(theta);
        MgcReal csp = MgcMath::Cos(phi);
        MgcReal snp = MgcMath::Sin(phi);
        point[i].x = pa + pr*cst*snp;
        point[i].y = pb + pr*snt*snp;
        point[i].z = pc + pr*csp;
    }

    bool bSuccess = MgcSphereFit(N,point,kCenter,fRadius);
}
//----------------------------------------------------------------------------
void TestSphereFit1 ()
{
    const int N = 128;
    MgcVector3 point[N];

    for (int i = 0; i < N; i++)
    {
        point[i].x = DRAND;
        point[i].y = DRAND;
        point[i].z = DRAND;
    }

    MgcVector3 kCenter;
    float fRadius;
    bool bSuccess = MgcSphereFit(N,point,kCenter,fRadius);
}
//----------------------------------------------------------------------------
void TestParaboloidFit ()
{
    // 9 data samples given on a 3x3 grid
    const int n = 9;
    MgcVector3 point[n];
    point[0].x = 0.0;  point[1].x = 0.5;  point[2].x = 1.0;
    point[3].x = 0.0;  point[4].x = 0.5;  point[5].x = 1.0;
    point[6].x = 0.0;  point[7].x = 0.5;  point[8].x = 1.0;

    point[0].y = 0.0;  point[1].y = 0.0;  point[2].y = 0.0;
    point[3].y = 0.5;  point[4].y = 0.5;  point[5].y = 0.5;
    point[6].y = 1.0;  point[7].y = 1.0;  point[8].y = 1.0;

    point[0].z = 1.0;  point[1].z = 1.0;  point[2].z = 1.0;
    point[3].z = 0.5;  point[4].z = 0.0;  point[5].z = 0.7;
    point[6].z = 0.9;  point[7].z = 1.1;  point[8].z = 0.8;

    MgcReal afCoeff[6];
    bool bSuccess = MgcParaboloidFit(n,point,afCoeff);
}
//----------------------------------------------------------------------------
void TestQuadraticFit2 ()
{
    const int n = 256;
    MgcVector2 p[n];

    // generate points on (x/a)^2 + (y/b)^2 - 1 = 0
    MgcReal a = 1.0;
    MgcReal b = 2.0;

    for (int i = 0; i < n; i++)
    {
        MgcReal angle = MgcMath::TWO_PI*MgcMath::UnitRandom();
        p[i].x = a*MgcMath::Cos(angle);
        p[i].y = b*MgcMath::Sin(angle);
    }

    MgcReal coeff[6];
    MgcReal minEnergy = MgcQuadraticFit(n,p,coeff);
}
//----------------------------------------------------------------------------
void TestQuadraticFit3 ()
{
    const int n = 256;
    MgcVector3 p[n];

    // generate points on (x/a)^2 + (y/b)^2 - (z/c)^2 - 1 = 0
    MgcReal a = 1.0;
    MgcReal b = 2.0;
    MgcReal c = 3.0;

    for (int i = 0; i < n; i++)
    {
        p[i].z = MgcMath::UnitRandom();
        MgcReal ratio = p[i].z/c;
        MgcReal r2 = 1.0+ratio*ratio;
        MgcReal r = MgcMath::Sqrt(r2);
        MgcReal angle = MgcMath::TWO_PI*MgcMath::UnitRandom();
        p[i].x = a*r*MgcMath::Cos(angle);
        p[i].y = b*r*MgcMath::Sin(angle);
    }

    MgcReal coeff[10];
    MgcReal minEnergy = MgcQuadraticFit(n,p,coeff);
}
//----------------------------------------------------------------------------
int main ()
{
    //TestLineFit2();
    //TestLineFit3();
    //TestPlaneFit();
    //TestCircleFit0();
    //TestCircleFit1();
    //TestSphereFit0();
    //TestSphereFit1();
    //TestParaboloidFit();
    TestQuadraticFit2();
    TestQuadraticFit3();
    return 0;
}
//----------------------------------------------------------------------------