Approximation.h

Go to the documentation of this file.

/***************************************************************************
 *   Copyright (c) 2005 Imetric 3D GmbH                                    *
 *                                                                         *
 *   This file is part of the FreeCAD CAx development system.              *
 *                                                                         *
 *   This library is free software; you can redistribute it and/or         *
 *   modify it under the terms of the GNU Library General Public           *
 *   License as published by the Free Software Foundation; either          *
 *   version 2 of the License, or (at your option) any later version.      *
 *                                                                         *
 *   This library  is distributed in the hope that it will be useful,      *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU Library General Public License for more details.                  *
 *                                                                         *
 *   You should have received a copy of the GNU Library General Public     *
 *   License along with this library; see the file COPYING.LIB. If not,    *
 *   write to the Free Software Foundation, Inc., 59 Temple Place,         *
 *   Suite 330, Boston, MA  02111-1307, USA                                *
 *                                                                         *
 ***************************************************************************/


#ifndef MESH_APPROXIMATION_H
#define MESH_APPROXIMATION_H

#include <Mod/Mesh/App/WildMagic4/Wm4Vector3.h>
#include <Mod/Mesh/App/WildMagic4/Wm4QuadricSurface.h>
#include <Mod/Mesh/App/WildMagic4/Wm4Eigen.h>
#include <list>
#include <set>
#include <vector>

#include <Base/Vector3D.h>
#include <Base/Matrix.h>

namespace MeshCore {


class MeshExport Approximation
{
public:
    Approximation() : _bIsFitted(false) { }
    virtual ~Approximation(){ Clear(); }
    void AddPoint(const Base::Vector3f &rcVector);
    void AddPoints(const std::vector<Base::Vector3f> &rvPointVect);
    void AddPoints(const std::set<Base::Vector3f> &rsPointSet);
    void AddPoints(const std::list<Base::Vector3f> &rsPointList);
    const std::list<Base::Vector3f>& GetPoints() const { return _vPoints; }
    Base::Vector3f GetGravity() const;
    unsigned long CountPoints() const;
    void Clear();
    float GetLastResult() const;
    virtual float Fit() = 0;
    bool Done() const;

protected:
    static void Convert( const Wm4::Vector3<float>&, Base::Vector3f&);
    static void Convert( const Base::Vector3f&, Wm4::Vector3<float>&);
    void GetMgcVectorArray( std::vector< Wm4::Vector3<float> >& rcPts ) const;

protected:
    std::list< Base::Vector3f > _vPoints; 
    bool _bIsFitted; 
    float _fLastResult; 
};

// -------------------------------------------------------------------------------

class MeshExport PlaneFit : public Approximation
{
public:
    PlaneFit(){};
    virtual ~PlaneFit(){};
    Base::Vector3f GetBase() const;
    Base::Vector3f GetDirU() const;
    Base::Vector3f GetDirV() const;
    Base::Vector3f GetNormal() const;
    float Fit();
    float GetDistanceToPlane(const Base::Vector3f &rcPoint) const;
    float GetStdDeviation() const;
    float GetSignedStdDeviation() const;
    void ProjectToPlane();

protected:
    Base::Vector3f _vBase; 
    Base::Vector3f _vDirU;
    Base::Vector3f _vDirV;
    Base::Vector3f _vDirW; 
};

// -------------------------------------------------------------------------------

class MeshExport QuadraticFit : public Approximation
{
public:
    QuadraticFit() : Approximation() {};
    virtual ~QuadraticFit(){};
    float GetCoeff(unsigned long ulIndex) const;
    const float& GetCoeffArray() const;
    float Fit();

    void CalcZValues(float x, float y, float &dZ1, float &dZ2) const;
    bool GetCurvatureInfo(float x, float y, float z,
                          float &rfCurv0, float &rfCurv1,
                          Base::Vector3f &rkDir0, Base::Vector3f &rkDir1, float &dDistance);

    bool GetCurvatureInfo(float x, float y, float z,
                          float &rfCurv0, float &rfcurv1);
    void CalcEigenValues(float &dLambda1, float &dLambda2, float &dLambda3,
                         Base::Vector3f &clEV1, Base::Vector3f &clEV2, Base::Vector3f &clEV3) const;

protected:
    float _fCoeff[ 10 ];  
};

// -------------------------------------------------------------------------------

class MeshExport SurfaceFit : public PlaneFit
{
public:
    SurfaceFit();
    virtual ~SurfaceFit(){};

    bool GetCurvatureInfo(float x, float y, float z, float &rfCurv0, float &rfCurv1,
                          Base::Vector3f &rkDir0, Base::Vector3f &rkDir1, float &dDistance);
    bool GetCurvatureInfo(float x, float y, float z, float &rfCurv0, float &rfcurv1);
    float Fit();
    float Value(float x, float y) const;

protected:
    float PolynomFit();
    float _fCoeff[ 10 ];  
};

// -------------------------------------------------------------------------------

class FunctionContainer
{
public:
    typedef float (*Function)(float,float,float);
    FunctionContainer(const float *pKoef)
    {
        Assign( pKoef );
/*
        Function oF;     
        Function aoDF[3];
        Function aoD2F[6];

        oF = &F;
        aoDF[0] = &Fx;
        aoDF[1] = &Fy;
        aoDF[2] = &Fz;
        aoD2F[0] = &Fxx;
        aoD2F[1] = &Fxy;
        aoD2F[2] = &Fxz;
        aoD2F[3] = &Fyy;
        aoD2F[4] = &Fyz;
        aoD2F[5] = &Fzz;

        pImplSurf = new Wm4::QuadricSurface<float>( oF, aoDF, aoD2F );*/

        pImplSurf = new Wm4::QuadricSurface<float>( dKoeff );
    }
    void Assign( const float *pKoef )
    {
        for (long ct=0; ct < 10; ct++)
            dKoeff[ ct ] = pKoef[ ct ];
    }
    virtual ~FunctionContainer(){ delete pImplSurf; }
    float& operator[](int idx){ return dKoeff[ idx ]; }
    bool CurvatureInfo(float x, float y, float z, 
                       float &rfCurv0, float &rfCurv1,
                       Wm4::Vector3<float> &rkDir0,  Wm4::Vector3<float> &rkDir1, float &dDistance)
    {
        return pImplSurf->ComputePrincipalCurvatureInfo( Wm4::Vector3<float>(x, y, z),rfCurv0, rfCurv1, rkDir0, rkDir1 );
    }

    Base::Vector3f GetGradient( float x, float y, float z ) const
    {
        Wm4::Vector3<float> grad = pImplSurf->GetGradient( Wm4::Vector3<float>(x, y, z) );
        return Base::Vector3f( grad.X(), grad.Y(), grad.Z() );
    }

    Base::Matrix4D GetHessian( float x, float y, float z ) const
    {
        Wm4::Matrix3<float> hess = pImplSurf->GetHessian( Wm4::Vector3<float>(x, y, z) );
        Base::Matrix4D cMat; cMat.setToUnity();
        cMat[0][0] = hess[0][0]; cMat[0][1] = hess[0][1]; cMat[0][2] = hess[0][2];
        cMat[1][0] = hess[1][0]; cMat[1][1] = hess[1][1]; cMat[1][2] = hess[1][2];
        cMat[2][0] = hess[2][0]; cMat[2][1] = hess[2][1]; cMat[2][2] = hess[2][2];
        return cMat;
    }

    bool CurvatureInfo(float x, float y, float z,
                       float &rfCurv0, float &rfCurv1)
    {
        float dQuot = Fz(x,y,z);
        float zx = - ( Fx(x,y,z) / dQuot );
        float zy = - ( Fy(x,y,z) / dQuot );
        
        float zxx = - ( 2.0f * ( dKoeff[5] + dKoeff[6] * zx * zx + dKoeff[8] * zx ) ) / dQuot;
        float zyy = - ( 2.0f * ( dKoeff[5] + dKoeff[6] * zy * zy + dKoeff[9] * zy ) ) / dQuot;
        float zxy = - ( dKoeff[6] * zx * zy + dKoeff[7] + dKoeff[8] * zy + dKoeff[9] * zx ) / dQuot;

        float dNen = 1 + zx*zx + zy*zy;
        float dNenSqrt = (float)sqrt( dNen );
        float K = ( zxx * zyy - zxy * zxy ) / ( dNen * dNen );
        float H = 0.5f * ( ( 1.0f+zx*zx - 2*zx*zy*zxy + (1.0f+zy*zy)*zxx ) / ( dNenSqrt * dNenSqrt * dNenSqrt ) ) ;

        float dDiscr = (float)sqrt(fabs(H*H-K));
        rfCurv0 = H - dDiscr;
        rfCurv1 = H + dDiscr;

        return true;
    }

    //+++++++++ Quadric +++++++++++++++++++++++++++++++++++++++
    static float F  ( float x, float y, float z ) 
    {
        return (dKoeff[0] + dKoeff[1]*x + dKoeff[2]*y + dKoeff[3]*z +
                dKoeff[4]*x*x + dKoeff[5]*y*y + dKoeff[6]*z*z +
                dKoeff[7]*x*y + dKoeff[8]*x*z + dKoeff[9]*y*z);
    }
  
    //+++++++++ 1. derivations ++++++++++++++++++++++++++++++++
    static float Fx ( float x, float y, float z )
    {
        return( dKoeff[1] + 2.0f*dKoeff[4]*x + dKoeff[7]*y + dKoeff[8]*z );
    }
    static float Fy ( float x, float y, float z ) 
    {
        return( dKoeff[2] + 2.0f*dKoeff[5]*y + dKoeff[7]*x + dKoeff[9]*z );
    }
    static float Fz ( float x, float y, float z ) 
    {
        return( dKoeff[3] + 2.0f*dKoeff[6]*z + dKoeff[8]*x + dKoeff[9]*y );
    }

    //+++++++++ 2. derivations ++++++++++++++++++++++++++++++++
    static float Fxx( float x, float y, float z ) 
    {
        return( 2.0f*dKoeff[4] );
    }
    static float Fxy( float x, float y, float z ) 
    {
        return( dKoeff[7] );
    }
    static float Fxz( float x, float y, float z ) 
    {
        return( dKoeff[8] );
    }
    static float Fyy( float x, float y, float z ) 
    {
        return( 2.0f*dKoeff[5] );
    }
    static float Fyz( float x, float y, float z ) 
    {
        return( dKoeff[9] );
    }
    static float Fzz( float x, float y, float z ) 
    {
        return( 2.0f*dKoeff[6] );
    }
   
protected:
    static float dKoeff[ 10 ];     
    Wm4::ImplicitSurface<float> *pImplSurf;  
private:
    FunctionContainer(){};
};

class MeshExport PolynomialFit : public Approximation
{
public:
    PolynomialFit();

    virtual ~PolynomialFit();
    float Fit();
    float Value(float x, float y) const;

protected:
    float _fCoeff[9];
};

} // namespace MeshCore

#endif // MESH_APPROXIMATION_H