From dfec257e74298096a8f480890da6c9b5d62ea33f Mon Sep 17 00:00:00 2001 From: Lionel Untereiner Date: Thu, 14 Feb 2013 15:19:51 +0100 Subject: [PATCH] laplacian smoothing for sqrt3 volumique --- include/Algo/Geometry/laplacian.h | 20 +++ include/Algo/Geometry/laplacian.hpp | 40 +++++ include/Algo/Modelisation/subdivision3.hpp | 177 ++++++++++++++++++++- 3 files changed, 233 insertions(+), 4 deletions(-) diff --git a/include/Algo/Geometry/laplacian.h b/include/Algo/Geometry/laplacian.h index a1d779d3..726f1ac7 100644 --- a/include/Algo/Geometry/laplacian.h +++ b/include/Algo/Geometry/laplacian.h @@ -86,6 +86,26 @@ void computeCotanWeightEdges( } // namespace Surface +namespace Volume +{ + +namespace Geometry +{ + +template +ATTR_TYPE computeLaplacianTopoVertex( typename PFP::MAP& map, Dart d, const VertexAttribute& attr) ; + +template +void computeLaplacianTopoVertices(typename PFP::MAP& map, const VertexAttribute& attr, + VertexAttribute& laplacian, const FunctorSelect& select = allDarts) ; + + + +} // namespace Geometry + +} // namespace Volume + + } // namespace Algo } // namespace CGoGN diff --git a/include/Algo/Geometry/laplacian.hpp b/include/Algo/Geometry/laplacian.hpp index b0538bb8..89351b23 100644 --- a/include/Algo/Geometry/laplacian.hpp +++ b/include/Algo/Geometry/laplacian.hpp @@ -154,6 +154,46 @@ void computeCotanWeightEdges( } // namespace Surface + +namespace Volume +{ + +namespace Geometry +{ + +template +ATTR_TYPE computeLaplacianTopoVertex(typename PFP::MAP& map, Dart d, const VertexAttribute& attr) +{ + ATTR_TYPE l(0) ; + ATTR_TYPE value = attr[d] ; + unsigned int wSum = 0 ; + + Traversor3VE t(map, d) ; + for(Dart it = t.begin(); it != t.end(); it = t.next()) + { + l += attr[map.phi1(it)] - value ; + ++wSum ; + } + + l /= wSum ; + return l ; +} + +template +void computeLaplacianTopoVertices(typename PFP::MAP& map, const VertexAttribute& attr, + VertexAttribute& laplacian, const FunctorSelect& select) +{ + TraversorV t(map, select) ; + for(Dart d = t.begin(); d != t.end(); d = t.next()) + laplacian[d] = computeLaplacianTopoVertex(map, d, attr) ; +} + + +} // namespace Geometry + +} // namespace Volume + + } // namespace Algo } // namespace CGoGN diff --git a/include/Algo/Modelisation/subdivision3.hpp b/include/Algo/Modelisation/subdivision3.hpp index 8939a52f..148f891d 100644 --- a/include/Algo/Modelisation/subdivision3.hpp +++ b/include/Algo/Modelisation/subdivision3.hpp @@ -26,6 +26,9 @@ #include "Algo/Modelisation/subdivision.h" #include "Algo/Modelisation/extrusion.h" #include "Geometry/intersection.h" +#include "NL/nl.h" +#include "Algo/LinearSolving/basic.h" +#include "Algo/Geometry/laplacian.h" namespace CGoGN { @@ -631,12 +634,35 @@ void catmullClarkVol(typename PFP::MAP& map, EMBV& attributs, const FunctorSelec // } } +inline double sqrt3_K(unsigned int n) +{ + switch(n) + { + case 1: return 0.333333 ; + case 2: return 0.555556 ; + case 3: return 0.5 ; + case 4: return 0.444444 ; + case 5: return 0.410109 ; + case 6: return 0.388889 ; + case 7: return 0.375168 ; + case 8: return 0.365877 ; + case 9: return 0.359328 ; + case 10: return 0.354554 ; + case 11: return 0.350972 ; + case 12: return 0.348219 ; + default: + double t = cos((2.0*M_PI)/double(n)) ; + return (4.0 - t) / 9.0 ; + } +} template void sqrt3Vol(typename PFP::MAP& map, VertexAttribute& position, const FunctorSelect& selected) { DartMarkerStore m(map); + DartMarkerStore newBoundaryV(map); + // // 1-4 flip of all tetrahedra // @@ -697,6 +723,8 @@ void sqrt3Vol(typename PFP::MAP& map, VertexAttribute& posit Dart dres = Volume::Modelisation::Tetrahedralization::flip1To3(map, dit); position[dres] = faceCenter; + + newBoundaryV.markOrbit(dres); } } @@ -714,6 +742,147 @@ void sqrt3Vol(typename PFP::MAP& map, VertexAttribute& posit } } + + TraversorV tVg(map,selected); + for(Dart dit = tVg.begin() ; dit != tVg.end() ; dit = tVg.next()) + { + if(map.isBoundaryVertex(dit) && !newBoundaryV.isMarked(dit)) + { + typename PFP::VEC3 P = position[dit] ; + typename PFP::VEC3 newP(0) ; + unsigned int val = 0 ; + Dart vit = dit ; + do + { + newP += position[map.phi_1(map.phi2(map.phi1(vit)))] ; + ++val ; + vit = map.phi2(map.phi_1(vit)) ; + } while(vit != dit) ; + typename PFP::REAL K = sqrt3_K(val) ; + newP *= typename PFP::REAL(3) ; + newP -= typename PFP::REAL(val) * P ; + newP *= K / typename PFP::REAL(2 * val) ; + newP += (typename PFP::REAL(1) - K) * P ; + position[dit] = newP ; + } + } + + //AutoVertexAttribute laplacian qui est une copie de position + +// TraversorV tVg2(map,selected); +// for(Dart dit = tVg2.begin() ; dit != tVg2.end() ; dit = tVg2.next()) +// { +// if(!map.isBoundaryVertex(dit)) +// { +// //modifie laplacian +// } +// } + + //echange lapaclian et position + + + VertexAutoAttribute diffCoord(map, "diffCoord"); + Algo::Volume::Geometry::computeLaplacianTopoVertices(map, position, diffCoord) ; + + VertexAutoAttribute vIndex(map, "vIndex"); + + unsigned int nb_vertices = map.template computeIndexCells(vIndex); + + + CellMarker lockingMarker(map); + + TraversorV tv(map); + for(Dart dit = tv.begin() ; dit != tv.end() ; dit = tv.next()) + { + if(!lockingMarker.isMarked(dit) && map.isBoundaryVertex(dit)) + lockingMarker.mark(dit); + } + + + NLContext nlContext = nlNewContext(); + nlSolverParameteri(NL_NB_VARIABLES, nb_vertices); + nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE); + nlSolverParameteri(NL_SOLVER, NL_CHOLMOD_EXT); + + + nlMakeCurrent(nlContext); + if(nlGetCurrentState() == NL_STATE_INITIAL) + nlBegin(NL_SYSTEM) ; + + for(int coord = 0; coord < 3; ++coord) + { + LinearSolving::setupVariables(map, vIndex, lockingMarker, position, coord) ; + nlBegin(NL_MATRIX) ; + LinearSolving::addRowsRHS_Laplacian_Topo(map, vIndex, diffCoord, coord) ; +// LinearSolving::addRowsRHS_Laplacian_Cotan(*map, perMap->vIndex, perMap->edgeWeight, perMap->vertexArea, perMap->diffCoord, coord) ; + nlEnd(NL_MATRIX) ; + nlEnd(NL_SYSTEM) ; + nlSolve() ; + LinearSolving::getResult(map, vIndex, position, coord) ; + nlReset(NL_TRUE) ; + } + + nlDeleteContext(nlContext); + + +// +// float weight = 1.0; + +// LinearSolving::LinearSolver ls(nbV); +// ls.set_least_squares(true); + +// for(unsigned int coord = 0 ; coord < 3 ; ++coord) +// { +// std::cout << "coord " << coord << std::flush; + +// TraversorV tv(map); +// for(Dart dit = tv.begin() ; dit != tv.end() ; dit = tv.next()) +// { +// ls.variable(indexV[dit]).set_value(position[dit][coord]); + +// if(map.isBoundaryVertex(dit)) +// ls.variable(indexV[dit]).lock(); +// } +// std::cout << "... variables set... " << std::flush; + +// ls.begin_system(); + +// TraversorV tv2(map); +// for(Dart dit = tv2.begin() ; dit != tv2.end() ; dit = tv2.next()) +// { +// if(!map.isBoundaryVertex(dit)) +// { +// ls.begin_row(); + +// float sum = 0; +// Traversor3VVaE tvvae(map, dit); +// for(Dart ditvvae = tvvae.begin() ; ditvvae != tvvae.end() ; ditvvae = tvvae.next()) +// { +// ls.add_coefficient(indexV[ditvvae],weight); +// sum += weight; +// } + +// ls.add_coefficient(indexV[dit],-sum); +// ls.normalize_row(); + +// ls.end_row(); +// } +// } + +// ls.end_system(); +// std::cout << "... system built... " << std::flush; +// ls.solve(); +// std::cout << "... system solved... " << std::flush; + +// TraversorV tv3(map); +// for(Dart dit = tv3.begin() ; dit != tv3.end() ; dit = tv3.next()) +// { +// position[dit][coord] = ls.variable(indexV[dit]).value(); +// } + +// ls.reset(false); +// std::cout << "... done" << std::endl; +// } } template @@ -733,11 +902,11 @@ void computeDual(typename PFP::MAP& map, VertexAttribute& po } -} //namespace Modelisation +} // namespace Modelisation -} // volume +} // namespace volume -} //namespace Algo +} // namespace Algo -} //namespace CGoGN +} // namespace CGoGN -- GitLab