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/Viewer3D.h>
58 #include <DGtal/math/linalg/EigenSupport.h>
63 using namespace DGtal;
64 using namespace Eigen;
76 return RealPoint3D(a.norm(), atan2(a[1], a[0]), acos(a[2]));
85 std::string error_output;
88 template <
typename Shape>
90 std::function<
double(
const RealPoint3D&)> input_function,
91 std::function<
double(
const RealPoint3D&)> target_function,
92 int argc,
char** argv)
99 digitizer.attach(shape);
105 trace.
beginBlock(
"Construct the Khalimsky space from the image domain." );
111 trace.
error() <<
"Error in the Khamisky space construction."<<std::endl;
119 MySurfelAdjacency surfAdj(
true );
124 MySetOfSurfels theSetOfSurfels( kspace, surfAdj );
130 trace.
info() <<
"Digital surface has " << digSurf.
size() <<
" surfels."
137 typedef std::map< MyDigitalSurface::Vertex, TriMesh::Index > VertexMap;
142 MeshHelpers::digitalSurface2DualTriangulatedSurface
143 ( digSurf, canonicCellembedder, trimesh, vmap );
145 trace.
info() <<
"Triangulated surface is " << trimesh << std::endl;
153 for(
auto v : vertices)
156 if(options.smooth == 1)
160 std::ofstream error_out(options.error_output, std::ofstream::app);
161 std::ofstream function_out(
"function.dat");
163 Eigen::VectorXd laplacian_result( trimesh.
nbVertices() );
164 Eigen::VectorXd error( trimesh.
nbVertices() );
165 Eigen::VectorXd error_faces( trimesh.
nbFaces() );
167 double total_area = 0.;
170 for(
auto v : vertices)
177 for(
auto a : out_arcs)
201 const double alpha = acos( v1.dot(v2) );
202 const double beta = acos( vv1.
dot(vv2) );
205 accum += (tan(M_PI_2 - alpha) + tan(M_PI_2 - beta)) * (input_function(p_j) - input_function(p_i));
208 double accum_area = 0.;
210 for(
auto f : faces_around)
218 const double faceArea = .5 * cross.
norm();
220 accum_area += faceArea / 3.;
223 total_area += accum_area;
225 (options.smooth == 1)
226 ? laplacian_result(i) = (1 / (2. * accum_area)) * accum
227 : laplacian_result(i) = .5 * accum;
230 const double real_laplacian_value = target_function(w_projected);
234 function_out << w_s[1] <<
" "
236 << laplacian_result(i) <<
" "
237 << real_laplacian_value <<
" "
238 << input_function(p_i) << std::endl;
240 error(i) = laplacian_result(i) - real_laplacian_value;
241 for(
auto f : faces_around)
242 error_faces( f ) += error(i) / faces_around.size();
249 trace.
info() <<
"Computed Area / Real Area : " << total_area <<
" " << 4 * M_PI << std::endl;
250 trace.
info() <<
"Mean Error / Max Error : "
251 << error.array().abs().mean() <<
" / " << error.array().abs().maxCoeff(&max_index) << std::endl;
252 error_out << options.h <<
" "
253 << error.array().abs().mean() <<
" "
254 << error.array().abs().maxCoeff() << std::endl;
260 MeshHelpers::triangulatedSurface2Mesh( trimesh, viewmesh );
261 trace.
info() <<
"Mesh has " << viewmesh.nbVertex()
262 <<
" vertices and " << viewmesh.nbFaces() <<
" faces." << std::endl;
265 for(
unsigned int k = 0; k < viewmesh.nbFaces(); k++)
266 viewmesh.setFaceColor(k, colormap_error( error_faces(k) ));
268 QApplication application(argc,argv);
272 viewer << Viewer3D<>::updateDisplay;
278 int main(
int argc,
char **argv)
282 options.function = 0;
284 options.error_output =
"error_out.dat";
293 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]))
294 + 2 * (sin(p_s[1]) * sin(p_s[1]) - cos(p_s[1]) * cos(p_s[1]));
301 return - 2 * cos(p_s[2]);
307 return exp(p_sphere[0]);
314 if(p_s[1] == 0 && p_s[2] == 0)
return 1.;
316 const double theta_derivative = (sin(p_s[1]) * sin(p_s[1]) * sin(p_s[2])
317 - cos(p_s[1])) * (1 / sin(p_s[2])) * exp(p_sphere[0]);
318 const double phi_derivative = ( cos(p_s[1]) * cos(p_s[2]) * cos(p_s[2])
320 + cos(p_s[2]) * cos(p_s[2]) / sin(p_s[2])) * cos(p_s[1]) * exp(p_sphere[0]);
322 return theta_derivative + phi_derivative;
325 std::function<double(
const RealPoint3D&)> input_function
326 = ( options.function == 0 ) ? xx_function : ( (options.function == 1) ? cos_function : exp_function );
327 std::function<double(
const RealPoint3D&)> target_function
328 = ( options.function == 0 ) ? xx_result : ( (options.function == 1) ? cos_result : exp_result );
330 laplacian<Ball>(ball, options, input_function, target_function, argc, argv);
RealPoint getLowerBound() const
RealPoint getUpperBound() const
Aim: Model of the concept StarShaped represents any circle in the plane.
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...
const Point & lowerBound() const
const Point & upperBound() const
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...
std::set< SCell > SurfelSet
Preferred type for defining a set of surfels (always signed cells).
bool init(const Point &lower, const Point &upper, bool isClosed)
Specifies the upper and lower bounds for the maximal cells in this space.
Aim: This class is defined to represent a surface mesh through a set of vertices and faces....
Aim: Implements basic operations that will be used in Point and Vector classes.
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
Aim: A model of CDigitalSurfaceContainer which defines the digital surface as connected surfels....
Aim: A utility class for constructing surfaces (i.e. set of (n-1)-cells).
void beginBlock(const std::string &keyword="")
Aim: Represents a triangulated surface. The topology is stored with a half-edge data structure....
HalfEdgeDataStructure::HalfEdgeIndex Arc
std::vector< Vertex > VertexRange
Vertex head(const Arc &a) const
FaceRange facesAroundVertex(const Vertex &v) const
Point & position(Vertex v)
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
std::vector< Arc > ArcRange
std::vector< Face > FaceRange
virtual void show()
Overload QWidget method in order to add a call to updateList() method (to ensure that the lists are w...
DigitalSurface< MyDigitalSurfaceContainer > MyDigitalSurface
MyDigitalSurface::SurfelSet SurfelSet
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)
int main(int argc, char **argv)
Z3i::RealVector RealVector
Z3i::RealPoint RealPoint3D
Aim: A trivial embedder for signed and unsigned cell, which corresponds to the canonic injection of c...
TriangulatedSurface< RealPoint > TriMesh