CGAL users discussion list

["Laurent Rineau (CGAL/GeometryFactory)" <>]

On Monday, October 28, 2019 10:18:15 PM CET Gaetan wrote:
> 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-t
> d952799.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?

You are using the Surface_mesh package of CGAL. When using that package, the 
initialization is just the insert of a few points in the triangulation. You 
were right with that. However, just one point per component is probably not 
enough to initialize the Delaunay refinement process.

I have modified your program, and have added eight extra points (four per 
sphere), and it was enough. See the attached .cpp file.



--  
Laurent Rineau, PhD
R&D Engineer at GeometryFactory           http://www.geometryfactory.com/
Release Manager of the CGAL Project       http://www.cgal.org/
#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<Tr>;
using Polyhedron          = CGAL::Polyhedron_3<Kernel>;
using ImplicitSurface_3 = CGAL::Implicit_surface_3<Kernel,TwoSpheres<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);
  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);
  tr.insert({-4.,  2.,  0.});
  tr.insert({-4., -2.,  0.});
  tr.insert({-4.,  0.,  2.});
  tr.insert({-4.,  0., -2.});
  tr.insert({ 4.,  2.,  0.});
  tr.insert({ 4., -2.,  0.});
  tr.insert({ 4.,  0.,  2.});
  tr.insert({ 4.,  0., -2.});

  // 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<Tr> criteria(30.,  // angular bound
                                                                                       
0.1,  // radius bound
                                                                                       
0.1); // distance bound
  // Mesh surface
  CGAL::make_surface_mesh(c2t3, surface, criteria,
CGAL::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;

}