Skip to Content.
Sympa Menu

cgal-discuss - [cgal-discuss] 3D Mesher: missing and mislabeled tetrahedra

Subject: CGAL users discussion list

List archive

[cgal-discuss] 3D Mesher: missing and mislabeled tetrahedra


Chronological Thread 
  • From: Maarten Moesen <>
  • To:
  • Subject: [cgal-discuss] 3D Mesher: missing and mislabeled tetrahedra
  • Date: Wed, 7 Nov 2012 16:09:55 +0100


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



  • [cgal-discuss] 3D Mesher: missing and mislabeled tetrahedra, Maarten Moesen, 11/07/2012

Archive powered by MHonArc 2.6.18.

Top of Page