Commit 306ac488 authored by Pierre Kraemer's avatar Pierre Kraemer

ajout Collector

parent 48e27569
......@@ -43,30 +43,29 @@ namespace Geometry
typedef CPULinearSolverTraits< SparseMatrix<double>, FullVector<double> > CPUSolverTraits ;
template <typename PFP>
void computeCurvatureVertices(
void computeCurvatureVertices_QuadraticFitting(
typename PFP::MAP& map,
const typename PFP::TVEC3& position,
const typename PFP::TVEC3& normal,
typename PFP::TREAL& k1,
typename PFP::TREAL& k2,
typename PFP::TVEC3& K1,
typename PFP::TVEC3& K2,
typename PFP::TREAL& kmax,
typename PFP::TREAL& kmin,
typename PFP::TVEC3& Kmax,
typename PFP::TVEC3& Kmin,
const FunctorSelect& select = SelectorTrue()) ;
template <typename PFP>
void computeCurvatureVertex(
void computeCurvatureVertex_QuadraticFitting(
typename PFP::MAP& map,
Dart dart,
const typename PFP::TVEC3& position,
const typename PFP::TVEC3& normal,
typename PFP::TREAL& k1,
typename PFP::TREAL& k2,
typename PFP::TVEC3& K1,
typename PFP::TVEC3& K2) ;
typename PFP::TREAL& kmax,
typename PFP::TREAL& kmin,
typename PFP::TVEC3& Kmax,
typename PFP::TVEC3& Kmin) ;
template <typename PFP>
void vertexQuadricFitting(
void vertexQuadraticFitting(
typename PFP::MAP& map,
Dart dart,
typename PFP::MATRIX33& localFrame,
......@@ -75,10 +74,10 @@ void vertexQuadricFitting(
float& a, float& b, float& c, float& d, float& e) ;
template <typename PFP>
void quadricFittingAddVertexPos(typename PFP::VEC3& v, typename PFP::VEC3& p, typename PFP::MATRIX33& localFrame, LinearSolver<CPUSolverTraits>& solver) ;
void quadraticFittingAddVertexPos(typename PFP::VEC3& v, typename PFP::VEC3& p, typename PFP::MATRIX33& localFrame, LinearSolver<CPUSolverTraits>& solver) ;
template <typename PFP>
void quadricFittingAddVertexNormal(typename PFP::VEC3& v, typename PFP::VEC3& n, typename PFP::VEC3& p, typename PFP::MATRIX33& localFrame, LinearSolver<CPUSolverTraits>& solver) ;
void quadraticFittingAddVertexNormal(typename PFP::VEC3& v, typename PFP::VEC3& n, typename PFP::VEC3& p, typename PFP::MATRIX33& localFrame, LinearSolver<CPUSolverTraits>& solver) ;
/*
template <typename PFP>
void vertexCubicFitting(Dart dart, typename PFP::VEC3& normal, float& a, float& b, float& c, float& d, float& e, float& f, float& g, float& h, float& i) ;
......@@ -90,6 +89,33 @@ template <typename PFP>
void cubicFittingAddVertexNormal(typename PFP::VEC3& v, typename PFP::VEC3& n, typename PFP::VEC3& p, typename PFP::MATRIX33& localFrame) ;
*/
template <typename PFP>
void computeCurvatureVertices_NormalCycles(
typename PFP::MAP& map,
typename PFP::REAL radius,
const typename PFP::TVEC3& position,
const typename PFP::TVEC3& normal,
const typename PFP::TREAL& angles,
typename PFP::TREAL& kmax,
typename PFP::TREAL& kmin,
typename PFP::TVEC3& Kmax,
typename PFP::TVEC3& Kmin,
typename PFP::TVEC3& Knormal,
const FunctorSelect& select = SelectorTrue()) ;
template <typename PFP>
void computeCurvatureVertex_NormalCycles(
typename PFP::MAP& map,
Dart dart,
const typename PFP::TVEC3& position,
const typename PFP::TVEC3& normal,
const typename PFP::TREAL& angles,
typename PFP::TREAL& kmax,
typename PFP::TREAL& kmin,
typename PFP::TVEC3& Kmax,
typename PFP::TVEC3& Kmin,
typename PFP::TVEC3& Knormal) ;
} // namespace Geoemtry
} // namespace Algo
......
......@@ -25,6 +25,7 @@
#include "Algo/Geometry/localFrame.h"
#include "Geometry/matrix.h"
#include "Topology/generic/cellmarker.h"
#include "Algo/Selection/collector.h"
extern "C"
{
......@@ -44,15 +45,15 @@ namespace Geometry
{
template <typename PFP>
void computeCurvatureVertices(
typename PFP::MAP& map,
const typename PFP::TVEC3& position,
const typename PFP::TVEC3& normal,
typename PFP::TREAL& k1,
typename PFP::TREAL& k2,
typename PFP::TVEC3& K1,
typename PFP::TVEC3& K2,
const FunctorSelect& select)
void computeCurvatureVertices_QuadraticFitting(
typename PFP::MAP& map,
const typename PFP::TVEC3& position,
const typename PFP::TVEC3& normal,
typename PFP::TREAL& kmax,
typename PFP::TREAL& kmin,
typename PFP::TVEC3& Kmax,
typename PFP::TVEC3& Kmin,
const FunctorSelect& select)
{
CellMarker marker(map, VERTEX);
for(Dart d = map.begin(); d != map.end(); map.next(d))
......@@ -60,61 +61,62 @@ void computeCurvatureVertices(
if(select(d) && !marker.isMarked(d))
{
marker.mark(d);
computeCurvatureVertex<PFP>(map, d, position, normal, k1, k2, K1, K2) ;
computeCurvatureVertex_QuadraticFitting<PFP>(map, d, position, normal, kmax, kmin, Kmax, Kmin) ;
}
}
}
template <typename PFP>
void computeCurvatureVertex(
typename PFP::MAP& map,
Dart dart,
const typename PFP::TVEC3& position,
const typename PFP::TVEC3& normal,
typename PFP::TREAL& k1,
typename PFP::TREAL& k2,
typename PFP::TVEC3& K1,
typename PFP::TVEC3& K2)
void computeCurvatureVertex_QuadraticFitting(
typename PFP::MAP& map,
Dart dart,
const typename PFP::TVEC3& position,
const typename PFP::TVEC3& normal,
typename PFP::TREAL& kmax,
typename PFP::TREAL& kmin,
typename PFP::TVEC3& Kmax,
typename PFP::TVEC3& Kmin)
{
typedef typename PFP::REAL REAL ;
typedef typename PFP::VEC3 VEC3 ;
typedef typename PFP::MATRIX33 MATRIX33 ;
VEC3 n = normal[dart] ;
MATRIX33 localFrame = Algo::Geometry::vertexLocalFrame<PFP>(map, dart, position, n) ;
MATRIX33 invLocalFrame ;
localFrame.invert(invLocalFrame) ;
REAL a, b, c, d, e;
//vertexCubicFitting(map,dart,localFrame,a,b,c,d,e,f,g,h,i) ;
vertexQuadricFitting<PFP>(map, dart, localFrame, position, normal, a, b, c, d, e) ;
vertexQuadraticFitting<PFP>(map, dart, localFrame, position, normal, a, b, c, d, e) ;
REAL k1_v, k2_v, K1x, K1y ;
REAL kmax_v, kmin_v, Kmax_x, Kmax_y ;
//int res = slaev2_(&e, &f, &g, &maxC, &minC, &dirX, &dirY) ;
/*int res = */slaev2_(&a, &b, &c, &k1_v, &k2_v, &K1x, &K1y) ;
/*int res = */slaev2_(&a, &b, &c, &kmax_v, &kmin_v, &Kmax_x, &Kmax_y) ;
VEC3 K1_v(K1x, K1y, 0.0f) ;
K1_v = invLocalFrame * K1_v ;
VEC3 K2_v = n ^ K1_v ;
VEC3 Kmax_v(Kmax_x, Kmax_y, 0.0f) ;
Kmax_v = invLocalFrame * Kmax_v ;
VEC3 Kmin_v = n ^ Kmax_v ;
if (k1_v < k2_v)
if (kmax_v < kmin_v)
{
k1[dart] = -k1_v ;
k2[dart] = -k2_v ;
K1[dart] = K1_v ;
K2[dart] = K2_v ;
kmax[dart] = -kmax_v ;
kmin[dart] = -kmin_v ;
Kmax[dart] = Kmax_v ;
Kmin[dart] = Kmin_v ;
}
else
{
k1[dart] = -k2_v ;
k2[dart] = -k1_v ;
K1[dart] = K2_v ;
K2[dart] = K1_v ;
kmax[dart] = -kmin_v ;
kmin[dart] = -kmax_v ;
Kmax[dart] = Kmin_v ;
Kmin[dart] = Kmax_v ;
}
}
template <typename PFP>
void vertexQuadricFitting(
void vertexQuadraticFitting(
typename PFP::MAP& map,
Dart dart,
typename PFP::MATRIX33& localFrame,
......@@ -122,11 +124,7 @@ void vertexQuadricFitting(
const typename PFP::TVEC3& normal,
float& a, float& b, float& c, float& d, float& e)
{
typedef typename PFP::REAL REAL ;
typedef typename PFP::VEC3 VEC3 ;
typedef typename PFP::MATRIX33 MATRIX33 ;
VEC3 p = position[dart] ;
typename PFP::VEC3 p = position[dart] ;
LinearSolver<CPUSolverTraits> solver(5) ;
solver.set_least_squares(true) ;
......@@ -135,10 +133,10 @@ void vertexQuadricFitting(
do
{
// 1-ring vertices
VEC3 v = position[map.phi2(it)] ;
quadricFittingAddVertexPos<PFP>(v, p, localFrame, solver) ;
VEC3 n = normal[map.phi2(it)] ;
quadricFittingAddVertexNormal<PFP>(v, n, p, localFrame, solver) ;
typename PFP::VEC3 v = position[map.phi2(it)] ;
quadraticFittingAddVertexPos<PFP>(v, p, localFrame, solver) ;
typename PFP::VEC3 n = normal[map.phi2(it)] ;
quadraticFittingAddVertexNormal<PFP>(v, n, p, localFrame, solver) ;
// 2-ring vertices
// Dart d2 = map.phi1(map.phi1(map.phi2(map.phi1(it)))) ;
// VEC3 v2 = position[d2] ;
......@@ -158,7 +156,7 @@ void vertexQuadricFitting(
}
template <typename PFP>
void quadricFittingAddVertexPos(typename PFP::VEC3& v, typename PFP::VEC3& p, typename PFP::MATRIX33& localFrame, LinearSolver<CPUSolverTraits>& solver)
void quadraticFittingAddVertexPos(typename PFP::VEC3& v, typename PFP::VEC3& p, typename PFP::MATRIX33& localFrame, LinearSolver<CPUSolverTraits>& solver)
{
typename PFP::VEC3 vec = v - p ;
vec = localFrame * vec ;
......@@ -175,33 +173,33 @@ void quadricFittingAddVertexPos(typename PFP::VEC3& v, typename PFP::VEC3& p, ty
}
template <typename PFP>
void quadricFittingAddVertexNormal(typename PFP::VEC3& v, typename PFP::VEC3& n, typename PFP::VEC3& p, typename PFP::MATRIX33& localFrame, LinearSolver<CPUSolverTraits>& solver)
void quadraticFittingAddVertexNormal(typename PFP::VEC3& v, typename PFP::VEC3& n, typename PFP::VEC3& p, typename PFP::MATRIX33& localFrame, LinearSolver<CPUSolverTraits>& solver)
{
typename PFP::VEC3 vec = v - p ;
vec = localFrame * vec ;
typename PFP::VEC3 norm = localFrame * n ;
solver.begin_row() ;
solver.add_coefficient(0, 2.0f*vec[0]) ;
solver.add_coefficient(0, 2.0f * vec[0]) ;
solver.add_coefficient(1, vec[1]) ;
solver.add_coefficient(2, 0) ;
solver.add_coefficient(3, 1.0f) ;
solver.add_coefficient(4, 0) ;
solver.set_right_hand_side(-1.0f*norm[0]/norm[2]) ;
solver.set_right_hand_side(-1.0f * norm[0] / norm[2]) ;
solver.end_row() ;
solver.begin_row() ;
solver.add_coefficient(0, 0) ;
solver.add_coefficient(1, vec[0]) ;
solver.add_coefficient(2, 2.0f*vec[1]) ;
solver.add_coefficient(2, 2.0f * vec[1]) ;
solver.add_coefficient(3, 0) ;
solver.add_coefficient(4, 1.0f) ;
solver.set_right_hand_side(-1.0f*norm[1]/norm[2]) ;
solver.set_right_hand_side(-1.0f * norm[1] / norm[2]) ;
solver.end_row() ;
}
/*
template <typename PFP>
void DifferentialProperties<PFP>::vertexCubicFitting(Dart dart, gmtl::Vec3f& normal, float& a, float& b, float& c, float& d, float& e, float& f, float& g, float& h, float& i)
void vertexCubicFitting(Dart dart, gmtl::Vec3f& normal, float& a, float& b, float& c, float& d, float& e, float& f, float& g, float& h, float& i)
{
gmtl::Matrix33f localFrame, invLocalFrame ;
Algo::Geometry::vertexLocalFrame<PFP>(m_map,dart,normal,localFrame) ;
......@@ -245,7 +243,7 @@ void DifferentialProperties<PFP>::vertexCubicFitting(Dart dart, gmtl::Vec3f& nor
}
template <typename PFP>
void DifferentialProperties<PFP>::cubicFittingAddVertexPos(gmtl::Vec3f& v, gmtl::Vec3f& p, gmtl::Matrix33f& localFrame)
void cubicFittingAddVertexPos(gmtl::Vec3f& v, gmtl::Vec3f& p, gmtl::Matrix33f& localFrame)
{
gmtl::Vec3f vec = v - p ;
vec = localFrame * vec ;
......@@ -266,7 +264,7 @@ void DifferentialProperties<PFP>::cubicFittingAddVertexPos(gmtl::Vec3f& v, gmtl:
}
template <typename PFP>
void DifferentialProperties<PFP>::cubicFittingAddVertexNormal(gmtl::Vec3f& v, gmtl::Vec3f& n, gmtl::Vec3f& p, gmtl::Matrix33f& localFrame)
void cubicFittingAddVertexNormal(gmtl::Vec3f& v, gmtl::Vec3f& n, gmtl::Vec3f& p, gmtl::Matrix33f& localFrame)
{
gmtl::Vec3f vec = v - p ;
vec = localFrame * vec ;
......@@ -299,6 +297,52 @@ void DifferentialProperties<PFP>::cubicFittingAddVertexNormal(gmtl::Vec3f& v, gm
solverC->end_row() ;
}
*/
template <typename PFP>
void computeCurvatureVertices_NormalCycles(
typename PFP::MAP& map,
typename PFP::REAL radius,
const typename PFP::TVEC3& position,
const typename PFP::TVEC3& normal,
const typename PFP::TREAL& angles,
typename PFP::TREAL& kmax,
typename PFP::TREAL& kmin,
typename PFP::TVEC3& Kmax,
typename PFP::TVEC3& Kmin,
typename PFP::TVEC3& Knormal,
const FunctorSelect& select = SelectorTrue())
{
CellMarker marker(map, VERTEX);
for(Dart d = map.begin(); d != map.end(); map.next(d))
{
if(select(d) && !marker.isMarked(d))
{
marker.mark(d);
computeCurvatureVertex_NormalCycles<PFP>(map, d, radius, position, normal, angles, kmax, kmin, Kmax, Kmin, Knormal) ;
}
}
}
template <typename PFP>
void computeCurvatureVertex_NormalCycles(
typename PFP::MAP& map,
Dart dart,
typename PFP::REAL radius,
const typename PFP::TVEC3& position,
const typename PFP::TVEC3& normal,
const typename PFP::TREAL& angles,
typename PFP::TREAL& kmax,
typename PFP::TREAL& kmin,
typename PFP::TVEC3& Kmax,
typename PFP::TVEC3& Kmin,
typename PFP::TVEC3& Knormal)
{
Collector_WithinSphere<PFP> neigh(map, position);
neigh.init(dart, radius) ;
neigh.collect() ;
neigh.computeArea() ;
}
} // namespace Geometry
} // namespace Algo
......
/*******************************************************************************
* CGoGN: Combinatorial and Geometric modeling with Generic N-dimensional Maps *
* version 0.1 *
* Copyright (C) 2009, IGG Team, LSIIT, University of Strasbourg *
* *
* This library is free software; you can redistribute it and/or modify it *
* under the terms of the GNU Lesser General Public License as published by the *
* Free Software Foundation; either version 2.1 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 Lesser General Public License *
* for more details. *
* *
* You should have received a copy of the GNU Lesser General Public License *
* along with this library; if not, write to the Free Software Foundation, *
* Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. *
* *
* Web site: https://iggservis.u-strasbg.fr/CGoGN/ *
* Contact information: cgogn@unistra.fr *
* *
*******************************************************************************/
#ifndef __COLLECTOR_H__
#define __COLLECTOR_H__
/*****************************************
* Class hierarchy :
* Collector (virtual)
* - Collector_WithinSphere
* - Collector_OneRing
****************************************/
namespace CGoGN
{
namespace Selection
{
/*********************************************************
* Collector
*********************************************************/
template <typename PFP>
class Collector
{
protected:
typedef typename PFP::VEC3 VEC3;
typedef typename PFP::REAL REAL;
typename PFP::MAP& map;
const typename PFP::TVEC3& position;
Dart centerDart;
REAL radius;
std::vector<Dart> insideVertices;
std::vector<Dart> insideEdges;
std::vector<Dart> insideFaces;
std::vector<Dart> border;
public:
Collector(typename PFP::MAP& mymap, const typename PFP::TVEC3& myposition);
virtual void init(Dart d, REAL r = 0) = 0;
virtual void collect() = 0;
void sort()
{
std::sort(insideVertices.begin(), insideVertices.end());
std::sort(insideEdges.begin(), insideEdges.end());
std::sort(insideFaces.begin(), insideFaces.end());
std::sort(border.begin(), border.end());
}
typename PFP::MAP& getMap() { return map; }
const AttributeHandler<typename PFP::VEC3>& getPosition() const { return position; }
Dart getCenter() const { return centerDart; }
REAL getRadius() const { return radius; }
const std::vector<Dart>& getInsideVertices() const { return insideVertices; }
const std::vector<Dart>& getInsideEdges() const { return insideEdges; }
const std::vector<Dart>& getInsideFaces() const { return insideFaces; }
const std::vector<Dart>& getBorder() const { return border; }
unsigned int getNbInsideVertices() const { return insideVertices.size(); }
unsigned int getNbInsideEdges() const { return insideEdges.size(); }
unsigned int getNbInsideFaces() const { return insideFaces.size(); }
unsigned int getNbBorder() const { return border.size(); }
template <typename PPFP>
friend std::ostream& operator<<(std::ostream &out, const Collector<PPFP>& c);
};
/*********************************************************
* Collector One Ring
*********************************************************/
/*
* insideVertices = centerDart
* insideEdges = star (edges incident to centerDart)
* insideFaces = triangles incident to centerDart
* border = vertices of 1-ring
* = edges of 1-ring
*/
template <typename PFP>
class Collector_OneRing : public Collector<PFP>
{
public:
Collector_OneRing(typename PFP::MAP& mymap, const typename PFP::TVEC3& myposition) :
Collector<PFP>(mymap, myposition)
{}
void init(Dart d, typename PFP::REAL r = 0);
void collect();
};
/*********************************************************
* Collector Within Sphere
*********************************************************/
/*
* collect all primitives of the connected component containing "centerDart"
* within the sphere of radius "radius" and center "position[centerDart]"
* (hopefully) it defines a 2-manifold (if inserting border-vertices along the border-edges)
*/
template <typename PFP>
class Collector_WithinSphere : public Collector<PFP>
{
protected:
typename PFP::REAL radius_2;
typename PFP::VEC3 centerPosition;
typename PFP::REAL area;
public:
Collector_WithinSphere(typename PFP::MAP& mymap, const typename PFP::TVEC3& myposition) :
Collector<PFP>(mymap,myposition)
{}
void init(Dart d, typename PFP::REAL r = 0);
void collect();
bool isInside(Dart d)
{
return (this->position[d] - centerPosition).norm2() < radius_2;
}
// alpha = coef d'interpolation dans [0 ,1] tel que v= (1-alpha)*pin + alpha*pout
// est le point d'intersection entre la sphère et le segment [pin, pout]
// avec pin = position[din] à l'intérieur de la sphère
// avec pout = position[dout] à l'extérieur de la sphère
typename PFP::REAL intersect_SphereEdge(const Dart din, const Dart dout);
void computeArea();
typename PFP::REAL getArea() const { return area; }
};
} // namespace Selection
} // namespace CGoGN
#include "Algo/Selection/collector.hpp"
#endif
/*******************************************************************************
* CGoGN: Combinatorial and Geometric modeling with Generic N-dimensional Maps *
* version 0.1 *
* Copyright (C) 2009, IGG Team, LSIIT, University of Strasbourg *
* *
* This library is free software; you can redistribute it and/or modify it *
* under the terms of the GNU Lesser General Public License as published by the *
* Free Software Foundation; either version 2.1 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 Lesser General Public License *
* for more details. *
* *
* You should have received a copy of the GNU Lesser General Public License *
* along with this library; if not, write to the Free Software Foundation, *
* Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. *
* *
* Web site: https://iggservis.u-strasbg.fr/CGoGN/ *
* Contact information: cgogn@unistra.fr *
* *
*******************************************************************************/
namespace CGoGN
{
namespace Selection
{
/********************************************
* GENERIC COLLECTOR
********************************************/
template <typename PFP>
Collector<PFP>::Collector(typename PFP::MAP& mymap, const typename PFP::TVEC3& myposition):
map(mymap),
position(myposition)
{}
template <typename PPFP>
std::ostream& operator<<(std::ostream &out, const Collector<PPFP>& c)
{
out << "Collector around " << c.centerDart << std::endl;
out << "insideVertices (" << c.insideVertices.size() << ") = {";
for (unsigned int i=0; i< c.insideVertices.size(); ++i) out << c.insideVertices[i] << " ";
out << "}" << std::endl ;
out << "insideEdges (" << c.insideEdges.size() << ") = {";
for (unsigned int i=0; i< c.insideEdges.size(); ++i) out << c.insideEdges[i] << " ";
out << "}" << std::endl ;
out << "insideFaces (" << c.insideFaces.size() << ") = {";
for (unsigned int i=0; i< c.insideFaces.size(); ++i) out << c.insideFaces[i] << " ";
out << "}" << std::endl ;
out << "border (" << c.border.size() << ") = {";
for (unsigned int i=0; i< c.border.size(); ++i) out << c.border[i] << " ";
out << "}" << std::endl;
return out;
}
/*********************************************************
* Collector One Ring
*********************************************************/
template <typename PFP>
void Collector_OneRing<PFP>::init(Dart d, typename PFP::REAL r)
{
this->centerDart = d;
this->radius = r;
this->insideVertices.clear();
this->insideEdges.clear();
this->insideFaces.clear();
this->border.clear();
}
template <typename PFP>
void Collector_OneRing<PFP>::collect()
{
const Dart end = this->centerDart;
this->insideVertices.push_back(end);
Dart dd = end;
do
{
this->insideEdges.push_back(dd);
this->insideFaces.push_back(dd);
this->border.push_back(this->map.phi1(dd));
dd = this->map.alpha1(dd);
} while(dd != end);
}
/*********************************************************
* Collector Within Sphere
*********************************************************/
template <typename PFP>
void Collector_WithinSphere<PFP>::init(Dart d, typename PFP::REAL r)
{
this->centerDart = d;
this->radius = r;
this->insideVertices.clear();
this->insideEdges.clear();