34#include <DGtal/helpers/StdDefs.h>
36#include <DGtal/base/BasicFunctors.h>
38#include <DGtal/topology/DigitalSurface.h>
39#include <DGtal/topology/SetOfSurfels.h>
40#include <DGtal/topology/CanonicCellEmbedder.h>
41#include <DGtal/topology/LightImplicitDigitalSurface.h>
43#include <DGtal/io/colormaps/ColorBrightnessColorMap.h>
45#include <DGtal/geometry/surfaces/estimation/LocalEstimatorFromSurfelFunctorAdapter.h>
46#include <DGtal/geometry/surfaces/estimation/IIGeometricFunctors.h>
47#include <DGtal/geometry/surfaces/estimation/IntegralInvariantCovarianceEstimator.h>
49#include <DGtal/shapes/parametric/Ball3D.h>
50#include <DGtal/shapes/Shapes.h>
51#include <DGtal/shapes/GaussDigitizer.h>
52#include <DGtal/shapes/Mesh.h>
53#include <DGtal/shapes/MeshHelpers.h>
54#include "DGtal/shapes/TriangulatedSurface.h"
56#include <DGtal/io/viewers/PolyscopeViewer.h>
58#include <DGtal/math/linalg/EigenSupport.h>
86template <
typename Shape>
88 std::function<
double(
const RealPoint3D&)> input_function,
89 std::function<
double(
const RealPoint3D&)> target_function,
90 int argc,
char** argv)
92 trace.beginBlock(
"Laplacian");
95 trace.beginBlock(
"Digitizing the shape");
97 digitizer.attach(shape);
103 trace.beginBlock(
"Construct the Khalimsky space from the image domain." );
105 bool space_ok = kspace.
init(
domain.lowerBound(),
106 domain.upperBound(),
true );
109 trace.error() <<
"Error in the Khamisky space construction."<<std::endl;
114 trace.beginBlock(
"Extracting boundary by scanning the space. " );
117 MySurfelAdjacency surfAdj(
true );
122 MySetOfSurfels theSetOfSurfels( kspace, surfAdj );
128 trace.info() <<
"Digital surface has " << digSurf.
size() <<
" surfels."
132 trace.beginBlock(
"Making triangulated surface. " );
135 typedef std::map< MyDigitalSurface::Vertex, TriMesh::Index > VertexMap;
141 ( digSurf, canonicCellembedder, trimesh, vmap );
143 trace.info() <<
"Triangulated surface is " << trimesh << std::endl;
146 trace.beginBlock(
"Computing Laplacian");
151 for(
auto v : vertices)
158 std::ofstream error_out(options.
error_output, std::ofstream::app);
159 std::ofstream function_out(
"function.dat");
161 Eigen::VectorXd laplacian_result( trimesh.
nbVertices() );
162 Eigen::VectorXd error( trimesh.
nbVertices() );
163 Eigen::VectorXd error_faces( trimesh.
nbFaces() );
165 double total_area = 0.;
168 for(
auto v : vertices)
175 for(
auto a : out_arcs)
194 const RealPoint3D v1 = (p1 - p2) / (p1 - p2).norm();
195 const RealPoint3D v2 = (p3 - p2) / (p3 - p2).norm();
196 const RealPoint3D vv1 = (pp1 - pp2) / (pp1 - pp2).norm();
197 const RealPoint3D vv2 = (pp3 - pp2) / (pp3 - pp2).norm();
199 const double alpha = acos( v1.
dot(v2) );
200 const double beta = acos( vv1.
dot(vv2) );
203 accum += (tan(M_PI_2 - alpha) + tan(M_PI_2 - beta)) * (input_function(p_j) - input_function(p_i));
206 double accum_area = 0.;
208 for(
auto f : faces_around)
216 const double faceArea = .5 * cross.
norm();
218 accum_area += faceArea / 3.;
221 total_area += accum_area;
224 ? laplacian_result(i) = (1 / (2. * accum_area)) * accum
225 : laplacian_result(i) = .5 * accum;
228 const double real_laplacian_value = target_function(w_projected);
232 function_out << w_s[1] <<
" "
234 << laplacian_result(i) <<
" "
235 << real_laplacian_value <<
" "
236 << input_function(p_i) << std::endl;
238 error(i) = laplacian_result(i) - real_laplacian_value;
239 for(
auto f : faces_around)
240 error_faces( f ) += error(i) / faces_around.size();
247 trace.info() <<
"Computed Area / Real Area : " << total_area <<
" " << 4 * M_PI << std::endl;
248 trace.info() <<
"Mean Error / Max Error : "
249 << error.array().abs().mean() <<
" / " << error.array().abs().maxCoeff(&max_index) << std::endl;
250 error_out << options.
h <<
" "
251 << error.array().abs().mean() <<
" "
252 << error.array().abs().maxCoeff() << std::endl;
259 trace.info() <<
"Mesh has " << viewmesh.nbVertex()
260 <<
" vertices and " << viewmesh.nbFaces() <<
" faces." << std::endl;
263 for(
unsigned int k = 0; k < viewmesh.nbFaces(); k++)
264 viewmesh.setFaceColor(k, colormap_error( error_faces(k) ));
282 Ball ball(
Point(0.0,0.0,0.0), 1.0);
288 return 2 * cos(p_s[1]) * cos(p_s[1]) * (2 * cos(p_s[2]) * cos(p_s[2]) - sin(p_s[2]) * sin(p_s[2]))
289 + 2 * (sin(p_s[1]) * sin(p_s[1]) - cos(p_s[1]) * cos(p_s[1]));
296 return - 2 * cos(p_s[2]);
302 return exp(p_sphere[0]);
309 if(p_s[1] == 0 && p_s[2] == 0)
return 1.;
311 const double theta_derivative = (sin(p_s[1]) * sin(p_s[1]) * sin(p_s[2])
312 - cos(p_s[1])) * (1 / sin(p_s[2])) * exp(p_sphere[0]);
313 const double phi_derivative = ( cos(p_s[1]) * cos(p_s[2]) * cos(p_s[2])
315 + cos(p_s[2]) * cos(p_s[2]) / sin(p_s[2])) * cos(p_s[1]) * exp(p_sphere[0]);
317 return theta_derivative + phi_derivative;
320 std::function<double(
const RealPoint3D&)> input_function
321 = ( options.
function == 0 ) ? xx_function : ( (options.
function == 1) ? cos_function : exp_function );
322 std::function<double(
const RealPoint3D&)> target_function
323 = ( options.
function == 0 ) ? xx_result : ( (options.
function == 1) ? cos_result : exp_result );
325 laplacian<Ball>(ball, options, input_function, target_function, argc, argv);
RealPoint getLowerBound() const
RealPoint getUpperBound() const
Aim: Model of the concept StarShaped3D represents any Sphere in the space.
Aim: This class template may be used to (linearly) convert scalar values in a given range into a colo...
Aim: Represents a set of n-1-cells in a nD space, together with adjacency relation between these cell...
Aim: A class for computing the Gauss digitization of some Euclidean shape, i.e. its intersection with...
std::set< SCell > SurfelSet
bool init(const Point &lower, const Point &upper, bool isClosed)
Specifies the upper and lower bounds for the maximal cells in this space.
static void digitalSurface2DualTriangulatedSurface(const DigitalSurface< DigitalSurfaceContainer > &dsurf, const CellEmbedder &cembedder, TriangulatedSurface< typename CellEmbedder::Value > &trisurf, VertexMap &vertexmap)
static void triangulatedSurface2Mesh(const TriangulatedSurface< Point > &trisurf, Mesh< Point > &mesh)
Aim: This class is defined to represent a surface mesh through a set of vertices and faces....
auto dot(const PointVector< dim, OtherComponent, OtherStorage > &v) const -> decltype(DGtal::dotProduct(*this, v))
Dot product with a PointVector.
double norm(const NormType type=L_2) const
void show() override
Starts the event loop and display of elements.
static void sMakeBoundary(SCellSet &aBoundary, const KSpace &aKSpace, const PointPredicate &pp, const Point &aLowerBound, const Point &aUpperBound)
Aim: Represent adjacencies between surfel elements, telling if it follows an interior to exterior ord...
HalfEdgeDataStructure::HalfEdgeIndex Arc
std::vector< Vertex > VertexRange
Vertex head(const Arc &a) const
FaceRange facesAroundVertex(const Vertex &v) const
Arc opposite(const Arc &a) const
VertexRange allVertices() const
Arc next(const Arc &a) const
ArcRange outArcs(const Vertex &v) const
VertexRange verticesAroundFace(const Face &f) const
Point & position(Vertex v)
std::vector< Arc > ArcRange
std::vector< Face > FaceRange
DigitalSurface< MyDigitalSurfaceContainer > MyDigitalSurface
MyDigitalSurface::SurfelSet SurfelSet
Z3i this namespace gathers the standard of types for 3D imagery.
HyperRectDomain< Space > Domain
KhalimskySpaceND< 3, Integer > KSpace
Space::RealVector RealVector
Space::RealPoint RealPoint
DGtal is the top-level namespace which contains all DGtal functions and types.
auto crossProduct(PointVector< 3, LeftEuclideanRing, LeftContainer > const &lhs, PointVector< 3, RightEuclideanRing, RightContainer > const &rhs) -> decltype(DGtal::constructFromArithmeticConversion(lhs, rhs))
Cross product of two 3D Points/Vectors.
void laplacian(Shape &shape, const Options &options, std::function< double(const RealPoint3D &)> input_function, std::function< double(const RealPoint3D &)> target_function, int argc, char **argv)
RealPoint cartesian_to_spherical(const RealPoint3D &a)
Z3i::RealPoint RealPoint3D
Z3i::RealVector RealVector3D
Aim: A trivial embedder for signed and unsigned cell, which corresponds to the canonic injection of c...
TriangulatedSurface< RealPoint > TriMesh