dgtalCalculus-bunny.cpp File Reference

#include <iostream>
#include <string>
#include <DGtal/base/Common.h>
#include <DGtal/helpers/StdDefs.h>
#include <DGtal/helpers/Shortcuts.h>
#include <DGtal/helpers/ShortcutsGeometry.h>
#include <DGtal/shapes/SurfaceMesh.h>
#include <DGtal/geometry/surfaces/DigitalSurfaceRegularization.h>
#include <DGtal/dec/PolygonalCalculus.h>
#include <DGtal/math/linalg/DirichletConditions.h>
#include <polyscope/polyscope.h>
#include <polyscope/surface_mesh.h>
#include <polyscope/point_cloud.h>
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include "ConfigExamples.h"

Include dependency graph for dgtalCalculus-bunny.cpp:

Go to the source code of this file.

Typedefs
typedef Shortcuts< Z3i::KSpace >	SH3
typedef ShortcutsGeometry< Z3i::KSpace >	SHG3
typedef SurfaceMesh< RealPoint, RealVector >	SurfMesh
typedef SurfMesh::Face	Face
typedef SurfMesh::Vertex	Vertex

Functions
double	phiVertex (const Vertex v)
PolygonalCalculus< SH3::RealPoint, SH3::RealVector >::Vector	phi (const Face f)
void	initPhi ()
void	initQuantities ()
void	initQuantitiesCached ()
void	computeLaplace ()
void	myCallback ()
int	main ()

Variables
polyscope::SurfaceMesh *	psMesh
SurfMesh	surfmesh
float	scale = 0.1
PolygonalCalculus< SH3::RealPoint, SH3::RealVector >::Vector	phiEigen

Detailed Description

This program 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 3 of the License, or (at your option) any later version.

This program 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 General Public License for more details.

You should have received a copy of the GNU General Public License along with this program. If not, see http://www.gnu.org/licenses/.

Author

David Coeurjolly (david.nosp@m..coe.nosp@m.urjol.nosp@m.ly@l.nosp@m.iris..nosp@m.cnrs.nosp@m..fr ) Laboratoire d'InfoRmatique en Image et Systemes d'information - LIRIS (CNRS, UMR 5205), CNRS, France

Date

2021/09/02

This file is part of the DGtal library.

Definition in file dgtalCalculus-bunny.cpp.

Typedef Documentation

Face

typedef SurfMesh::Face Face

Definition at line 53 of file dgtalCalculus-bunny.cpp.

SH3

typedef Shortcuts<Z3i::KSpace> SH3

Definition at line 49 of file dgtalCalculus-bunny.cpp.

SHG3

typedef ShortcutsGeometry<Z3i::KSpace> SHG3

Definition at line 50 of file dgtalCalculus-bunny.cpp.

SurfMesh

typedef SurfaceMesh< RealPoint, RealVector > SurfMesh

Examples

shapes/exampleSurfaceMesh.cpp.

Definition at line 52 of file dgtalCalculus-bunny.cpp.

Vertex

typedef SurfMesh::Vertex Vertex

Definition at line 54 of file dgtalCalculus-bunny.cpp.

Function Documentation

computeLaplace()

void computeLaplace

(

)

Definition at line 167 of file dgtalCalculus-bunny.cpp.

{
  PolygonalCalculus<SH3::RealPoint,SH3::RealVector> calculus(surfmesh);
  trace.beginBlock("Operator construction...");
  PolygonalCalculus<SH3::RealPoint,SH3::RealVector>::SparseMatrix L = calculus.globalLaplaceBeltrami();
  trace.endBlock();
   
  const auto nbv = surfmesh.nbVertices();
  PolygonalCalculus<SH3::RealPoint,SH3::RealVector>::Vector g = PolygonalCalculus<SH3::RealPoint,SH3::RealVector>::Vector::Zero(nbv);
 
  //Setting some random sources
  typedef DirichletConditions< EigenLinearAlgebraBackend >  DC;
  DC::IntegerVector p = DC::nullBoundaryVector( g );
  for(auto cpt=0; cpt< 10;++cpt)
    {
      int idx  = rand() % nbv;
      g( idx ) = rand() % 100;
      p( idx ) = 1.0;
    }
  
  //Solve Δu=0 with g as boundary conditions
  PolygonalCalculus<SH3::RealPoint,SH3::RealVector>::Solver solver;
  
  trace.beginBlock("Prefactorization...");
  DC::SparseMatrix L_dirichlet = DC::dirichletOperator( L, p );
  solver.compute( L_dirichlet );  
  ASSERT(solver.info()==Eigen::Success);
  trace.endBlock();
  
  trace.beginBlock("Solve...");
  PolygonalCalculus<SH3::RealPoint,SH3::RealVector>::Vector g_dirichlet = DC::dirichletVector( L, g, p, g );
  PolygonalCalculus<SH3::RealPoint,SH3::RealVector>::Vector x_dirichlet = solver.solve( g_dirichlet );
  PolygonalCalculus<SH3::RealPoint,SH3::RealVector>::Vector u = DC::dirichletSolution( x_dirichlet, p, g );
  ASSERT(solver.info()==Eigen::Success);
  trace.endBlock();
  
  psMesh->addVertexScalarQuantity("g", g);
  psMesh->addVertexScalarQuantity("u", u);
}

References DGtal::Trace::beginBlock(), calculus, DGtal::Trace::endBlock(), DGtal::PolygonalCalculus< TRealPoint, TRealVector >::globalLaplaceBeltrami(), DGtal::SurfaceMesh< TRealPoint, TRealVector >::nbVertices(), psMesh, surfmesh, and DGtal::trace.

Referenced by myCallback().

initPhi()

void initPhi

(

)

Definition at line 83 of file dgtalCalculus-bunny.cpp.

{
  phiEigen.resize(surfmesh.nbVertices());
  for(auto i = 0; i < surfmesh.nbVertices(); ++i)
    phiEigen(i) = phiVertex(i);
  psMesh->addVertexScalarQuantity("Phi", phiEigen);
}

References DGtal::SurfaceMesh< TRealPoint, TRealVector >::nbVertices(), phiEigen, phiVertex(), psMesh, and surfmesh.

Referenced by myCallback().

initQuantities()

void initQuantities

(

)

Definition at line 91 of file dgtalCalculus-bunny.cpp.

{
  trace.beginBlock("Basic quantities");
  PolygonalCalculus<SH3::RealPoint,SH3::RealVector> calculus(surfmesh);
  
  std::vector<PolygonalCalculus<SH3::RealPoint,SH3::RealVector>::Vector> gradients;
  std::vector<PolygonalCalculus<SH3::RealPoint,SH3::RealVector>::Vector> cogradients;
  std::vector<PolygonalCalculus<SH3::RealPoint,SH3::RealVector>::Real3dPoint> normals;
  std::vector<PolygonalCalculus<SH3::RealPoint,SH3::RealVector>::Real3dPoint> vectorArea;
  std::vector<PolygonalCalculus<SH3::RealPoint,SH3::RealVector>::Real3dPoint> centroids;
  
  std::vector<double> faceArea;
  
  for(auto f=0; f < surfmesh.nbFaces(); ++f)
  {
    PolygonalCalculus<SH3::RealPoint,SH3::RealVector>::Vector grad = calculus.gradient(f) * phi(f);
    gradients.push_back( grad );
    
    PolygonalCalculus<SH3::RealPoint,SH3::RealVector>::Vector cograd =  calculus.coGradient(f) * phi(f);
    cogradients.push_back( cograd );
    
    normals.push_back(calculus.faceNormalAsDGtalVector(f));
    
    PolygonalCalculus<SH3::RealPoint,SH3::RealVector>::Vector vA = calculus.vectorArea(f);
    vectorArea.push_back({vA(0) , vA(1), vA(2)});
    
    faceArea.push_back( calculus.faceArea(f));
  }
  trace.endBlock();
  
  psMesh->addFaceVectorQuantity("Gradients", gradients);
  psMesh->addFaceVectorQuantity("co-Gradients", cogradients);
  psMesh->addFaceVectorQuantity("Normals", normals);
  psMesh->addFaceScalarQuantity("Face area", faceArea);
  psMesh->addFaceVectorQuantity("Vector area", vectorArea);
}

References DGtal::Trace::beginBlock(), calculus, DGtal::PolygonalCalculus< TRealPoint, TRealVector >::coGradient(), DGtal::Trace::endBlock(), DGtal::PolygonalCalculus< TRealPoint, TRealVector >::faceArea(), DGtal::PolygonalCalculus< TRealPoint, TRealVector >::faceNormalAsDGtalVector(), DGtal::PolygonalCalculus< TRealPoint, TRealVector >::gradient(), DGtal::SurfaceMesh< TRealPoint, TRealVector >::nbFaces(), phi(), psMesh, surfmesh, DGtal::trace, and DGtal::PolygonalCalculus< TRealPoint, TRealVector >::vectorArea().

Referenced by myCallback().

initQuantitiesCached()

void initQuantitiesCached

(

)

Definition at line 129 of file dgtalCalculus-bunny.cpp.

{
  trace.beginBlock("Basic quantities (cached)");
  PolygonalCalculus<SH3::RealPoint,SH3::RealVector> calculus(surfmesh,true);
  
  std::vector<PolygonalCalculus<SH3::RealPoint,SH3::RealVector>::Vector> gradients;
  std::vector<PolygonalCalculus<SH3::RealPoint,SH3::RealVector>::Vector> cogradients;
  std::vector<PolygonalCalculus<SH3::RealPoint,SH3::RealVector>::Real3dPoint> normals;
  std::vector<PolygonalCalculus<SH3::RealPoint,SH3::RealVector>::Real3dPoint> vectorArea;
  std::vector<PolygonalCalculus<SH3::RealPoint,SH3::RealVector>::Real3dPoint> centroids;
  
  std::vector<double> faceArea;
  
  for(auto f=0; f < surfmesh.nbFaces(); ++f)
  {
    PolygonalCalculus<SH3::RealPoint,SH3::RealVector>::Vector grad = calculus.gradient(f) * phi(f);
    gradients.push_back( grad );
    
    PolygonalCalculus<SH3::RealPoint,SH3::RealVector>::Vector cograd =  calculus.coGradient(f) * phi(f);
    cogradients.push_back( cograd );
    
    normals.push_back(calculus.faceNormalAsDGtalVector(f));
    
    PolygonalCalculus<SH3::RealPoint,SH3::RealVector>::Vector vA = calculus.vectorArea(f);
    vectorArea.push_back({vA(0) , vA(1), vA(2)});
    
    faceArea.push_back( calculus.faceArea(f));
  }
  trace.endBlock();
  
  psMesh->addFaceVectorQuantity("Gradients", gradients);
  psMesh->addFaceVectorQuantity("co-Gradients", cogradients);
  psMesh->addFaceVectorQuantity("Normals", normals);
  psMesh->addFaceScalarQuantity("Face area", faceArea);
  psMesh->addFaceVectorQuantity("Vector area", vectorArea);
}

References DGtal::Trace::beginBlock(), calculus, DGtal::PolygonalCalculus< TRealPoint, TRealVector >::coGradient(), DGtal::Trace::endBlock(), DGtal::PolygonalCalculus< TRealPoint, TRealVector >::faceArea(), DGtal::PolygonalCalculus< TRealPoint, TRealVector >::faceNormalAsDGtalVector(), DGtal::PolygonalCalculus< TRealPoint, TRealVector >::gradient(), DGtal::SurfaceMesh< TRealPoint, TRealVector >::nbFaces(), phi(), psMesh, surfmesh, DGtal::trace, and DGtal::PolygonalCalculus< TRealPoint, TRealVector >::vectorArea().

Referenced by myCallback().

main()

int main

(

void

)

Definition at line 224 of file dgtalCalculus-bunny.cpp.

{
  auto params = SH3::defaultParameters() | SHG3::defaultParameters() |  SHG3::parametersGeometryEstimation();
  params("surfaceComponents", "All");
 
  std::string filename = examplesPath + std::string("/samples/bunny-32.vol");
  auto binary_image    = SH3::makeBinaryImage(filename, params );
  auto K               = SH3::getKSpace( binary_image, params );
  auto surface         = SH3::makeDigitalSurface( binary_image, K, params );
  SH3::Cell2Index c2i;
  auto primalSurface   = SH3::makePrimalSurfaceMesh(surface);
  
  //Need to convert the faces
  std::vector<std::vector<SH3::SurfaceMesh::Vertex>> faces;
  std::vector<RealPoint> positions;
  
  for(auto face= 0 ; face < primalSurface->nbFaces(); ++face)
    faces.push_back(primalSurface->incidentVertices( face ));
  
  //Recasting to vector of vertices
  positions = primalSurface->positions();
   
  surfmesh = SurfMesh(positions.begin(),
                      positions.end(),
                      faces.begin(),
                      faces.end());
 
  std::cout<<"number of non-manifold Edges = " << surfmesh.computeNonManifoldEdges().size()<<std::endl;
  psMesh = polyscope::registerSurfaceMesh("digital surface", positions, faces);
 
  // Initialize polyscope
  polyscope::init();
  
  // Set the callback function
  polyscope::state::userCallback = myCallback;
  polyscope::show();
  return EXIT_SUCCESS;
}

References binary_image, DGtal::SurfaceMesh< TRealPoint, TRealVector >::computeNonManifoldEdges(), DGtal::Shortcuts< TKSpace >::defaultParameters(), DGtal::ShortcutsGeometry< TKSpace >::defaultParameters(), DGtal::Shortcuts< TKSpace >::getKSpace(), K, DGtal::Shortcuts< TKSpace >::makeBinaryImage(), DGtal::Shortcuts< TKSpace >::makeDigitalSurface(), DGtal::Shortcuts< TKSpace >::makePrimalSurfaceMesh(), myCallback(), DGtal::ShortcutsGeometry< TKSpace >::parametersGeometryEstimation(), psMesh, surface, and surfmesh.

myCallback()

void myCallback

(

)

Definition at line 207 of file dgtalCalculus-bunny.cpp.

{
  ImGui::SliderFloat("Phi scale", &scale, 0., 1.);
  if (ImGui::Button("Phi and basic operators"))
  {
    initPhi();
    initQuantities();
  }
  if (ImGui::Button("Phi and basic operators (cached)"))
  {
    initPhi();
    initQuantitiesCached();
  }
  if(ImGui::Button("Compute Laplace problem"))
    computeLaplace();
}

References computeLaplace(), initPhi(), initQuantities(), and initQuantitiesCached().

Referenced by main().

phi()

PolygonalCalculus< SH3::RealPoint, SH3::RealVector >::Vector phi

(

const Face

f

)

Examples

dec/exampleDECSurface.cpp, and geometry/volumes/checkFullConvexityTheorems.cpp.

Definition at line 69 of file dgtalCalculus-bunny.cpp.

{
  auto vertices = surfmesh.incidentVertices(f);
  auto nf = vertices.size();
  Eigen::VectorXd ph(nf);
  size_t cpt=0;
  for(auto v: vertices)
  {
    ph(cpt) =  phiVertex(v);
    ++cpt;
  }
  return  ph;
}

References DGtal::SurfaceMesh< TRealPoint, TRealVector >::incidentVertices(), phiVertex(), and surfmesh.

Referenced by checkFullConvexityCharacterization(), checkProjectionFullConvexity(), initQuantities(), initQuantitiesCached(), DGtal::functors::SphericalHoughNormalVectorEstimator< TSurfel, TEmbedder >::randomRotation(), TEST_CASE(), and TEST_CASE().

phiVertex()

double phiVertex

(

const Vertex

v

)

Definition at line 63 of file dgtalCalculus-bunny.cpp.

{
  return  cos(scale*(surfmesh.position(v)[0]))*sin(scale*surfmesh.position(v)[1]);
}

References DGtal::SurfaceMesh< TRealPoint, TRealVector >::position(), and surfmesh.

Referenced by initPhi(), and phi().

Variable Documentation

phiEigen

PolygonalCalculus<SH3::RealPoint,SH3::RealVector>::Vector phiEigen

Definition at line 60 of file dgtalCalculus-bunny.cpp.

Referenced by initPhi().

psMesh

polyscope::SurfaceMesh* psMesh

Definition at line 57 of file dgtalCalculus-bunny.cpp.

Referenced by computeLaplace(), initPhi(), initQuantities(), initQuantitiesCached(), and main().

scale

float scale = 0.1

Definition at line 59 of file dgtalCalculus-bunny.cpp.

surfmesh

SurfMesh surfmesh

Examples

geometry/volumes/digitalPolyhedronBuilder3D.cpp, geometry/volumes/fullConvexitySphereGeodesics.cpp, and geometry/volumes/standardDigitalPolyhedronBuilder3D.cpp.

Definition at line 58 of file dgtalCalculus-bunny.cpp.

Referenced by DGtal::GeodesicsInHeat< TPolygonalCalculus >::compute(), DGtal::VectorsInHeat< TPolygonalCalculus >::compute(), computeLaplace(), DGtal::functors::EmbedderFromNormalVectors< TRealPoint, TRealVector >::EmbedderFromNormalVectors(), DGtal::GeodesicsInHeat< TPolygonalCalculus >::init(), DGtal::VectorsInHeat< TPolygonalCalculus >::init(), initPhi(), initQuantities(), initQuantitiesCached(), main(), main(), phi(), phiVertex(), saveToObj(), and TEST_CASE().