3D surface mesh generation of two implicit spheres

classic Classic list List threaded Threaded
5 messages Options
Reply | Threaded
Open this post in threaded view
|

3D surface mesh generation of two implicit spheres

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/

--
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://sympa.inria.fr/sympa/info/cgal-discuss


Reply | Threaded
Open this post in threaded view
|

Re: 3D surface mesh generation of two implicit spheres

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/

--
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://sympa.inria.fr/sympa/info/cgal-discuss



main.cpp (4K) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: 3D surface mesh generation of two implicit spheres

Gaetan
Hi Laurent,

Thanks for your prompt answer, I see your example works fine.
Therefore my next question: how can we know the minimum number of initial
points required per component?
Is it more related to the initial cells (one initial cell per component) or
is it something else?




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

--
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://sympa.inria.fr/sympa/info/cgal-discuss


Reply | Threaded
Open this post in threaded view
|

Re: 3D surface mesh generation of two implicit spheres

Laurent Rineau (CGAL/GeometryFactory)
On Tuesday, October 29, 2019 3:48:01 PM CET Gaetan wrote:
> Hi Laurent,
>
> Thanks for your prompt answer, I see your example works fine.
> Therefore my next question: how can we know the minimum number of initial
> points required per component?
> Is it more related to the initial cells (one initial cell per component) or
> is it something else?

Actually, there is no minimal number of points, in general.

In theory, three points per component is sufficient, if the size of the
triangle is sufficiently small compared to the they are close enough. See that
article for the details:

  Boissonnat, Jean-Daniel, and Steve Oudot. "Provably good sampling and
meshing of surfaces." Graphical Models 67.5 (2005): 405-451.

https://geometrica.saclay.inria.fr/team/Steve.Oudot/papers/bo-pgsms-05/bo-pgsms-05.pdf


Imagine you are meshing a torus, and you insert N points on one circle of the
torus, then the triangulation of those N points will be degenerated, and for
all triangles of that triangulation, its dual will not intersect the torus at
all. And that whatever the initial number of points.

Practically, the Surface_mesher package of CGAL tries to find 20 random points
on the surface, and that was enough in all cases. SO I suggest you choose 20
random points per component.


If you tell me more about you real need, I might find a better solution that
random points, for the initialization.

--
Laurent Rineau, PhD
R&D Engineer at GeometryFactory           http://www.geometryfactory.com/
Release Manager of the CGAL Project       http://www.cgal.org/




--
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://sympa.inria.fr/sympa/info/cgal-discuss


Reply | Threaded
Open this post in threaded view
|

Re: 3D surface mesh generation of two implicit spheres

Gaetan
What you describe makes perfect sense.
Our objective is to perform a Poisson surface reconstruction of several
connected components (see my old post @
http://cgal-discuss.949826.n4.nabble.com/Poisson-reconstruction-of-multiple-connected-components-td4664371.html#a4664383).
I did a couple of tests using 20 initial points per component and it seems
to work fine.
Indeed, it is very unlikely that we will end up with 20 degenerate points
within one component.

Thanks again for your help,
G.



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

--
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://sympa.inria.fr/sympa/info/cgal-discuss