PointDistributionTensor.hpp

#ifndef __FIBER_BASEOP_POINTDISTRIBUTIONTENSOR_HPP
#define __FIBER_BASEOP_POINTDISTRIBUTIONTENSOR_HPP
 
#include "PointCloudIpol.hpp"
//#include "elementary/eagle/QuadraticMatrix.hpp"
#include "elementary/eagle/EigenVectors.hpp"
 
namespace Fiber
{
 
namespace PointCloud
{
 
 
template<int N>
inline Eagle::PhysicalSpace::tvector getWeightedVec(const Eagle::PhysicalSpace::tvector& vec, double radius )
{
        return vec;
}
 
template<>
inline Eagle::PhysicalSpace::tvector getWeightedVec<LINEAR>(const Eagle::PhysicalSpace::tvector& vec, double radius )
{
const double len = norm(vec); 
        return vec * (1 - len/radius );
}
template<>
inline Eagle::PhysicalSpace::tvector getWeightedVec<SQUARE>(const Eagle::PhysicalSpace::tvector& vec, double radius )
{
const double len = norm(vec); 
        if(len > 0.0)
                return vec/(len*len); 
 
        return  Eagle::PhysicalSpace::tvector(0,0,0);
}
template<>
inline Eagle::PhysicalSpace::tvector getWeightedVec<SSQUARE>(const Eagle::PhysicalSpace::tvector& vec, double radius )
{
const double len = norm(vec); 
                if(len > 0.0)
                        return vec/(len*len); 
 
                return Eagle::PhysicalSpace::tvector(0,0,0);
}
template<>
inline Eagle::PhysicalSpace::tvector getWeightedVec<SPHCUBIC>(const Eagle::PhysicalSpace::tvector& vec, double radius )
{
const double len = norm(vec); 
                return vec * PointCloud::sph_cubic( len / radius * 2.0 ); 
}
template<>
inline Eagle::PhysicalSpace::tvector getWeightedVec<SPHQUARTIC>(const Eagle::PhysicalSpace::tvector& vec, double radius )
{
const double len = norm(vec); 
        return vec *  PointCloud::sph_quartic( len / radius * 2.5 ); 
}
template<>
inline Eagle::PhysicalSpace::tvector getWeightedVec<SPHQUINTIC>(const Eagle::PhysicalSpace::tvector& vec, double radius )
{
const double len = norm(vec); 
        return vec *  PointCloud::sph_quintic( len / radius * 3.0 ); 
}
 
 
 
 
 
 
template<int N,class C>
Eagle::metric33 computePointDistributionTensor( const std::multimap<double, index_t>& verts,
                                                const C& crds,
                                                const MultiIndex<1>& mi,
                                                const size_t& size,
                                                const double& radius )
{
size_t j = 0; 
Eagle::metric33 tmp;
        tmp(0,0) = 0;
        tmp(1,0) = 0;
        tmp(2,0) = 0;
        tmp(1,1) = 0;
        tmp(2,2) = 0;
        tmp(1,2) = 0;
index_t ind = mi[0];
 
        for (std::multimap<double, index_t>::const_iterator it = verts.begin(); it != verts.end(); it++)
        {
                if( std::fabs((*it).first < 1.0e-4) ) continue;
                if( j == size ) break;
 
        Eagle::PhysicalSpace::tvector vec = crds[ (*it).second ] - crds[ind];
 
                vec = getWeightedVec<N>(vec, radius);
 
                tmp += dyadic(vec, vec);
 
                j++;
        } 
 
        if (j > 0)
                tmp *= 1.0 / double(j);
 
        return tmp;
}
 
template<int N,class C>
Eagle::metric33 computePointDistributionTensor(const C& crds,
                                               const MultiIndex<1>& mi,
                                               const unsigned& size,
                                               const double& radius)
{
Eagle::metric33 tmp;
        tmp(0,0) = 0;
        tmp(1,0) = 0;
        tmp(2,0) = 0;
        tmp(1,1) = 0;
        tmp(2,2) = 0;
        tmp(1,2) = 0;
 
        for (unsigned i = 0; i < size; i++)
        {
        Eagle::PhysicalSpace::tvector vec = crds[ i ] - crds[ mi[0] ];
 
                if(std::fabs(vec[0]) < 1e-4 && std::fabs(vec[1]) < 1e-4 && std::fabs(vec[2]) < 1e-4)
                        continue;
 
                vec = getWeightedVec<N>(vec, radius);
 
                tmp += dyadic(vec, vec);
        }
 
        tmp *= 1.0 / double(size);
 
        return tmp;
}
 
template<int N,class C>
Eagle::metric33 computePointDistributionTensor(const C& crds,
                                               const Eagle::PhysicalSpace::point& p,
                                               const unsigned& size,
                                               const double& radius)
{
Eagle::metric33 tmp;
        tmp(0,0) = 0;
        tmp(1,0) = 0;
        tmp(2,0) = 0;
        tmp(1,1) = 0;
        tmp(2,2) = 0;
        tmp(1,2) = 0;
 
        for (unsigned i = 0; i < size; i++)
        {
        Eagle::PhysicalSpace::tvector vec = crds[ i ] - p;
 
                if(std::fabs(vec[0]) < 1e-4 && std::fabs(vec[1]) < 1e-4 && std::fabs(vec[2]) < 1e-4)
                        continue;
 
                vec = getWeightedVec<N>(vec, radius);
 
                tmp += dyadic(vec, vec);
        }
 
        tmp *= 1.0 / double(size);
 
        return tmp;
}
template <int N, class value>
double getPlanarity(Eagle::Row<N,value>& Evalues )
{
double denominator = Evalues[0]+Evalues[1]+Evalues[2];
 
        if (std::fabs(denominator) < 1.0e-6)
                return 0.0;
        else
                return  2*(Evalues[1] - Evalues[0]) / denominator;
}
 
template <int N, class value>
double getLinearity(Eagle::Row<N,value>& Evalues )
{
double denominator = Evalues[0]+Evalues[1]+Evalues[2];
 
        if (std::fabs(denominator) < 1.0e-6)
                return 0.0;
        else
                return  (Evalues[2] - Evalues[1]) / denominator;
}
 
template <int N, class value>
double getSphericity(Eagle::Row<N,value>& Evalues )
{
double denominator = Evalues[0]+Evalues[1]+Evalues[2];
 
        if (std::fabs(denominator) < 1.0e-6)
                return 0.0;
        else
                return  3*(Evalues[0]) / denominator;
}
 
 
double computePlanarity(const std::vector<Eagle::PhysicalSpace::point>& v,
                        const Eagle::PhysicalSpace::point& p,
                        const double& radius);
 
 
double computeLinearity(const std::vector<Eagle::PhysicalSpace::point>& v,
                        const Eagle::PhysicalSpace::point& p,
                        const double& radius);
 
 
double computeSphericity(const std::vector<Eagle::PhysicalSpace::point>& v,
                         const Eagle::PhysicalSpace::point& p,
                         const double& radius);
 
} //PointCloud
} //Fiber
 
#endif