Commit c0fb7beb authored by Pierre Kraemer's avatar Pierre Kraemer
Browse files

generic laplacian computation + cotangent weights

parent 3ea0e63a
...@@ -103,7 +103,7 @@ typename PFP::REAL vertexVoronoiArea(typename PFP::MAP& map, Dart d, const typen ...@@ -103,7 +103,7 @@ typename PFP::REAL vertexVoronoiArea(typename PFP::MAP& map, Dart d, const typen
Dart it = d ; Dart it = d ;
do do
{ {
area += convexFaceArea<PFP>(map, it, position) ; area += convexFaceArea<PFP>(map, it, position) / 3 ;
it = map.alpha1(it) ; it = map.alpha1(it) ;
} while(it != d) ; } while(it != d) ;
return area ; return area ;
......
...@@ -40,34 +40,34 @@ namespace Algo ...@@ -40,34 +40,34 @@ namespace Algo
namespace Geometry namespace Geometry
{ {
template <typename PFP> template <typename PFP, typename ATTR_TYPE>
typename PFP::VEC3 computeLaplacianTopoVertex( ATTR_TYPE computeLaplacianTopoVertex(
typename PFP::MAP& map, typename PFP::MAP& map,
Dart d, Dart d,
const typename PFP::TVEC3& position) ; const AttributeHandler<ATTR_TYPE>& attr) ;
template <typename PFP> template <typename PFP, typename ATTR_TYPE>
typename PFP::VEC3 computeLaplacianCotanVertex( ATTR_TYPE computeLaplacianCotanVertex(
typename PFP::MAP& map, typename PFP::MAP& map,
Dart d, Dart d,
const typename PFP::TVEC3& position,
const typename PFP::TREAL& edgeWeight, const typename PFP::TREAL& edgeWeight,
const typename PFP::TREAL& vertexArea) ; const typename PFP::TREAL& vertexArea,
const AttributeHandler<ATTR_TYPE>& attr) ;
template <typename PFP> template <typename PFP, typename ATTR_TYPE>
void computeLaplacianTopoVertices( void computeLaplacianTopoVertices(
typename PFP::MAP& map, typename PFP::MAP& map,
const typename PFP::TVEC3& position, const AttributeHandler<ATTR_TYPE>& attr,
typename PFP::TVEC3& laplacian, AttributeHandler<ATTR_TYPE>& laplacian,
const FunctorSelect& select = SelectorTrue()) ; const FunctorSelect& select = SelectorTrue()) ;
template <typename PFP> template <typename PFP, typename ATTR_TYPE>
void computeLaplacianCotanVertices( void computeLaplacianCotanVertices(
typename PFP::MAP& map, typename PFP::MAP& map,
const typename PFP::TVEC3& position,
const typename PFP::TREAL& edgeWeight, const typename PFP::TREAL& edgeWeight,
const typename PFP::TREAL& vertexArea, const typename PFP::TREAL& vertexArea,
typename PFP::TVEC3& laplacian, const AttributeHandler<ATTR_TYPE>& attr,
AttributeHandler<ATTR_TYPE>& laplacian,
const FunctorSelect& select = SelectorTrue()) ; const FunctorSelect& select = SelectorTrue()) ;
template <typename PFP> template <typename PFP>
......
...@@ -33,45 +33,42 @@ namespace Algo ...@@ -33,45 +33,42 @@ namespace Algo
namespace Geometry namespace Geometry
{ {
template <typename PFP> template <typename PFP, typename ATTR_TYPE>
typename PFP::VEC3 computeLaplacianTopoVertex( ATTR_TYPE computeLaplacianTopoVertex(
typename PFP::MAP& map, typename PFP::MAP& map,
Dart d, Dart d,
const typename PFP::TVEC3& position) const AttributeHandler<ATTR_TYPE>& attr)
{ {
typedef typename PFP::VEC3 VEC3 ; ATTR_TYPE l(0) ;
VEC3 l(0) ;
unsigned int val = 0 ; unsigned int val = 0 ;
Dart dd = d ; Dart it = d ;
do do
{ {
l += vectorOutOfDart<PFP>(map, dd, position) ; l += attr[map.phi1(it)] - attr[it] ;
val++ ; val++ ;
dd = map.alpha1(dd) ; it = map.alpha1(it) ;
} while(dd != d) ; } while(it != d) ;
l /= val ; l /= val ;
return l ; return l ;
} }
template <typename PFP> template <typename PFP, typename ATTR_TYPE>
typename PFP::VEC3 computeLaplacianCotanVertex( ATTR_TYPE computeLaplacianCotanVertex(
typename PFP::MAP& map, typename PFP::MAP& map,
Dart d, Dart d,
const typename PFP::TVEC3& position,
const typename PFP::TREAL& edgeWeight, const typename PFP::TREAL& edgeWeight,
const typename PFP::TREAL& vertexArea) const typename PFP::TREAL& vertexArea,
const AttributeHandler<ATTR_TYPE>& attr)
{ {
typedef typename PFP::VEC3 VEC3 ;
typedef typename PFP::REAL REAL ; typedef typename PFP::REAL REAL ;
VEC3 l(0) ; ATTR_TYPE l(0) ;
Dart it = d ; Dart it = d ;
REAL vArea = vertexArea[d] ; REAL vArea = vertexArea[d] ;
REAL val = 0 ; REAL val = 0 ;
do do
{ {
REAL w = edgeWeight[it] / vArea ; REAL w = edgeWeight[it] / vArea ;
VEC3 v = vectorOutOfDart<PFP>(map, it, position) * w ; l += (attr[map.phi1(it)] - attr[it]) * w ;
l += v ;
val += w ; val += w ;
it = map.alpha1(it) ; it = map.alpha1(it) ;
} while(it != d) ; } while(it != d) ;
...@@ -79,12 +76,12 @@ typename PFP::VEC3 computeLaplacianCotanVertex( ...@@ -79,12 +76,12 @@ typename PFP::VEC3 computeLaplacianCotanVertex(
return l ; return l ;
} }
template <typename PFP> template <typename PFP, typename ATTR_TYPE>
void computeLaplacianTopoVertices( void computeLaplacianTopoVertices(
typename PFP::MAP& map, typename PFP::MAP& map,
const typename PFP::TVEC3& position, const AttributeHandler<ATTR_TYPE>& attr,
typename PFP::TVEC3& laplacian, AttributeHandler<ATTR_TYPE>& laplacian,
const FunctorSelect& select) const FunctorSelect& select)
{ {
CellMarker marker(map, VERTEX); CellMarker marker(map, VERTEX);
for(Dart d = map.begin(); d != map.end(); map.next(d)) for(Dart d = map.begin(); d != map.end(); map.next(d))
...@@ -92,19 +89,19 @@ void computeLaplacianTopoVertices( ...@@ -92,19 +89,19 @@ void computeLaplacianTopoVertices(
if(select(d) && !marker.isMarked(d)) if(select(d) && !marker.isMarked(d))
{ {
marker.mark(d); marker.mark(d);
laplacian[d] = computeLaplacianTopoVertex<PFP>(map, d, position) ; laplacian[d] = computeLaplacianTopoVertex<PFP, ATTR_TYPE>(map, d, attr) ;
} }
} }
} }
template <typename PFP> template <typename PFP, typename ATTR_TYPE>
void computeLaplacianCotanVertices( void computeLaplacianCotanVertices(
typename PFP::MAP& map, typename PFP::MAP& map,
const typename PFP::TVEC3& position, const typename PFP::TREAL& edgeWeight,
const typename PFP::TREAL& edgeWeight, const typename PFP::TREAL& vertexArea,
const typename PFP::TREAL& vertexArea, const AttributeHandler<ATTR_TYPE>& attr,
typename PFP::TVEC3& laplacian, AttributeHandler<ATTR_TYPE>& laplacian,
const FunctorSelect& select) const FunctorSelect& select)
{ {
CellMarker marker(map, VERTEX); CellMarker marker(map, VERTEX);
for(Dart d = map.begin(); d != map.end(); map.next(d)) for(Dart d = map.begin(); d != map.end(); map.next(d))
...@@ -112,7 +109,7 @@ void computeLaplacianCotanVertices( ...@@ -112,7 +109,7 @@ void computeLaplacianCotanVertices(
if(select(d) && !marker.isMarked(d)) if(select(d) && !marker.isMarked(d))
{ {
marker.mark(d); marker.mark(d);
laplacian[d] = computeLaplacianCotanVertex<PFP>(map, d, position, edgeWeight, vertexArea) ; laplacian[d] = computeLaplacianCotanVertex<PFP, ATTR_TYPE>(map, d, edgeWeight, vertexArea, attr) ;
} }
} }
} }
......
...@@ -26,6 +26,7 @@ ...@@ -26,6 +26,7 @@
#define __LINEAR_SOLVING_BASIC__ #define __LINEAR_SOLVING_BASIC__
#include "OpenNL/linear_solver.h" #include "OpenNL/linear_solver.h"
#include "Algo/LinearSolving/variablesSetup.h"
#include "Algo/LinearSolving/matrixSetup.h" #include "Algo/LinearSolving/matrixSetup.h"
namespace CGoGN namespace CGoGN
...@@ -34,112 +35,9 @@ namespace CGoGN ...@@ -34,112 +35,9 @@ namespace CGoGN
namespace LinearSolving namespace LinearSolving
{ {
template <typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS> /*******************************************************************************
class FunctorMeshToSolver_Scalar : public FunctorMap<typename PFP::MAP> * SOLVER INIT
{ *******************************************************************************/
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 ;
}
} ;
template <class SOLVER_TRAITS> template <class SOLVER_TRAITS>
void initSolver(LinearSolver<SOLVER_TRAITS>* s, unsigned int nb_variables, bool least_squares, bool direct) void initSolver(LinearSolver<SOLVER_TRAITS>* s, unsigned int nb_variables, bool least_squares, bool direct)
...@@ -159,81 +57,208 @@ void initSolver(LinearSolver<SOLVER_TRAITS>* s, unsigned int nb_variables, bool ...@@ -159,81 +57,208 @@ void initSolver(LinearSolver<SOLVER_TRAITS>* s, unsigned int nb_variables, bool
s->set_direct(direct) ; s->set_direct(direct) ;
} }
/*******************************************************************************
* VARIABLES SETUP
*******************************************************************************/
template <typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS> 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) ; m.foreach_orbit(VERTEX, fmts) ;
} }
template <typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS> 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) ; m.foreach_orbit(VERTEX, fmts) ;
} }
/*******************************************************************************
* START MATRIX DEFINITION
*******************************************************************************/
template <class SOLVER_TRAITS> template <class SOLVER_TRAITS>
void startMatrix(LinearSolver<SOLVER_TRAITS>* s) void startMatrix(LinearSolver<SOLVER_TRAITS>* s)
{ {
s->begin_system() ; s->begin_system() ;
} }
/*******************************************************************************
* 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> template <typename PFP, class SOLVER_TRAITS>
void addRows_Laplacian_Topo(typename PFP::MAP& m, LinearSolver<SOLVER_TRAITS>* s, const AttributeHandler<unsigned int> index) void addRows_Laplacian_Topo(
typename PFP::MAP& m,
LinearSolver<SOLVER_TRAITS>* s,
const AttributeHandler<unsigned int> index)
{
FunctorLaplacianTopo<PFP, SOLVER_TRAITS> flt(m, s, index) ;
m.foreach_orbit(VERTEX, flt) ;
}
template <typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS>
void addRowsRHS_Laplacian_Topo(
typename PFP::MAP& m,
LinearSolver<SOLVER_TRAITS>* s,
const AttributeHandler<unsigned int> index,
const AttributeHandler<ATTR_TYPE>& attr)
{ {
FunctorLaplacianTopo<PFP, SOLVER_TRAITS> lt(m, s, index) ; FunctorLaplacianTopoRHS_Scalar<PFP, ATTR_TYPE, SOLVER_TRAITS> flt(m, s, index, attr) ;
m.foreach_orbit(VERTEX, lt) ; m.foreach_orbit(VERTEX, flt) ;
} }
template <typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS>
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)
{
FunctorLaplacianTopoRHS_Vector<PFP, ATTR_TYPE, SOLVER_TRAITS> flt(m, s, index, attr, coord) ;
m.foreach_orbit(VERTEX, flt) ;
}
/*******************************************************************************
* MATRIX SETUP : LAPLACIAN COTAN
*******************************************************************************/
template <typename PFP, class SOLVER_TRAITS> template <typename PFP, class SOLVER_TRAITS>
void addRows_Laplacian_Cotan(typename PFP::MAP& m, LinearSolver<SOLVER_TRAITS>* s, const AttributeHandler<unsigned int> index, const typename PFP::TREAL& edgeWeight, const typename PFP::TREAL& vertexArea) void addRows_Laplacian_Cotan(
typename PFP::MAP& m,
LinearSolver<SOLVER_TRAITS>* s,
const AttributeHandler<unsigned int> index,
const typename PFP::TREAL& edgeWeight,
const typename PFP::TREAL& vertexArea)
{ {
FunctorLaplacianCotan<PFP, SOLVER_TRAITS> lt(m, s, index, edgeWeight, vertexArea) ; FunctorLaplacianCotan<PFP, SOLVER_TRAITS> flc(m, s, index, edgeWeight, vertexArea) ;
m.foreach_orbit(VERTEX, lt) ; m.foreach_orbit(VERTEX, flc) ;
} }
template <typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS> template <typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS>
void addRows_Equality(typename PFP::MAP& m, LinearSolver<SOLVER_TRAITS>* s, const AttributeHandler<ATTR_TYPE>& attr, const AttributeHandler<unsigned int> index, float amount) void addRowsRHS_Laplacian_Cotan(
typename PFP::MAP& m,
LinearSolver<SOLVER_TRAITS>* s,
const AttributeHandler<unsigned int> index,
const typename PFP::TREAL& edgeWeight,
const typename PFP::TREAL& vertexArea,
const AttributeHandler<ATTR_TYPE>& attr)
{ {
setupEqualityMatrix<PFP, ATTR_TYPE, SOLVER_TRAITS>(m, s, attr, index, amount) ; FunctorLaplacianCotanRHS_Scalar<PFP, ATTR_TYPE, SOLVER_TRAITS> flc(m, s, index, edgeWeight, vertexArea, attr) ;
m.foreach_orbit(VERTEX, flc) ;
} }
template <typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS> template <typename PFP, typename ATTR_TYPE, class SOLVER_TRAITS>
void addRows_Equality(typename PFP::MAP& m, LinearSolver<SOLVER_TRAITS>* s, const AttributeHandler<ATTR_TYPE>& attr, unsigned int coord, const AttributeHandler<unsigned int> index, float amount) void addRowsRHS_Laplacian_Cotan(
typename PFP::MAP& m,
LinearSolver<SOLVER_TRAITS>* s,
const AttributeHandler<unsigned int> index,
const typename PFP::TREAL& edgeWeight,
const typename PFP::TREAL& vertexArea,
const AttributeHandler<ATTR_TYPE>& attr,
unsigned int coord)
{ {
setupEqualityMatrix<PFP, ATTR_TYPE, SOLVER_TRAITS>(m, s, attr, coord, index, amount) ; FunctorLaplacianCotanRHS_Vector<PFP, ATTR_TYPE, SOLVER_TRAITS> flc(m, s, index, edgeWeight, vertexArea, attr, coord) ;
m.foreach_orbit(VERTEX, flc) ;
} }
/*******************************************************************************
* END MATRIX DEFINITION
*******************************************************************************/
template <class SOLVER_TRAITS> template <class SOLVER_TRAITS>
void endMatrix(LinearSolver<SOLVER_TRAITS>* s) void endMatrix(LinearSolver<SOLVER_TRAITS>* s)
{ {
s->end_system() ; s->end_system() ;
} }
/*******************************************************************************
* SOLVE SYSTEM
*******************************************************************************/
template <class SOLVER_TRAITS> template <class SOLVER_TRAITS>
void solve(LinearSolver<SOLVER_TRAITS>* s) void solve(LinearSolver<SOLVER_TRAITS>* s)
{ {
s->solve() ; s->solve() ;
} }
/*******************************************************************************
* RESET SOLVER
*******************************************************************************/
template <class SOLVER_TRAITS> template <class SOLVER_TRAITS>
void resetSolver(LinearSolver<SOLVER_TRAITS>* s, bool keepMatrix) void resetSolver(LinearSolver<SOLVER_TRAITS>* s, bool keepMatrix)
{