AnimaL

Tutorial

Documentation

---

Documentation

• Overview
• Modules
• Namespace List
• Class Hierarchy
• Alphabetical List
• Compound List
• File List
• Namespace Members
• Compound Members
• File Members
• Examples

---

ParametersCalculus.cpp

Go to the documentation of this file.

#include "ParametersCalculus.h"
#include "SFVec3fCellConstrained.h"
#include "ConstrainedVertex.h"


#include <animal/vectorVec.h>
#include "octree_instanciation.h"

namespace animal
{
namespace octree
{

ParametersCalculus::ParametersCalculus( Cell *cell, HermiteFunction hf ) :
    _cell(cell), _hf(hf)
{
}
ParametersCalculus::~ParametersCalculus()
{
}
    
/*
 * Get alpha, beta, gamma in the cell given the 3d points Q
 */
Vec3d ParametersCalculus::getParameters( Vec3d Q, Vec3d initialParameters, unsigned int maxIterations, FloatingPointType threshold )
{

    Vec3d res;

    
    /*  
    res = initialParameters;

    
    unsigned int it=0;
    CellFrames cf( _cell );
    Vec3d Fx;
    FloatingPointType precNorm = INFINITY;
    do
    {   
        FloatingPointType fAlpha = _hf.compute( res[0] );
        FloatingPointType fBeta = _hf.compute( res[1] );
        FloatingPointType fGamma = _hf.compute( res[2] );

        FloatingPointType fAlpha_ = _hf.computeDerivative( res[0] );
        FloatingPointType fBeta_ = _hf.computeDerivative( res[1] );
        FloatingPointType fGamma_ = _hf.computeDerivative( res[2] );
        

        //dAlphaCoeffs stores the i coefficients associated to
        //f'({1-,}alpha) . f({1-,}beta) . f({1-,}gamma).
        //They're the same for x, y and z
        //As we know that f(1-x) = 1-f(x) so we can say that f'(1-x) = (1-f(x))' = -f'(x)
        FloatingPointType dAlphaCoeffs[8], dBetaCoeffs[8], dGammaCoeffs[8];
        for( unsigned int i=0 ; i<8 ; ++i )
        {
            // (((i%2)==0)?1:-1) will give 1 if we have alpha and -1 if we have 1-alpha
            // And so on
            dAlphaCoeffs[i] = (((i%2)==0)?1:-1)*fAlpha_ * LINEARFACTORBETA(i,fBeta) * LINEARFACTORGAMMA(i,fGamma);
            dBetaCoeffs[i] = LINEARFACTORALPHA(i,fAlpha) * (((i%4)<=1)?1:-1)*fBeta_ * LINEARFACTORGAMMA(i,fGamma);
            dGammaCoeffs[i] = LINEARFACTORALPHA(i,fAlpha) * LINEARFACTORBETA(i,fBeta) * ((i<=3)?1:-1)*fGamma_;
        }

        // Ok, now these are the second part coefficients namely f({1-,}alpha) . f({1-,}beta) . f({1-,}gamma)
        // They're constant for d alpha, d beta and d gamma and the 3 dimensions
        FloatingPointType secondCoeffs[8];
        for( unsigned int i=0 ; i<8 ; ++i )
        {
            secondCoeffs[i] = LINEARFACTORALPHA(i,fAlpha) * LINEARFACTORBETA(i,fBeta) * LINEARFACTORGAMMA(i,fGamma);
        }

        Vec3d Qi[8];
        for( unsigned int i=0 ; i<8 ; ++i )
        {
            Qi[i] = _cell->vertex(i)->getPosition() + LINEARFACTORALPHA(i,res[0])*cf.getFrame(i)[0] + LINEARFACTORBETA(i,res[1])*cf.getFrame(i)[1] + LINEARFACTORGAMMA(i,res[2])*cf.getFrame(i)[2];
        }

        // Now, d Qi(alpha,beta,gamma) / d alpha is quite easily
        // +- ui
        // where ui is the frame's vector associated to x
        // so ((i%2)==0)?1:-1) will give the good sign and cf.getFrame(i)[0] will give ui
        // and so on

        Vec3d jacOnAlpha, jacOnBeta, jacOnGamma;
        for( unsigned int i=0 ; i<8 ; ++i )
        {
            jacOnAlpha += dAlphaCoeffs[i]*Qi[i] + secondCoeffs[i] * (((i%2)==0)?1:-1) * cf.getFrame(i,false)[0];
            jacOnBeta += dBetaCoeffs[i]*Qi[i] + secondCoeffs[i] * (((i%4)<=1)?1:-1) * cf.getFrame(i,false)[1];
            jacOnGamma += dGammaCoeffs[i]*Qi[i] + secondCoeffs[i] * ((i<=3)?1:-1) * cf.getFrame(i,false)[2];
        }
        
        Mat Jx(3,3);
        for(unsigned int i=0 ; i<3 ; ++i )
        {
            Jx[i][0] = jacOnAlpha[i];
            Jx[i][1] = jacOnBeta[i];
            Jx[i][2] = jacOnGamma[i];
        }

        Vec3d pos = computePosition( res );

        Fx = pos - Q;

        
        // Compute delta
        Mat JxtJx(3,3), LU(3,3), LU_1(3,3);
        m_eq_AtB(JxtJx,Jx,Jx);
        m_eq_ludcmp( LU, JxtJx );
    animal::Vec3 aux;
    m_eq_luinv( LU_1, LU, aux );
        
        Mat B(3,3);
        m_eq_ABt(B,LU_1,Jx);
        
        animal::Vec3 Fxv; Fxv[0] = Fx[0]; Fxv[1] = Fx[1]; Fxv[2] = Fx[2];
        m_teq( B, -1.0 );
        
        animal::Vec3 deltaResv;
        v_eq_Ab(deltaResv,B,Fxv);
        
        Vec3d deltaRes; deltaRes[0] = deltaResv[0]; deltaRes[1] = deltaResv[1]; deltaRes[2] = deltaResv[2];
                
        res += deltaRes/2;
        it++;
        
        //Require( Fx.norm() < precNorm );
        precNorm = Fx.norm();
    }
    while( (Fx.norm() > threshold) && (it<maxIterations) );
    
    if( it == maxIterations )
    {
        std::cerr << "MAx iterations reached !\n";
    }
    
    //std::cerr << "Final params are " << res << "\n";
    */
    return res;
}


Vec3d ParametersCalculus::computePosition( Vec3d parameters )
{

Vec3d vRes;
/*
                CellFrames cf( _cell );
                Vec3d P[8];
                for( unsigned short i=0 ; i<8 ; ++i )
                {
                    Vec3d localParams( LINEARFACTORALPHA(i,parameters[0]), LINEARFACTORBETA(i,parameters[1]), LINEARFACTORGAMMA(i,parameters[2]) );
                
                    P[i] = cf.computePosition( i, false, localParams );
                }

                FloatingPointType weightSum = 0.0;
                for( unsigned short i=0 ; i<8 ; ++i )
                {
                    Vec3d localParams( LINEARFACTORALPHA(i,parameters[0]), LINEARFACTORBETA(i,parameters[1]), LINEARFACTORGAMMA(i,parameters[2]) );

                    FloatingPointType weight = hf.compute(localParams[0]) * hf.compute(localParams[1]) * hf.compute(localParams[2]);
                                                    
                    vRes += weight * P[i];
                    weightSum += weight;
                }
                
                vRes /= weightSum;
                
                return vRes;
                
*/
}



}
}

---