DGtal  1.2.0
Typedefs | Functions
sphereCotangentLaplaceOperator.cpp File Reference
#include <DGtal/helpers/StdDefs.h>
#include <DGtal/base/BasicFunctors.h>
#include <DGtal/topology/DigitalSurface.h>
#include <DGtal/topology/SetOfSurfels.h>
#include <DGtal/topology/CanonicCellEmbedder.h>
#include <DGtal/topology/LightImplicitDigitalSurface.h>
#include <DGtal/io/colormaps/ColorBrightnessColorMap.h>
#include <DGtal/geometry/surfaces/estimation/LocalEstimatorFromSurfelFunctorAdapter.h>
#include <DGtal/geometry/surfaces/estimation/IIGeometricFunctors.h>
#include <DGtal/geometry/surfaces/estimation/IntegralInvariantCovarianceEstimator.h>
#include <DGtal/shapes/parametric/Ball3D.h>
#include <DGtal/shapes/Shapes.h>
#include <DGtal/shapes/GaussDigitizer.h>
#include <DGtal/shapes/Mesh.h>
#include <DGtal/shapes/MeshHelpers.h>
#include "DGtal/shapes/TriangulatedSurface.h"
#include <DGtal/io/viewers/Viewer3D.h>
#include <DGtal/math/linalg/EigenSupport.h>
Include dependency graph for sphereCotangentLaplaceOperator.cpp:

Go to the source code of this file.

Typedefs

typedef Z3i::Space Space
 
typedef Z3i::KSpace KSpace
 
typedef Z3i::RealVector RealVector
 
typedef Z3i::RealPoint RealPoint3D
 

Functions

RealPoint cartesian_to_spherical (const RealPoint3D &a)
 
template<typename Shape >
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)
 
int main (int argc, char **argv)
 

Detailed Description

This program is free software : you can redistribuate 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 you option) any later version.

This program is distributed int 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, set http://www.gnu.org/license/.

Author
Thomas Caissard (thoma.nosp@m.s.ca.nosp@m.issar.nosp@m.d@li.nosp@m.ris.c.nosp@m.nrs..nosp@m.fr ) Laboratoire d'InfoRmatique en Image et Systèmes d'information - LIRIS (CNRS, UMR 5205), INSA-Lyon, France LAboratoire de MAthématiques - LAMA (CNRS, UMR 5127), Université de Savoie, France
Date
2017/06/20

An example file named exampleTriangulatedSurface.

This file is part of the DGtal library.

Definition in file sphereCotangentLaplaceOperator.cpp.

Typedef Documentation

◆ KSpace

Definition at line 70 of file sphereCotangentLaplaceOperator.cpp.

◆ RealPoint3D

Definition at line 72 of file sphereCotangentLaplaceOperator.cpp.

◆ RealVector

◆ Space

typedef Z3i::Space Space

Definition at line 69 of file sphereCotangentLaplaceOperator.cpp.

Function Documentation

◆ cartesian_to_spherical()

RealPoint cartesian_to_spherical ( const RealPoint3D a)
Examples
dec/exampleHeatLaplace.cpp.

Definition at line 74 of file sphereCotangentLaplaceOperator.cpp.

75 {
76  return RealPoint3D(a.norm(), atan2(a[1], a[0]), acos(a[2]));
77 }
Z3i::RealPoint RealPoint3D

Referenced by laplacian(), and main().

◆ laplacian()

template<typename Shape >
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 
)

Definition at line 89 of file sphereCotangentLaplaceOperator.cpp.

93 {
94  trace.beginBlock("Laplacian");
95  typedef GaussDigitizer<Z3i::Space, Shape> Digitizer;
96 
97  trace.beginBlock("Digitizing the shape");
98  Digitizer digitizer;
99  digitizer.attach(shape);
100  digitizer.init(shape.getLowerBound() + Z3i::Vector(-1,-1,-1), shape.getUpperBound() + Z3i::Vector(1,1,1), options.h);
101 
102  Z3i::Domain domain = digitizer.getDomain();
103  trace.endBlock();
104 
105  trace.beginBlock( "Construct the Khalimsky space from the image domain." );
106  KSpace kspace;
107  bool space_ok = kspace.init( domain.lowerBound(),
108  domain.upperBound(), true );
109  if (!space_ok)
110  {
111  trace.error() << "Error in the Khamisky space construction."<<std::endl;
112  return;
113  }
114  trace.endBlock();
115 
116  trace.beginBlock( "Extracting boundary by scanning the space. " );
117 
118  typedef SurfelAdjacency<KSpace::dimension> MySurfelAdjacency;
119  MySurfelAdjacency surfAdj( true ); // interior in all directions.
120 
122  typedef SetOfSurfels< KSpace, SurfelSet > MySetOfSurfels;
124  MySetOfSurfels theSetOfSurfels( kspace, surfAdj );
125  Surfaces<KSpace>::sMakeBoundary( theSetOfSurfels.surfelSet(),
126  kspace, digitizer,
127  domain.lowerBound(),
128  domain.upperBound() );
129  MyDigitalSurface digSurf( theSetOfSurfels );
130  trace.info() << "Digital surface has " << digSurf.size() << " surfels."
131  << std::endl;
132  trace.endBlock();
133 
134  trace.beginBlock( "Making triangulated surface. " );
137  typedef std::map< MyDigitalSurface::Vertex, TriMesh::Index > VertexMap;
138  TriMesh trimesh;
139  CanonicCellEmbedder canonicCellembedder(kspace);
140 
141  VertexMap vmap; // stores the map Vertex -> Index
142  MeshHelpers::digitalSurface2DualTriangulatedSurface
143  ( digSurf, canonicCellembedder, trimesh, vmap );
144 
145  trace.info() << "Triangulated surface is " << trimesh << std::endl;
146 
147  trace.endBlock();
148  trace.beginBlock("Computing Laplacian");
151 
152  //Grid scaling factor and sphere projection
153  for(auto v : vertices)
154  {
155  trimesh.position( v ) *= options.h;
156  if(options.smooth == 1)
157  trimesh.position( v ) /= trimesh.position( v ).norm();
158  }
159 
160  std::ofstream error_out(options.error_output, std::ofstream::app);
161  std::ofstream function_out("function.dat");
162 
163  Eigen::VectorXd laplacian_result( trimesh.nbVertices() );
164  Eigen::VectorXd error( trimesh.nbVertices() );
165  Eigen::VectorXd error_faces( trimesh.nbFaces() );
166  int i = 0;
167  double total_area = 0.;
168 
169  // Iteration over all vertices
170  for(auto v : vertices)
171  {
172  const RealPoint3D p_i = trimesh.position(v);
173  const TriangulatedSurface::ArcRange out_arcs = trimesh.outArcs(v);
174  double accum = 0.;
175 
176  // We compute here \Delta f(p_i) by iteration over arcs going out from p_i
177  for(auto a : out_arcs)
178  {
179  // The point p_i -----> p_j
180  const RealPoint3D p_j = trimesh.position( trimesh.head(a) );
181 
182  const TriangulatedSurface::Arc next_left_arc = trimesh.next( a );
183  const TriangulatedSurface::Arc next_right_arc = trimesh.next( trimesh.opposite(a) );
184 
185  // Three points of the left triangle
186  const RealPoint3D p1 = p_j;
187  const RealPoint3D p2 = trimesh.position( trimesh.head( next_left_arc ) );
188  const RealPoint3D p3 = p_i;
189 
190  // Three points of the right triangle
191  const RealPoint3D pp1 = p_j;
192  const RealPoint3D pp2 = trimesh.position(trimesh.head( next_right_arc ));
193  const RealPoint3D pp3 = p_i;
194 
195  // Left and right angles
196  const RealPoint3D v1 = (p1 - p2) / (p1 - p2).norm();
197  const RealPoint3D v2 = (p3 - p2) / (p3 - p2).norm();
198  const RealPoint3D vv1 = (pp1 - pp2) / (pp1 - pp2).norm();
199  const RealPoint3D vv2 = (pp3 - pp2) / (pp3 - pp2).norm();
200 
201  const double alpha = acos( v1.dot(v2) );
202  const double beta = acos( vv1.dot(vv2) );
203 
204  // Cotan accumulator
205  accum += (tan(M_PI_2 - alpha) + tan(M_PI_2 - beta)) * (input_function(p_j) - input_function(p_i));
206  }
207 
208  double accum_area = 0.;
209  const TriangulatedSurface::FaceRange faces_around = trimesh.facesAroundVertex( v );
210  for(auto f : faces_around)
211  {
213  RealPoint3D p = trimesh.position(vr[0]);
214  RealPoint3D q = trimesh.position(vr[1]);
215  RealPoint3D r = trimesh.position(vr[2]);
216 
217  const RealPoint3D cross = (r - p).crossProduct(r - q);
218  const double faceArea = .5 * cross.norm();
219 
220  accum_area += faceArea / 3.;
221  }
222 
223  total_area += accum_area;
224 
225  (options.smooth == 1)
226  ? laplacian_result(i) = (1 / (2. * accum_area)) * accum
227  : laplacian_result(i) = .5 * accum;
228 
229  const RealPoint3D w_projected = p_i / p_i.norm();
230  const double real_laplacian_value = target_function(w_projected);
231 
232  const RealPoint3D w_s = cartesian_to_spherical(w_projected);
233 
234  function_out << w_s[1] << " "
235  << w_s[2] << " "
236  << laplacian_result(i) << " "
237  << real_laplacian_value << " "
238  << input_function(p_i) << std::endl;
239 
240  error(i) = laplacian_result(i) - real_laplacian_value;
241  for(auto f : faces_around)
242  error_faces( f ) += error(i) / faces_around.size();
243 
244  i++;
245  }
246 
247  int max_index;
248 
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;
255 
256  trace.endBlock();
257 
258  typedef Mesh< CanonicCellEmbedder::Value > ViewMesh;
259  ViewMesh viewmesh;
260  MeshHelpers::triangulatedSurface2Mesh( trimesh, viewmesh );
261  trace.info() << "Mesh has " << viewmesh.nbVertex()
262  << " vertices and " << viewmesh.nbFaces() << " faces." << std::endl;
263 
264  DGtal::ColorBrightnessColorMap<float> colormap_error(error_faces.minCoeff(), error_faces.maxCoeff(), DGtal::Color::Red);
265  for(unsigned int k = 0; k < viewmesh.nbFaces(); k++)
266  viewmesh.setFaceColor(k, colormap_error( error_faces(k) ));
267 
268  QApplication application(argc,argv);
269  Viewer3D<> viewer;
270  viewer.show();
271  viewer << viewmesh;
272  viewer << Viewer3D<>::updateDisplay;
273  application.exec();
274 
275  trace.endBlock();
276 }
RealPoint getLowerBound() const
Definition: Astroid2D.h:122
RealPoint getUpperBound() const
Definition: Astroid2D.h:131
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:392
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....
Definition: Mesh.h:92
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
Aim: A utility class for constructing surfaces (i.e. set of (n-1)-cells).
Definition: Surfaces.h:79
std::ostream & error()
void beginBlock(const std::string &keyword="")
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
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
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
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
MessageStream error
std::pair< typename graph_traits< DGtal::DigitalSurface< TDigitalSurfaceContainer > >::vertex_iterator, typename graph_traits< DGtal::DigitalSurface< TDigitalSurfaceContainer > >::vertex_iterator > vertices(const DGtal::DigitalSurface< TDigitalSurfaceContainer > &digSurf)
RealPoint cartesian_to_spherical(const RealPoint3D &a)
Aim: A trivial embedder for signed and unsigned cell, which corresponds to the canonic injection of c...
Domain domain
TriangulatedSurface< RealPoint > TriMesh

References DGtal::TriangulatedSurface< TPoint >::allVertices(), DGtal::Trace::beginBlock(), cartesian_to_spherical(), DGtal::crossProduct(), domain, DGtal::PointVector< dim, TEuclideanRing, TContainer >::dot(), DGtal::Trace::endBlock(), DGtal::Trace::error(), DGtal::TriangulatedSurface< TPoint >::facesAroundVertex(), DGtal::Astroid2D< TSpace >::getLowerBound(), DGtal::Astroid2D< TSpace >::getUpperBound(), DGtal::TriangulatedSurface< TPoint >::head(), DGtal::Trace::info(), DGtal::KhalimskySpaceND< dim, TInteger >::init(), DGtal::HyperRectDomain< TSpace >::lowerBound(), DGtal::TriangulatedSurface< TPoint >::nbFaces(), DGtal::TriangulatedSurface< TPoint >::nbVertices(), DGtal::TriangulatedSurface< TPoint >::next(), DGtal::PointVector< dim, TEuclideanRing, TContainer >::norm(), DGtal::TriangulatedSurface< TPoint >::opposite(), DGtal::TriangulatedSurface< TPoint >::outArcs(), DGtal::TriangulatedSurface< TPoint >::position(), DGtal::Color::Red, DGtal::Viewer3D< TSpace, TKSpace >::show(), DGtal::DigitalSurface< TDigitalSurfaceContainer >::size(), DGtal::trace, DGtal::HyperRectDomain< TSpace >::upperBound(), and DGtal::TriangulatedSurface< TPoint >::verticesAroundFace().

◆ main()

int main ( int  argc,
char **  argv 
)

Definition at line 278 of file sphereCotangentLaplaceOperator.cpp.

279 {
280  Options options;
281  options.h = 0.1;
282  options.function = 0;
283  options.smooth = 0;
284  options.error_output = "error_out.dat";
285 
286  typedef Ball3D<Z3i::Space> Ball;
287  Ball ball(Point(0.0,0.0,0.0), 1.0);
288 
289  std::function<double(const RealPoint3D&)> xx_function = [](const RealPoint3D& p) {return p[0] * p[0];};
290  std::function<double(const RealPoint3D&)> xx_result = [](const RealPoint3D& p)
291  {
292  const RealPoint3D p_s = cartesian_to_spherical(p);
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]));
295  };
296 
297  std::function<double(const RealPoint3D&)> cos_function = [](const RealPoint3D& p) {return p[2];};
298  std::function<double(const RealPoint3D&)> cos_result = [](const RealPoint3D& p)
299  {
300  const RealPoint3D p_s = cartesian_to_spherical(p);
301  return - 2 * cos(p_s[2]);
302  };
303 
304  std::function<double(const RealPoint3D&)> exp_function = [](const RealPoint3D& p)
305  {
306  const RealPoint3D p_sphere = p / p.norm();
307  return exp(p_sphere[0]);
308  };
309  std::function<double(const RealPoint3D&)> exp_result = [](const RealPoint3D& p)
310  {
311  const RealPoint3D p_sphere = p / p.norm();
312  const RealPoint3D p_s = cartesian_to_spherical(p);
313 
314  if(p_s[1] == 0 && p_s[2] == 0) return 1.;
315 
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])
319  - sin(p_s[2])
320  + cos(p_s[2]) * cos(p_s[2]) / sin(p_s[2])) * cos(p_s[1]) * exp(p_sphere[0]);
321 
322  return theta_derivative + phi_derivative;
323  };
324 
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 );
329 
330  laplacian<Ball>(ball, options, input_function, target_function, argc, argv);
331 
332  return 0;
333 }
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
MyPointD Point
Definition: testClone2.cpp:383
Ball2D< Space > Ball

References cartesian_to_spherical(), and DGtal::PointVector< dim, TEuclideanRing, TContainer >::norm().