CGAL users discussion list

[Maarten Moesen <>]

Dear CGAL,

To try the 3D mesher in CGAL 4.1, I
was generating a mesh for a sphere that intersects a box using a labelled
domain. Different labels are given to the sphere, the box, and the exterior
of the box. The exterior of the box is to be discarded, the inside mesh
is what matters to me.

I noticed that there are almost always
some tetrahedra missing on the sides of the box. Sometimes the tetrahedra
are wrongly labeled, i.e. tetrahedra that are clearly outside the sphere
are labeled as inside and vice-versa. I personally think that the wrong
labeling/selection is at the heart of the problem. Anyway, the resulting
holes and pits are quite large, too large in fact :(.

My question:
* How can I solve this problem most
easily?

In the attachment you can find a screenshot
and below you can find the source code.

Specs: CGAL 4.1 on a 64-bit Windows
7 machine with a 64-bit MinGW GCC 4.7.1 compiler (x86_64-w64-mingw32).
With this configuration, the Mesh_3 examples work fine. Compiled in release
mode.

Best regards,

Maarten


#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>

#include <CGAL/Mesh_triangulation_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>

#include <CGAL/Mesh_domain_with_polyline_features_3.h>
#include <CGAL/Implicit_mesh_domain_3.h>
#include <CGAL/Mesh_3/Labeled_mesh_domain_3.h>

#include <CGAL/make_mesh_3.h>
#include <CGAL/Timer.h>

namespace BUBS
{
const double PI = 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798;
}

// Domain
typedef CGAL::Exact_predicates_inexact_constructions_kernel
K;
typedef K::FT FT;
typedef K::Point_3 Point;

template<class BGT,
         typename
Return_type = int >
class Cubic_labeled_function
{
  public:
    // Types
    typedef Return_type return_type;
    typedef typename BGT::FT
       FT;
    typedef typename BGT::Point_3
  Point_3;
    typedef typename BGT::Sphere_3
 Sphere_3;

  /// Constructor
    Cubic_labeled_function(FT
const& radius = 0.6): m_rad(radius), m_rad_squared(radius*radius) {}

  // Default copy constructor and
assignment operator are ok

  /// Destructor
  ~Cubic_labeled_function() {}

  /**
   * Returns an int corresponding
to the label at point \c p
   * @param p the input point
   * @return the label at
point \c p
   */
  return_type operator()(const
Point_3& p, const bool = true) const
  {
    if (p.x() < 0 || p.x()
> 1 || p.y() < 0 || p.y() > 1 || p.z() < 0 || p.z() > 1)
      return 0;

    const Point_3 center(0.5,
0.5, 0.5); 
    if (squared_distance(p,
center) < squared_radius())
      return 2;

    return 1;
  }

  FT x_min() const { return 0;
}
  FT x_max() const { return 1;
}
  FT y_min() const { return 0;
}
  FT y_max() const { return 1;
}
  FT z_min() const { return 0;
}
  FT z_max() const { return 1;
}

   Sphere_3 bounding_sphere()
const
  {
    return Sphere_3(Point_3(0.04,
0.04, 0.04), 3.0); 
  }

  FT const& radius() const
{ return m_rad; }
  FT const& squared_radius()
const { return m_rad_squared; }

private:
  FT m_rad, m_rad_squared;
};

class Cubic_function : public std::unary_function<Point,FT>

{
public:
  Cubic_function(double const&
radius = 0.6): m_x0(0), m_y0(0), m_z0(0), m_x1(1), m_y1(1), m_z1(1), m_xc(0.5),
m_yc(0.5), m_zc(0.5), m_r(radius) {}

  FT const& x_min() const {
return m_x0; }
  FT const& x_max() const {
return m_x1; }
  FT const& y_min() const {
return m_y0; }
  FT const& y_max() const {
return m_y1; }
  FT const& z_min() const {
return m_z0; }
  FT const& z_max() const {
return m_z1; }
  FT const& x_cen() const {
return m_xc; }
  FT const& y_cen() const {
return m_yc; }
  FT const& z_cen() const {
return m_zc; }

  FT const& radius() const
{ return m_r; }
  
  FT square(FT const& x) const
  {
    return x*x;
  }

  FT operator()(const Point&
p) const
  {
    return std::max(std::max(std::max(std::max(x_min()-p.x(),
p.x()-x_max()), std::max(y_min()-p.y(), p.y()-y_max())), std::max(z_min()-p.z(),
p.z()-z_max())),
           
        square(radius()) - square(x_cen()-p.x()) -
square(y_cen()-p.y()) - square(z_cen()-p.z()) ); 
  }

  K::Sphere_3 bounding_sphere()
const
  {
    return K::Sphere_3(Point(0.04,
0.04, 0.04), 3.0);
    
  }

private:
  FT m_x0, m_x1, m_y0, m_y1, m_z0,
m_z1, m_xc, m_yc, m_zc, m_r;
};

//typedef FT (Function)(const Point&);
typedef Cubic_labeled_function<K>
CLF;
typedef CGAL::Implicit_mesh_domain_3<Cubic_function,K>
Implicit_mesh_domain;
typedef CGAL::Mesh_3::Labeled_mesh_domain_3<CLF,K>
Labeled_mesh_domain;
//typedef CGAL::Mesh_domain_with_polyline_features_3<Implicit_mesh_domain>
Mesh_domain;
typedef CGAL::Mesh_domain_with_polyline_features_3<Labeled_mesh_domain>
Mesh_domain;

// Triangulation
typedef CGAL::Mesh_triangulation_3<Mesh_domain>::type
Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr,Mesh_domain::Corner_index,Mesh_domain::Curve_segment_index>
C3t3;

// Criteria
typedef CGAL::Mesh_criteria_3<Tr>
Mesh_criteria;

typedef std::vector<Point>  
     Polyline;
typedef std::vector<Polyline>
    Polylines;

int main()
{
  CGAL::Timer timer;

  CLF fun(0.6);
  Mesh_domain domain(fun, fun.bounding_sphere(),
1e-6);

  double edge_length = 0.05;
  // Mesh criteria
  Mesh_criteria criteria(CGAL::parameters::facet_angle=30,
    CGAL::parameters::edge_size=edge_length,
    CGAL::parameters::facet_size=edge_length*2,

    CGAL::parameters::facet_distance=0.01,
    CGAL::parameters::cell_radius_edge_ratio=2,
    CGAL::parameters::cell_size=edge_length*2);
  

  // The ribs of the cube
  Point p000(fun.x_min(), fun.y_min(),
fun.z_min());
  Point p001(fun.x_min(), fun.y_min(),
fun.z_max());  
  Point p010(fun.x_min(), fun.y_max(),
fun.z_min());
  Point p011(fun.x_min(), fun.y_max(),
fun.z_max());
  Point p100(fun.x_max(), fun.y_min(),
fun.z_min());
  Point p101(fun.x_max(), fun.y_min(),
fun.z_max());  
  Point p110(fun.x_max(), fun.y_max(),
fun.z_min());
  Point p111(fun.x_max(), fun.y_max(),
fun.z_max());
  
  // The edges of the cube
  Polylines polylines(18);
  polylines[0].push_back(p000);
polylines[0].push_back(p100);
  polylines[1].push_back(p001);
polylines[1].push_back(p101);
  polylines[2].push_back(p010);
polylines[2].push_back(p110);
  polylines[3].push_back(p011);
polylines[3].push_back(p111);
  polylines[4].push_back(p000);
polylines[4].push_back(p010);
  polylines[5].push_back(p001);
polylines[5].push_back(p011);
  polylines[6].push_back(p100);
polylines[6].push_back(p110);
  polylines[7].push_back(p101);
polylines[7].push_back(p111);
  polylines[8].push_back(p000);
polylines[8].push_back(p001);
  polylines[9].push_back(p010);
polylines[9].push_back(p011);
  polylines[10].push_back(p100);
polylines[10].push_back(p101);
  polylines[11].push_back(p110);
polylines[11].push_back(p111);
 
  // The intersections of the sphere
with the cube
  double rad = CGAL::to_double(fun.radius());
  rad = std::sqrt(rad*rad - 0.25);

  int n = std::max(8, static_cast<int>(std::ceil(2*BUBS::PI*rad/edge_length)));

  double* cs = new double[n+1];
  double* ss = new double[n+1];
  for (int ci = 0; ci < n; ++ci)
  {
    double angle = 2*BUBS::PI*ci/n;
    cs[ci] = 0.5 + rad*std::cos(angle);
    ss[ci] = 0.5 + rad*std::sin(angle);
  }
  cs[n] = cs[0];
  ss[n] = ss[0];

  std::cout << "Radius:
" << rad << "\nNumber of points: " <<
n << std::endl;

  for (int ci = 0; ci <= n;
++ci)
  {
    //std::cout << fun(Point(0.0,
cs[ci], ss[ci])) << "   " << fun(Point(1.0,
cs[ci], ss[ci])) << std::endl;
    polylines[12].push_back(Point(0.0,
cs[ci], ss[ci]));
    polylines[13].push_back(Point(1.0,
cs[ci], ss[ci]));
    polylines[14].push_back(Point(cs[ci],
0.0, ss[ci]));
    polylines[15].push_back(Point(cs[ci],
1.0, ss[ci]));
    polylines[16].push_back(Point(cs[ci],
ss[ci], 0.0));
    polylines[17].push_back(Point(cs[ci],
ss[ci], 1.0));
  }
  
  delete [] cs;
  delete [] ss;

  // Add the constraint segments
to the domain
  domain.add_features(polylines.begin(),
polylines.end());
  
  // Mesh generation
  std::cout << "Generating
the mesh: ";
  timer.reset();timer.start();
  C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain,
criteria, CGAL::parameters::no_exude(), CGAL::parameters::no_perturb());
  timer.stop(); 
  std::cout << timer.time()
<< " s." << std::endl;

  std::cout << "ODT
Smoothing: ";
  timer.reset();timer.start();
  CGAL::Mesh_optimization_return_code
code = CGAL::odt_optimize_mesh_3(c3t3, domain);
  timer.stop(); 
  std::cout << timer.time()
<< " s. Return code: " << code << std::endl;

  std::cout << "Perturbation:
";
  timer.reset();timer.start();
  code = CGAL::perturb_mesh_3(c3t3,
domain);
  timer.stop(); 
  std::cout << timer.time()
<< " s. Return code: " << code << std::endl;

  std::cout << "Sliver
exudation: ";
  timer.reset();timer.start();
  code = CGAL::exude_mesh_3(c3t3);
  timer.stop(); 
  std::cout << timer.time()
<< " s. Return code: " << code << std::endl;
  
  // Output
  std::ofstream medit_file("sphere_in_box.mesh");
  c3t3.output_to_medit(medit_file);
  medit_file.close();
  std::cout << "Mesh
saved to sphere_in_box.mesh" << std::endl;

  std::cout << "Number
of cells in complex: " << c3t3.number_of_cells() << std::endl;
  std::cout << "Number
of cells in triangulation: " << c3t3.triangulation().number_of_cells()
<< std::endl;
  std::cout << "Number
of finite cells in triangulation: " << c3t3.triangulation().number_of_finite_cells()
<< std::endl;
  std::cout << "Number
of vertices in triangulation: " << c3t3.triangulation().number_of_vertices()
<< std::endl;
  
  return 0;
}

Attachment: sphere_in_box_missing_wrongly_labeled_tetrahedra.png
Description: Binary data