CGAL users discussion list

[Gaetan <>]

Hi all,

I want to be able to mesh 3d surfaces from an implicit function representing
several smooth closed surfaces.
From what I understand by reading the doc and a previous post
(http://cgal-discuss.949826.n4.nabble.com/About-the-surface-mesh-generator-td952799.html)
I need to add one initial sample point for each connected component to the
initial triangulation.

I did a test with two spheres (inspired from CGAL tests
@https://github.com/CGAL/cgal/blob/69fad29842b8a3a7316ea478b3a4da3a47d6bf6a/Surface_mesher/test/Surface_mesher/implicit_surface_mesher_test.cpp#L138),
but I can't get the final mesh with both the spheres in it.

<http://cgal-discuss.949826.n4.nabble.com/file/t376135/Capture.png> 

Is there something else I need to do after adding the initial points to the
triangulation or am I doing something wrong?

Thanks for your help,

What I did (see code below):
1 . create two implicit spheres
   - left  -> center = (-4,0,0), radius = 2
   - right -> center = (4, 0,0), radius = 2
2. Insert one seed point per sphere to the triangulation
  -   left   -> seed = (-2,0,0);
  -   right -> seed = (2,0,0);
3. Mesh the implicit surface representing both spheres
4. Output mesh to .off file


#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Surface_mesh_default_triangulation_3.h>
#include <CGAL/Complex_2_in_triangulation_3.h>
#include <CGAL/make_surface_mesh.h>
#include <CGAL/Implicit_surface_3.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/IO/facets_in_complex_2_to_triangle_mesh.h>

#include <string>
#include <iostream>
#include <fstream>

template <typename K>
class Sphere {
public:
  typedef typename K::FT FT;
  typedef typename K::Point_3 Point_3;
  typedef typename K::Sphere_3 Sphere_3;

  Sphere(Sphere_3 sphere)
    : sphere(sphere)
  {
  }

  FT operator()(const Point_3& p) const
  {
    FT x = p.x();
    FT y = p.y();
    FT z = p.z();
    x-=sphere.center().x();
    y-=sphere.center().y();
    z-=sphere.center().z();
    return x*x+y*y+z*z-sphere.squared_radius();
  }
private:
  const Sphere_3 sphere;
};

template <typename K>
class TwoSpheres 
{
public:
  typedef typename K::FT FT;
  typedef typename K::Point_3 Point_3;
  
  TwoSpheres(Sphere<K> sphere1, Sphere<K> sphere2)
    :  sphere1(sphere1), sphere2(sphere2)
  {
  }

  FT operator()(const Point_3& p) const
  {
    return sphere1(p)*sphere2(p);
  }
private:
  const Sphere<K> sphere1;
  const Sphere<K> sphere2;
};


// Aliases
using Kernel                 =
CGAL::Exact_predicates_inexact_constructions_kernel;
using Sphere_3            = Kernel::Sphere_3;
using Point_3               = Kernel::Point_3;
using FT                      = Kernel::FT;
using Tr                       = CGAL::Surface_mesh_default_triangulation_3;
using C2t3                   = CGAL::Complex_2_in_triangulation_3
;
using Polyhedron          = CGAL::Polyhedron_3<Kernel>;
using ImplicitSurface_3 = CGAL::Implicit_surface_3<Kernel,
TwoSpheres&lt;Kernel>>;


int main() {

  // Construct the spheres
  //  - left  -> center = (-4,0,0), radius = 2
  //  - right -> center = (4, 0,0), radius = 2
  const Sphere<Kernel> sphere_left (Sphere_3(Point_3(-4., 0., 0.), 4.));
  const Sphere<Kernel> sphere_right(Sphere_3(Point_3(4.,  0., 0.), 4.));

  // Assemble the two spheres
  const TwoSpheres<Kernel> two_spheres(sphere_left, sphere_right);

  // Create the seeds for the 3D-Delaunay triangulation
  //   - one seed per connected component
  const Point_3 point_left (-2,0,0);
  const Point_3 point_right(2,0,0);
  if(two_spheres(point_left) != 0)
    std::cerr << "ERROR: Left seed is not on sphere!" << std::endl;
  if(two_spheres(point_right) != 0)
    std::cerr << "ERROR: Right seed is not on sphere!" << std::endl;

  // Set the 3D-Delaunay triangulation
  Tr tr;
  tr.insert(point_left);
  tr.insert(point_right);

  // Set the 2D-complex in 3D-Delaunay triangulation
  C2t3 c2t3 (tr);   

  // Define the implicit surface
  const ImplicitSurface_3 surface(two_spheres,
                                                Sphere_3(CGAL::ORIGIN,
100));

  // Define meshing criteria
  const CGAL::Surface_mesh_default_criteria_3
 criteria(30.,  // angular bound
                                                                              
         
0.1,  // radius bound
                                                                              
         
0.1); // distance bound
  // Mesh surface
  CGAL::make_surface_mesh(c2t3, surface, criteria,
CGAL::Non_manifold_tag());
  std::cout << "Final number of points: " << tr.number_of_vertices() <<
std::endl;

  // Converts to polyhedron
  Polyhedron output_mesh;
  CGAL::facets_in_complex_2_to_triangle_mesh(c2t3, output_mesh);

  // Save to .off file
  const std::string output_filename = "test.off";
  std::cout << "Write file " << output_filename << std::endl;
  std::ofstream out(output_filename.c_str());
  out << output_mesh;

}






--
Sent from: http://cgal-discuss.949826.n4.nabble.com/