DGtal 1.3.0
Loading...
Searching...
No Matches
sphereCotangentLaplaceOperator.cpp
Go to the documentation of this file.
1
31
33
34#include <DGtal/helpers/StdDefs.h>
35
36#include <DGtal/base/BasicFunctors.h>
37
38#include <DGtal/topology/DigitalSurface.h>
39#include <DGtal/topology/SetOfSurfels.h>
40#include <DGtal/topology/CanonicCellEmbedder.h>
41#include <DGtal/topology/LightImplicitDigitalSurface.h>
42
43#include <DGtal/io/colormaps/ColorBrightnessColorMap.h>
44
45#include <DGtal/geometry/surfaces/estimation/LocalEstimatorFromSurfelFunctorAdapter.h>
46#include <DGtal/geometry/surfaces/estimation/IIGeometricFunctors.h>
47#include <DGtal/geometry/surfaces/estimation/IntegralInvariantCovarianceEstimator.h>
48
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"
55
56#include <DGtal/io/viewers/Viewer3D.h>
57
58#include <DGtal/math/linalg/EigenSupport.h>
59
61
62using namespace std;
63using namespace DGtal;
64using namespace Eigen;
65using namespace Z3i;
66
68
71
73{
74 return RealPoint3D(a.norm(), atan2(a[1], a[0]), acos(a[2]));
75}
76
77struct Options
78{
79 double h;
80 int function;
81 int smooth;
82
83 std::string error_output;
84};
85
86template <typename Shape>
87void laplacian(Shape& shape, const Options& options,
88 std::function<double(const RealPoint3D&)> input_function,
89 std::function<double(const RealPoint3D&)> target_function,
90 int argc, char** argv)
91{
92 trace.beginBlock("Laplacian");
93 typedef GaussDigitizer<Z3i::Space, Shape> Digitizer;
94
95 trace.beginBlock("Digitizing the shape");
96 Digitizer digitizer;
97 digitizer.attach(shape);
98 digitizer.init(shape.getLowerBound() + Z3i::Vector(-1,-1,-1), shape.getUpperBound() + Z3i::Vector(1,1,1), options.h);
99
100 Z3i::Domain domain = digitizer.getDomain();
101 trace.endBlock();
102
103 trace.beginBlock( "Construct the Khalimsky space from the image domain." );
104 Z3i::KSpace kspace;
105 bool space_ok = kspace.init( domain.lowerBound(),
106 domain.upperBound(), true );
107 if (!space_ok)
108 {
109 trace.error() << "Error in the Khamisky space construction."<<std::endl;
110 return;
111 }
112 trace.endBlock();
113
114 trace.beginBlock( "Extracting boundary by scanning the space. " );
115
116 typedef SurfelAdjacency<Z3i::KSpace::dimension> MySurfelAdjacency;
117 MySurfelAdjacency surfAdj( true ); // interior in all directions.
118
120 typedef SetOfSurfels< KSpace, SurfelSet > MySetOfSurfels;
122 MySetOfSurfels theSetOfSurfels( kspace, surfAdj );
123 Surfaces<Z3i::KSpace>::sMakeBoundary( theSetOfSurfels.surfelSet(),
124 kspace, digitizer,
126 domain.upperBound() );
127 MyDigitalSurface digSurf( theSetOfSurfels );
128 trace.info() << "Digital surface has " << digSurf.size() << " surfels."
129 << std::endl;
130 trace.endBlock();
131
132 trace.beginBlock( "Making triangulated surface. " );
135 typedef std::map< MyDigitalSurface::Vertex, TriMesh::Index > VertexMap;
136 TriMesh trimesh;
137 CanonicCellEmbedder canonicCellembedder(kspace);
138
139 VertexMap vmap; // stores the map Vertex -> Index
141 ( digSurf, canonicCellembedder, trimesh, vmap );
142
143 trace.info() << "Triangulated surface is " << trimesh << std::endl;
144
145 trace.endBlock();
146 trace.beginBlock("Computing Laplacian");
149
150 //Grid scaling factor and sphere projection
151 for(auto v : vertices)
152 {
153 trimesh.position( v ) *= options.h;
154 if(options.smooth == 1)
155 trimesh.position( v ) /= trimesh.position( v ).norm();
156 }
157
158 std::ofstream error_out(options.error_output, std::ofstream::app);
159 std::ofstream function_out("function.dat");
160
161 Eigen::VectorXd laplacian_result( trimesh.nbVertices() );
162 Eigen::VectorXd error( trimesh.nbVertices() );
163 Eigen::VectorXd error_faces( trimesh.nbFaces() );
164 int i = 0;
165 double total_area = 0.;
166
167 // Iteration over all vertices
168 for(auto v : vertices)
169 {
170 const RealPoint3D p_i = trimesh.position(v);
171 const TriangulatedSurface::ArcRange out_arcs = trimesh.outArcs(v);
172 double accum = 0.;
173
174 // We compute here \Delta f(p_i) by iteration over arcs going out from p_i
175 for(auto a : out_arcs)
176 {
177 // The point p_i -----> p_j
178 const RealPoint3D p_j = trimesh.position( trimesh.head(a) );
179
180 const TriangulatedSurface::Arc next_left_arc = trimesh.next( a );
181 const TriangulatedSurface::Arc next_right_arc = trimesh.next( trimesh.opposite(a) );
182
183 // Three points of the left triangle
184 const RealPoint3D p1 = p_j;
185 const RealPoint3D p2 = trimesh.position( trimesh.head( next_left_arc ) );
186 const RealPoint3D p3 = p_i;
187
188 // Three points of the right triangle
189 const RealPoint3D pp1 = p_j;
190 const RealPoint3D pp2 = trimesh.position(trimesh.head( next_right_arc ));
191 const RealPoint3D pp3 = p_i;
192
193 // Left and right angles
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();
198
199 const double alpha = acos( v1.dot(v2) );
200 const double beta = acos( vv1.dot(vv2) );
201
202 // Cotan accumulator
203 accum += (tan(M_PI_2 - alpha) + tan(M_PI_2 - beta)) * (input_function(p_j) - input_function(p_i));
204 }
205
206 double accum_area = 0.;
207 const TriangulatedSurface::FaceRange faces_around = trimesh.facesAroundVertex( v );
208 for(auto f : faces_around)
209 {
211 RealPoint3D p = trimesh.position(vr[0]);
212 RealPoint3D q = trimesh.position(vr[1]);
213 RealPoint3D r = trimesh.position(vr[2]);
214
215 const RealPoint3D cross = (r - p).crossProduct(r - q);
216 const double faceArea = .5 * cross.norm();
217
218 accum_area += faceArea / 3.;
219 }
220
221 total_area += accum_area;
222
223 (options.smooth == 1)
224 ? laplacian_result(i) = (1 / (2. * accum_area)) * accum
225 : laplacian_result(i) = .5 * accum;
226
227 const RealPoint3D w_projected = p_i / p_i.norm();
228 const double real_laplacian_value = target_function(w_projected);
229
230 const RealPoint3D w_s = cartesian_to_spherical(w_projected);
231
232 function_out << w_s[1] << " "
233 << w_s[2] << " "
234 << laplacian_result(i) << " "
235 << real_laplacian_value << " "
236 << input_function(p_i) << std::endl;
237
238 error(i) = laplacian_result(i) - real_laplacian_value;
239 for(auto f : faces_around)
240 error_faces( f ) += error(i) / faces_around.size();
241
242 i++;
243 }
244
245 int max_index;
246
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;
253
254 trace.endBlock();
255
256 typedef Mesh< CanonicCellEmbedder::Value > ViewMesh;
257 ViewMesh viewmesh;
258 MeshHelpers::triangulatedSurface2Mesh( trimesh, viewmesh );
259 trace.info() << "Mesh has " << viewmesh.nbVertex()
260 << " vertices and " << viewmesh.nbFaces() << " faces." << std::endl;
261
262 DGtal::ColorBrightnessColorMap<float> colormap_error(error_faces.minCoeff(), error_faces.maxCoeff(), DGtal::Color::Red);
263 for(unsigned int k = 0; k < viewmesh.nbFaces(); k++)
264 viewmesh.setFaceColor(k, colormap_error( error_faces(k) ));
265
266 QApplication application(argc,argv);
267 Viewer3D<> viewer;
268 viewer.show();
269 viewer << viewmesh;
270 viewer << Viewer3D<>::updateDisplay;
271 application.exec();
272
273 trace.endBlock();
274}
275
276int main(int argc, char **argv)
277{
278 Options options;
279 options.h = 0.1;
280 options.function = 0;
281 options.smooth = 0;
282 options.error_output = "error_out.dat";
283
284 typedef Ball3D<Z3i::Space> Ball;
285 Ball ball(Point(0.0,0.0,0.0), 1.0);
286
287 std::function<double(const RealPoint3D&)> xx_function = [](const RealPoint3D& p) {return p[0] * p[0];};
288 std::function<double(const RealPoint3D&)> xx_result = [](const RealPoint3D& p)
289 {
290 const RealPoint3D p_s = cartesian_to_spherical(p);
291 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]))
292 + 2 * (sin(p_s[1]) * sin(p_s[1]) - cos(p_s[1]) * cos(p_s[1]));
293 };
294
295 std::function<double(const RealPoint3D&)> cos_function = [](const RealPoint3D& p) {return p[2];};
296 std::function<double(const RealPoint3D&)> cos_result = [](const RealPoint3D& p)
297 {
298 const RealPoint3D p_s = cartesian_to_spherical(p);
299 return - 2 * cos(p_s[2]);
300 };
301
302 std::function<double(const RealPoint3D&)> exp_function = [](const RealPoint3D& p)
303 {
304 const RealPoint3D p_sphere = p / p.norm();
305 return exp(p_sphere[0]);
306 };
307 std::function<double(const RealPoint3D&)> exp_result = [](const RealPoint3D& p)
308 {
309 const RealPoint3D p_sphere = p / p.norm();
310 const RealPoint3D p_s = cartesian_to_spherical(p);
311
312 if(p_s[1] == 0 && p_s[2] == 0) return 1.;
313
314 const double theta_derivative = (sin(p_s[1]) * sin(p_s[1]) * sin(p_s[2])
315 - cos(p_s[1])) * (1 / sin(p_s[2])) * exp(p_sphere[0]);
316 const double phi_derivative = ( cos(p_s[1]) * cos(p_s[2]) * cos(p_s[2])
317 - sin(p_s[2])
318 + cos(p_s[2]) * cos(p_s[2]) / sin(p_s[2])) * cos(p_s[1]) * exp(p_sphere[0]);
319
320 return theta_derivative + phi_derivative;
321 };
322
323 std::function<double(const RealPoint3D&)> input_function
324 = ( options.function == 0 ) ? xx_function : ( (options.function == 1) ? cos_function : exp_function );
325 std::function<double(const RealPoint3D&)> target_function
326 = ( options.function == 0 ) ? xx_result : ( (options.function == 1) ? cos_result : exp_result );
327
328 laplacian<Ball>(ball, options, input_function, target_function, argc, argv);
329
330 return 0;
331}
332
RealPoint getLowerBound() const
Definition: Astroid2D.h:122
RealPoint getUpperBound() const
Definition: Astroid2D.h:131
Aim: Model of the concept StarShaped represents any circle in the plane.
Definition: Ball2D.h:61
Aim: Model of the concept StarShaped3D represents any Sphere in the space.
Definition: Ball3D.h:61
Aim: This class template may be used to (linearly) convert scalar values in a given range into a colo...
static const Color Red
Definition: Color.h:416
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.
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....
Definition: Mesh.h:92
Aim: Implements basic operations that will be used in Point and Vector classes.
Definition: PointVector.h:593
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....
Definition: SetOfSurfels.h:74
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...
void beginBlock(const std::string &keyword="")
std::ostream & error()
std::ostream & info()
double endBlock()
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
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)
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
Space::Point Point
Definition: StdDefs.h:168
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.
Trace trace
Definition: Common.h:154
STL namespace.
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...
int main()
Definition: testBits.cpp:56
Domain domain
Ball2D< Space > Ball
TriangulatedSurface< RealPoint > TriMesh