exampleHarmonicParametrization.cpp

 
 
#include <iostream>
 
#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/dec/PolygonalCalculus.h>
#include "DGtal/math/linalg/DirichletConditions.h"
 
#include <polyscope/polyscope.h>
#include <polyscope/surface_mesh.h>
#include "ConfigExamples.h"
 
#include <DGtal/io/boards/Board2D.h>
 
using namespace DGtal;
 
// Using standard 3D digital space.
typedef Shortcuts<Z3i::KSpace>         SH3;
typedef ShortcutsGeometry<Z3i::KSpace> SHG3;
// The following typedefs are useful
typedef SurfaceMesh< Z3i::RealPoint, Z3i::RealVector >  SurfMesh;
typedef SurfMesh::Vertices                    Vertices;
typedef SurfMesh::RealPoint                   RealPoint;
typedef SurfMesh::Face                        Face;
typedef SurfMesh::Vertex                      Vertex;
typedef PolygonalCalculus<SH3::RealPoint,SH3::RealVector> PC;
typedef PC::SparseMatrix SparseMatrix;
typedef PC::DenseMatrix DenseMatrix;
typedef PC::Solver Solver;
typedef PC::Vector Vector;
typedef PC::Triplet Triplet;
typedef std::vector<Vertex> chain;
typedef DirichletConditions< PC::LinAlg > Conditions;
typedef Conditions::IntegerVector IntegerVector;
 
 
//Polyscope global
polyscope::SurfaceMesh *psMesh;
polyscope::SurfaceMesh *psParam;
 
//DEC
PC *calculus;
 
//surface mesh global
SurfMesh surfmesh;
std::vector<std::vector<size_t>> faces;
std::vector<RealPoint> positions;
 
 
std::pair<Vector,Vector> FixBoundaryParametrization(const std::vector<Vertex>& boundary)
{
    auto nb = boundary.size();
    auto n = surfmesh.nbVertices();
    Vector u = Vector::Zero(n),v = Vector::Zero(n);
    double totalBoundaryLength = 0;
    for (Vertex i = 0;i<nb;i++)
        totalBoundaryLength += surfmesh.distance(boundary[(i+1)%nb],boundary[i]);
 
    double partialSum = 0;
    for (Vertex i = 0;i<nb;i++)
    {
        double th = 2*M_PI*partialSum/totalBoundaryLength;
        auto vi = boundary[i];
        auto vj = boundary[(i+1)%nb];
        u(vi) = std::cos(th);
        v(vj) = std::sin(th);
        partialSum += surfmesh.distance(vi,vj);
    }
    return {u,v};
}
 
void VisualizeParametrizationOnCircle(const DenseMatrix& UV)
{
  auto n= surfmesh.nbVertices();
  std::vector<RealPoint> pos(n);
  double scale = 0;
  RealPoint avg = {0.,0.,0.};
  for (size_t v = 0;v<n;v++)
  {
    auto p = surfmesh.position(v);
    avg += p;
    if (p.norm() > scale)
      scale = p.norm();
  }
  avg /= surfmesh.nbVertices();
  scale /= 2;
  
  for (size_t v = 0;v<n;v++)
  {
    pos[v] = RealPoint{scale*UV(v,0),scale*UV(v,1),0.} + avg;
  }
  polyscope::registerSurfaceMesh("On circle parametrization", pos, faces)->setEnabled(false);
  
  Board2D board;
  auto pairs = surfmesh.allEdgeVertices();
  for(auto p: pairs)
  {
    if ((std::pow(UV(p.first,0)-UV(p.second,0),2.0)+std::pow(UV(p.first,1)- UV(p.second,1),2.0)) > 0.0)
        board.drawLine(UV(p.first,0), UV(p.first,1), UV(p.second,0), UV(p.second,1));
  }
  board.saveSVG("layout.svg");
  
}
 
DenseMatrix HarmonicParametrization()
{
    auto n = surfmesh.nbVertices();
    std::cout<<"Nb boundary edges = "<< surfmesh.computeManifoldBoundaryEdges().size()<<std::endl;
    std::vector<chain> chains = surfmesh.computeManifoldBoundaryChains();
    //choose longest chain as boundary of the parametrization
    std::cout<<"Nb boundaries  = "<< chains.size() << std::endl;
  
    //choose longest chain as boundary of the parametrization
    auto B = *std::max_element(chains.begin(),chains.end(),[] (const chain& A,const chain& B) {return A.size() < B.size();});
 
    IntegerVector boundary = IntegerVector::Zero(n);
    for (Vertex v : B)
        boundary(v) = 1;
 
    std::pair<Vector,Vector> uv_b = FixBoundaryParametrization(B);//maps boundary to circle
 
    calculus = new PC(surfmesh);
 
    //Impose dirichlet boundary condition to laplace problem
    Vector Z = Vector::Zero(n);
    SparseMatrix L = calculus->globalLaplaceBeltrami();
    SparseMatrix L_d = Conditions::dirichletOperator( L, boundary );
 
    PC::Solver solver;
    solver.compute(L_d);
 
    Vector b_u = Conditions::dirichletVector( L, Z,boundary, uv_b.first );
    Vector b_v = Conditions::dirichletVector( L, Z,boundary, uv_b.second );
 
    Vector rslt_u_d = solver.solve(b_u);
    Vector rslt_v_d = solver.solve(b_v);
 
    Vector rslt_u = Conditions::dirichletSolution(rslt_u_d,boundary,uv_b.first);
    Vector rslt_v = Conditions::dirichletSolution(rslt_v_d,boundary,uv_b.second);
 
    DenseMatrix uv(n,2);
    uv.col(0) = rslt_u;
    uv.col(1) = rslt_v;
 
    return uv;
}
 
int main(int argc, char **argv)
{
    //Import Voxel Model
    auto params = SH3::defaultParameters() | SHG3::defaultParameters() |  SHG3::parametersGeometryEstimation();
    params("surfaceComponents", "All");
 
    //load .vol
    std::string inputFilename(examplesPath + "samples/lucy-64.vol" );
 
    auto binary_image    = SH3::makeBinaryImage(inputFilename, params );
 
    //offset K space to create boundary
    auto K               = SH3::getKSpace( binary_image, params );
    K.init(K.lowerBound()+SH3::Point(5,0,0),K.upperBound(),true);
 
    auto surface         = SH3::makeDigitalSurface( binary_image, K, params );
    auto primalSurface   = SH3::makePrimalSurfaceMesh(surface);
 
 
    //Need to convert the faces
    for(size_t 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());
 
 
    // Initialize polyscope
    polyscope::init();
 
    //compute parametrization
    DenseMatrix UV = HarmonicParametrization();
 
    //visualize parametrization on 2D circle
    VisualizeParametrizationOnCircle(UV);
 
    psMesh = polyscope::registerSurfaceMesh("Digital Surface", positions, faces);
    psMesh->addVertexParameterizationQuantity("Harmonic parametrization",UV)->setEnabled(true);
 
    //set correct view for voxel mesh
    polyscope::view::upDir = polyscope::view::UpDir::XUp;
 
    polyscope::show();
    return EXIT_SUCCESS;
}