Commit 418f98f6 authored by Sylvain Thery's avatar Sylvain Thery

Merge branch 'master' of cgogn:~cgogn/CGoGN

parents 8fb58c43 c0fb7beb
......@@ -35,17 +35,29 @@ namespace Geometry
{
template <typename PFP>
typename PFP::REAL triangleArea(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position);
typename PFP::REAL triangleArea(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position) ;
template <typename PFP>
typename PFP::REAL convexFaceArea(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position);
typename PFP::REAL convexFaceArea(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position) ;
template <typename PFP>
typename PFP::REAL totalArea(typename PFP::MAP& map, const typename PFP::TVEC3& position, const FunctorSelect& select = SelectorTrue(), unsigned int th=0) ;
template <typename PFP>
typename PFP::REAL vertexOneRingArea(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position) ;
template <typename PFP>
typename PFP::REAL vertexVoronoiArea(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position) ;
template <typename PFP>
void computeAreaFaces(typename PFP::MAP& map, const typename PFP::TVEC3& position, typename PFP::TREAL& face_area, const FunctorSelect& select = SelectorTrue()) ;
template <typename PFP>
void computeOneRingAreaVertices(typename PFP::MAP& map, const typename PFP::TVEC3& position, typename PFP::TREAL& vertex_area, const FunctorSelect& select = SelectorTrue()) ;
template <typename PFP>
void computeVoronoiAreaVertices(typename PFP::MAP& map, const typename PFP::TVEC3& position, typename PFP::TREAL& vertex_area, const FunctorSelect& select = SelectorTrue()) ;
} // namespace Geometry
} // namespace Algo
......
......@@ -70,7 +70,7 @@ typename PFP::REAL convexFaceArea(typename PFP::MAP& map, Dart d, const typename
template <typename PFP>
typename PFP::REAL totalArea(typename PFP::MAP& map, const typename PFP::TVEC3& position, const FunctorSelect& select, unsigned int th)
{
float area = 0.0f ;
typename PFP::REAL area(0) ;
DartMarker mark(map,th) ;
for(Dart d = map.begin(); d != map.end(); map.next(d))
{
......@@ -83,6 +83,32 @@ typename PFP::REAL totalArea(typename PFP::MAP& map, const typename PFP::TVEC3&
return area ;
}
template <typename PFP>
typename PFP::REAL vertexOneRingArea(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position)
{
typename PFP::REAL area(0) ;
Dart it = d ;
do
{
area += convexFaceArea<PFP>(map, it, position) ;
it = map.alpha1(it) ;
} while(it != d) ;
return area ;
}
template <typename PFP>
typename PFP::REAL vertexVoronoiArea(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position)
{
typename PFP::REAL area(0) ;
Dart it = d ;
do
{
area += convexFaceArea<PFP>(map, it, position) / 3 ;
it = map.alpha1(it) ;
} while(it != d) ;
return area ;
}
template <typename PFP>
void computeAreaFaces(typename PFP::MAP& map, const typename PFP::TVEC3& position, typename PFP::TREAL& face_area, const FunctorSelect& select)
{
......@@ -97,6 +123,34 @@ void computeAreaFaces(typename PFP::MAP& map, const typename PFP::TVEC3& positio
}
}
template <typename PFP>
void computeOneRingAreaVertices(typename PFP::MAP& map, const typename PFP::TVEC3& position, typename PFP::TREAL& vertex_area, const FunctorSelect& select)
{
CellMarker marker(map, VERTEX);
for(Dart d = map.begin(); d != map.end(); map.next(d))
{
if(select(d) && !marker.isMarked(d))
{
marker.mark(d);
vertex_area[d] = vertexOneRingArea<PFP>(map, d, position) ;
}
}
}
template <typename PFP>
void computeVoronoiAreaVertices(typename PFP::MAP& map, const typename PFP::TVEC3& position, typename PFP::TREAL& vertex_area, const FunctorSelect& select)
{
CellMarker marker(map, VERTEX);
for(Dart d = map.begin(); d != map.end(); map.next(d))
{
if(select(d) && !marker.isMarked(d))
{
marker.mark(d);
vertex_area[d] = vertexVoronoiArea<PFP>(map, d, position) ;
}
}
}
} // namespace Geometry
} // namespace Algo
......
......@@ -40,25 +40,48 @@ namespace Algo
namespace Geometry
{
enum LaplacianType
{
TOPOLOGICAL
};
template <typename PFP, typename ATTR_TYPE>
ATTR_TYPE computeLaplacianTopoVertex(
typename PFP::MAP& map,
Dart d,
const AttributeHandler<ATTR_TYPE>& attr) ;
template <typename PFP>
void computeLaplacianVertices(
template <typename PFP, typename ATTR_TYPE>
ATTR_TYPE computeLaplacianCotanVertex(
typename PFP::MAP& map,
LaplacianType type,
const typename PFP::TVEC3& position,
typename PFP::TVEC3& laplacian,
Dart d,
const typename PFP::TREAL& edgeWeight,
const typename PFP::TREAL& vertexArea,
const AttributeHandler<ATTR_TYPE>& attr) ;
template <typename PFP, typename ATTR_TYPE>
void computeLaplacianTopoVertices(
typename PFP::MAP& map,
const AttributeHandler<ATTR_TYPE>& attr,
AttributeHandler<ATTR_TYPE>& laplacian,
const FunctorSelect& select = SelectorTrue()) ;
template <typename PFP, typename ATTR_TYPE>
void computeLaplacianCotanVertices(
typename PFP::MAP& map,
const typename PFP::TREAL& edgeWeight,
const typename PFP::TREAL& vertexArea,
const AttributeHandler<ATTR_TYPE>& attr,
AttributeHandler<ATTR_TYPE>& laplacian,
const FunctorSelect& select = SelectorTrue()) ;
template <typename PFP>
void computeLaplacianVertex_Topo(
typename PFP::REAL computeCotanWeightEdge(
typename PFP::MAP& map,
Dart d,
const typename PFP::TVEC3& position) ;
template <typename PFP>
void computeCotanWeightEdges(
typename PFP::MAP& map,
const typename PFP::TVEC3& position,
typename PFP::TVEC3& laplacian) ;
typename PFP::TREAL& edgeWeight,
const FunctorSelect& select = SelectorTrue()) ;
} // namespace Geoemtry
......
......@@ -33,13 +33,55 @@ namespace Algo
namespace Geometry
{
template <typename PFP>
void computeLaplacianVertices(
typename PFP::MAP& map,
LaplacianType type,
const typename PFP::TVEC3& position,
typename PFP::TVEC3& laplacian,
const FunctorSelect& select)
template <typename PFP, typename ATTR_TYPE>
ATTR_TYPE computeLaplacianTopoVertex(
typename PFP::MAP& map,
Dart d,
const AttributeHandler<ATTR_TYPE>& attr)
{
ATTR_TYPE l(0) ;
unsigned int val = 0 ;
Dart it = d ;
do
{
l += attr[map.phi1(it)] - attr[it] ;
val++ ;
it = map.alpha1(it) ;
} while(it != d) ;
l /= val ;
return l ;
}
template <typename PFP, typename ATTR_TYPE>
ATTR_TYPE computeLaplacianCotanVertex(
typename PFP::MAP& map,
Dart d,
const typename PFP::TREAL& edgeWeight,
const typename PFP::TREAL& vertexArea,
const AttributeHandler<ATTR_TYPE>& attr)
{
typedef typename PFP::REAL REAL ;
ATTR_TYPE l(0) ;
Dart it = d ;
REAL vArea = vertexArea[d] ;
REAL val = 0 ;
do
{
REAL w = edgeWeight[it] / vArea ;
l += (attr[map.phi1(it)] - attr[it]) * w ;
val += w ;
it = map.alpha1(it) ;
} while(it != d) ;
l /= val ;
return l ;
}
template <typename PFP, typename ATTR_TYPE>
void computeLaplacianTopoVertices(
typename PFP::MAP& map,
const AttributeHandler<ATTR_TYPE>& attr,
AttributeHandler<ATTR_TYPE>& laplacian,
const FunctorSelect& select)
{
CellMarker marker(map, VERTEX);
for(Dart d = map.begin(); d != map.end(); map.next(d))
......@@ -47,35 +89,59 @@ void computeLaplacianVertices(
if(select(d) && !marker.isMarked(d))
{
marker.mark(d);
switch(type)
{
case TOPOLOGICAL : {
computeLaplacianVertex_Topo<PFP>(map, d, position, laplacian) ;
break ; }
}
laplacian[d] = computeLaplacianTopoVertex<PFP, ATTR_TYPE>(map, d, attr) ;
}
}
}
template <typename PFP, typename ATTR_TYPE>
void computeLaplacianCotanVertices(
typename PFP::MAP& map,
const typename PFP::TREAL& edgeWeight,
const typename PFP::TREAL& vertexArea,
const AttributeHandler<ATTR_TYPE>& attr,
AttributeHandler<ATTR_TYPE>& laplacian,
const FunctorSelect& select)
{
CellMarker marker(map, VERTEX);
for(Dart d = map.begin(); d != map.end(); map.next(d))
{
if(select(d) && !marker.isMarked(d))
{
marker.mark(d);
laplacian[d] = computeLaplacianCotanVertex<PFP, ATTR_TYPE>(map, d, edgeWeight, vertexArea, attr) ;
}
}
}
template <typename PFP>
void computeLaplacianVertex_Topo(
typename PFP::REAL computeCotanWeightEdge(
typename PFP::MAP& map,
Dart d,
const typename PFP::TVEC3& position)
{
typename PFP::REAL alpha = angle<PFP>(map, map.phi_1(d), map.phi2(map.phi1(d)), position) ;
Dart dd = map.phi2(d) ;
typename PFP::REAL beta = angle<PFP>(map, map.phi_1(dd), map.phi2(map.phi1(dd)), position) ;
return 0.5 * ( 1 / tan(alpha) + 1 / tan(beta) ) ;
}
template <typename PFP>
void computeCotanWeightEdges(
typename PFP::MAP& map,
const typename PFP::TVEC3& position,
typename PFP::TVEC3& laplacian)
typename PFP::TREAL& edgeWeight,
const FunctorSelect& select = SelectorTrue())
{
typedef typename PFP::VEC3 VEC3 ;
VEC3 l(0) ;
unsigned int val = 0 ;
Dart dd = d ;
do
CellMarker marker(map, EDGE);
for(Dart d = map.begin(); d != map.end(); map.next(d))
{
l += vectorOutOfDart<PFP>(map, dd, position) ;
val++ ;
dd = map.alpha1(dd) ;
} while(dd != d) ;
l /= val ;
laplacian[d] = l ;
if(select(d) && !marker.isMarked(d))
{
marker.mark(d);
edgeWeight[d] = computeCotanWeightEdge<PFP>(map, d, position) ;
}
}
}
} // namespace Geometry
......
......@@ -50,24 +50,22 @@ typename PFP::VEC3 triangleNormal(typename PFP::MAP& map, Dart d, const typename
return N ;
}
template<typename PFP>
typename PFP::VEC3 newellNormal(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position)
{
Dart e=d;
Dart e = d;
typename PFP::VEC3 normal(0);
do
{
const typename PFP::VEC3& P = position[e];
e = map.phi1(e);
const typename PFP::VEC3& Q = position[e];
normal[0] += (P[1]-Q[1])*(P[2]+Q[2]);
normal[1] += (P[2]-Q[2])*(P[0]+Q[0]);
normal[2] += (P[0]-Q[0])*(P[1]+Q[1]);
}while (e !=d);
normal[0] += (P[1] - Q[1]) * (P[2] + Q[2]);
normal[1] += (P[2] - Q[2]) * (P[0] + Q[0]);
normal[2] += (P[0] - Q[0]) * (P[1] + Q[1]);
} while (e != d);
normal.normalize();
return normal;
}
......
......@@ -26,6 +26,7 @@
#define __LINEAR_SOLVING_BASIC__
#include "OpenNL/linear_solver.h"
#include "Algo/LinearSolving/variablesSetup.h"
#include "Algo/LinearSolving/matrixSetup.h"
namespace CGoGN
......@@ -34,112 +35,9 @@ namespace CGoGN
namespace LinearSolving
{
template <typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS>
class FunctorMeshToSolver_Scalar : public FunctorMap<typename PFP::MAP>
{
protected:
LinearSolver<SOLVER_TRAITS>* solver ;
CellMarker& lockingMarker ;
const AttributeHandler<ATTR_TYPE>& attrTable ;
const AttributeHandler<unsigned int>& indexTable ;
bool lockedVertices ;
public:
typedef typename PFP::MAP MAP ;
FunctorMeshToSolver_Scalar(MAP& m, LinearSolver<SOLVER_TRAITS>* s, CellMarker& lm, const AttributeHandler<ATTR_TYPE>& attr, const AttributeHandler<unsigned int> index) :
FunctorMap<MAP>(m), solver(s), lockingMarker(lm), attrTable(attr), indexTable(index), lockedVertices(false)
{}
bool operator()(Dart d)
{
solver->variable(indexTable[d]).set_value(attrTable[d]) ;
if(lockingMarker.isMarked(d))
{
solver->variable(indexTable[d]).lock() ;
lockedVertices = true ;
}
return false ;
}
bool hasLockedVertices() { return lockedVertices ; }
} ;
template <typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS>
class FunctorMeshToSolver_Vector : public FunctorMap<typename PFP::MAP>
{
protected:
LinearSolver<SOLVER_TRAITS>* solver ;
CellMarker& lockingMarker ;
const AttributeHandler<ATTR_TYPE>& attrTable ;
const AttributeHandler<unsigned int>& indexTable ;
unsigned int coord ;
bool lockedVertices ;
public:
typedef typename PFP::MAP MAP ;
FunctorMeshToSolver_Vector(MAP& m, LinearSolver<SOLVER_TRAITS>* s, CellMarker& lm, const AttributeHandler<ATTR_TYPE>& attr, unsigned int c, const AttributeHandler<unsigned int> index) :
FunctorMap<MAP>(m), solver(s), lockingMarker(lm), attrTable(attr), indexTable(index), coord(c), lockedVertices(false)
{}
bool operator()(Dart d)
{
solver->variable(indexTable[d]).set_value((attrTable[d])[coord]) ;
if(lockingMarker.isMarked(d))
{
solver->variable(indexTable[d]).lock() ;
lockedVertices = true ;
}
return false ;
}
bool hasLockedVertices() { return lockedVertices ; }
} ;
template <typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS>
class FunctorSolverToMesh_Scalar : public FunctorMap<typename PFP::MAP>
{
protected:
LinearSolver<SOLVER_TRAITS>* solver ;
AttributeHandler<ATTR_TYPE>& attrTable ;
const AttributeHandler<unsigned int>& indexTable ;
public:
typedef typename PFP::MAP MAP ;
FunctorSolverToMesh_Scalar(MAP& m, LinearSolver<SOLVER_TRAITS>* s, AttributeHandler<ATTR_TYPE>& attr, const AttributeHandler<unsigned int> index) :
FunctorMap<MAP>(m), solver(s), attrTable(attr), indexTable(index)
{}
bool operator()(Dart d)
{
attrTable[d] = solver->variable(indexTable[d]).value() ;
return false ;
}
} ;
template <typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS>
class FunctorSolverToMesh_Vector : public FunctorMap<typename PFP::MAP>
{
protected:
LinearSolver<SOLVER_TRAITS>* solver ;
AttributeHandler<ATTR_TYPE>& attrTable ;
const AttributeHandler<unsigned int>& indexTable ;
unsigned int coord ;
public:
typedef typename PFP::MAP MAP ;
FunctorSolverToMesh_Vector(MAP& m, LinearSolver<SOLVER_TRAITS>* s, AttributeHandler<ATTR_TYPE>& attr, unsigned int c, const AttributeHandler<unsigned int> index) :
FunctorMap<MAP>(m), solver(s), attrTable(attr), indexTable(index), coord(c)
{}
bool operator()(Dart d)
{
(attrTable[d])[coord] = solver->variable(indexTable[d]).value() ;
return false ;
}
} ;
/*******************************************************************************
* SOLVER INIT
*******************************************************************************/
template <class SOLVER_TRAITS>
void initSolver(LinearSolver<SOLVER_TRAITS>* s, unsigned int nb_variables, bool least_squares, bool direct)
......@@ -159,96 +57,208 @@ void initSolver(LinearSolver<SOLVER_TRAITS>* s, unsigned int nb_variables, bool
s->set_direct(direct) ;
}
/*******************************************************************************
* VARIABLES SETUP
*******************************************************************************/
template <typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS>
void setupVariables(typename PFP::MAP& m, LinearSolver<SOLVER_TRAITS>* s, CellMarker& lm, const AttributeHandler<ATTR_TYPE>& attr, const AttributeHandler<unsigned int> index)
void setupVariables(
typename PFP::MAP& m,
LinearSolver<SOLVER_TRAITS>* s,
const AttributeHandler<unsigned int> index,
CellMarker& lm,
const AttributeHandler<ATTR_TYPE>& attr)
{
FunctorMeshToSolver_Scalar<PFP, ATTR_TYPE, SOLVER_TRAITS> fmts(m, s, lm, attr, index) ;
FunctorMeshToSolver_Scalar<PFP, ATTR_TYPE, SOLVER_TRAITS> fmts(s, index, lm, attr) ;
m.foreach_orbit(VERTEX, fmts) ;
}
template <typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS>
void setupVariables(typename PFP::MAP& m, LinearSolver<SOLVER_TRAITS>* s, CellMarker& lm, const AttributeHandler<ATTR_TYPE>& attr, unsigned int coord, const AttributeHandler<unsigned int> index)
void setupVariables(
typename PFP::MAP& m,
LinearSolver<SOLVER_TRAITS>* s,
const AttributeHandler<unsigned int> index,
CellMarker& lm,
const AttributeHandler<ATTR_TYPE>& attr,
unsigned int coord)
{
FunctorMeshToSolver_Vector<PFP, ATTR_TYPE, SOLVER_TRAITS> fmts(m, s, lm, attr, coord, index) ;
FunctorMeshToSolver_Vector<PFP, ATTR_TYPE, SOLVER_TRAITS> fmts(s, index, lm, attr, coord) ;
m.foreach_orbit(VERTEX, fmts) ;
}
/*******************************************************************************
* START MATRIX DEFINITION
*******************************************************************************/
template <class SOLVER_TRAITS>
void startMatrix(LinearSolver<SOLVER_TRAITS>* s)
{
s->begin_system() ;
}
enum ConstraintType
/*******************************************************************************
* MATRIX SETUP : EQUALITY
*******************************************************************************/
template <typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS>
void addRowsRHS_Equality(
typename PFP::MAP& m,
LinearSolver<SOLVER_TRAITS>* s,
const AttributeHandler<unsigned int> index,
const AttributeHandler<ATTR_TYPE>& attr,
float amount)
{
FunctorEquality_Scalar<PFP, ATTR_TYPE, SOLVER_TRAITS> ec(s, index, attr, amount) ;
m.foreach_orbit(VERTEX, ec) ;
}
template <typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS>
void addRowsRHS_Equality(
typename PFP::MAP& m,
LinearSolver<SOLVER_TRAITS>* s,
const AttributeHandler<unsigned int> index,
const AttributeHandler<ATTR_TYPE>& attr,
float amount,
unsigned int coord)
{
FunctorEquality_Vector<PFP, ATTR_TYPE, SOLVER_TRAITS> ec(s, index, attr, amount, coord) ;
m.foreach_orbit(VERTEX, ec) ;
}
/*******************************************************************************
* MATRIX SETUP : LAPLACIAN TOPO
*******************************************************************************/
template <typename PFP, class SOLVER_TRAITS>
void addRows_Laplacian_Topo(
typename PFP::MAP& m,
LinearSolver<SOLVER_TRAITS>* s,
const AttributeHandler<unsigned int> index)
{
LAPLACIAN_TOPO,
LAPLACIAN_COTWEIGHT,
EQUALITY
};
FunctorLaplacianTopo<PFP, SOLVER_TRAITS> flt(m, s, index) ;
m.foreach_orbit(VERTEX, flt) ;
}
template <typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS>
void addMatrixRows(typename PFP::MAP& m, LinearSolver<SOLVER_TRAITS>* s, const AttributeHandler<ATTR_TYPE>& attr, ConstraintType ct, const AttributeHandler<unsigned int> index, float* params = NULL)
void addRowsRHS_Laplacian_Topo(
typename PFP::MAP& m,
LinearSolver<SOLVER_TRAITS>* s,
const AttributeHandler<unsigned int> index,
const AttributeHandler<ATTR_TYPE>& attr)
{
switch(ct)
{
case LAPLACIAN_TOPO :
setupLaplacianMatrix<PFP, SOLVER_TRAITS>(m, s, TOPOLOGICAL, index) ;
break ;
case LAPLACIAN_COTWEIGHT :
setupLaplacianMatrix<PFP, SOLVER_TRAITS>(m, s, COTWEIGHT, index) ;
break ;
case EQUALITY :
setupEqualityMatrix<PFP, ATTR_TYPE, SOLVER_TRAITS>(m, s, attr, index, params[0]) ;
break ;
}
FunctorLaplacianTopoRHS_Scalar<PFP, ATTR_TYPE, SOLVER_TRAITS> flt(m, s, index, attr) ;
m.foreach_orbit(VERTEX, flt) ;
}
template <typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS>
void addMatrixRows(typename PFP::MAP& m, LinearSolver<SOLVER_TRAITS>* s, const AttributeHandler<ATTR_TYPE>& attr, ConstraintType ct, unsigned int coord, const AttributeHandler<unsigned int> index, float* params = NULL)
void addRowsRHS_Laplacian_Topo(
typename PFP::MAP& m,
LinearSolver<SOLVER_TRAITS>* s,
const AttributeHandler<unsigned int> index,
const AttributeHandler<ATTR_TYPE>& attr,
unsigned int coord)
{
switch(ct)
{
case LAPLACIAN_TOPO :
setupLaplacianMatrix<PFP, SOLVER_TRAITS>(m, s, TOPOLOGICAL, index) ;
break ;
case LAPLACIAN_COTWEIGHT :
setupLaplacianMatrix<PFP, SOLVER_TRAITS>(m, s, COTWEIGHT, index) ;
break ;
case EQUALITY :
setupEqualityMatrix<PFP, ATTR_TYPE, SOLVER_TRAITS>(m, s, attr, coord, index, params[0]) ;
break ;
}
FunctorLaplacianTopoRHS_Vector<PFP, ATTR_TYPE, SOLVER_TRAITS> flt(m, s, index, attr, coord) ;
m.foreach_orbit(VERTEX, flt) ;
}
/*******************************************************************************
* MATRIX SETUP : LAPLACIAN COTAN
*******************************************************************************/