CGAL users discussion list

["Lehtonen, Matti/HIIT" <>]

Hi!

I have written a new mesh domain (labeled point set). How ever, it generates
invalid meshes (some triangles won't satisfy criterias like minimum angle of 
30
degrees and triangles are kind of 'soap'). I believe that the wrapper function
or CGAL::Mesh_3::Labeled_mesh_domain_3 aren't faulty, but the conversion from
C3T3 to HDS fails for some reason. Any help is more than welcome.


Function wrapper at CGAL::Point_set_to_labeled_function_wrapper (yes, I am
planning to donate this, if this works):
  return_type
  operator()
  (
    Point_3 const &    p,
    bool const = true
  )
  const
  {
    // Bounding box of kd -tree
    TreeBbox const &    bbox = _set.bounding_box();

    if( ( bbox.min_coord( 0 ) > p.x() ) || ( p.x() > bbox.max_coord( 0 ) ) ||
        ( bbox.min_coord( 1 ) > p.y() ) || ( p.y() > bbox.max_coord( 1 ) ) ||
        ( bbox.min_coord( 2 ) > p.z() ) || ( p.z() > bbox.max_coord( 2 ) ) )
    {
      // Outside of domain
      return CGAL::NEVER_CLASSIFIED;
    }

    // Initialize the search structure, and search nearest point
    Query_item const &     q = Query_item( Point_with_normal( p ) );
    Neighbor_search        search( _set, q );
    Labeled_point const &  n1 = search.begin()->first;

    double const projectionLengthOnNormal = ( p.x() - n1.position().x() ) *
n1.normal().x()
        + ( p.y() - n1.position().y() ) * n1.normal().y()
        + ( p.z() - n1.position().z() ) * n1.normal().z();
    if( projectionLengthOnNormal > 0.0 )
    {
      // Assumption: points are from surface with well defined normals
      // (pointing [mostly] up) and above of surface is AIR -domain
      return CGAL::ARTIFICAL_DOMAIN_AIR;
    }
    else
    {
      return n1.label();
    }
  }

I believe that actually the C3T3 to HDS conversion fails for some reason,
I just cannot Identify why ...

typedef struct BuilderInfo
{
  size_t          number_of_faces;
  PointIndex      index;
  HalfedgeDS *    hds;
  Pib *           builder;

  BuilderInfo()
  : number_of_faces( 0 ), index( PointIndex() ), hds( NULL ), builder( NULL )
  {}

  ~BuilderInfo()
  {
    delete builder;
    delete hds;
  }
} BuilderInfo;

typedef class Mesh_domain : public
CGAL::Labeled_point_set_mesh_domain_3<PointSet,TreeBoundingBox,
Neighbor_search,Kernel> {...}    Mesh_domain;
typedef CGAL::Mesh_triangulation_3<Mesh_domain>::type    Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr>    C3t3;
typedef CGAL::Polyhedron_incremental_builder_3<HalfedgeDS>    Pib;

static
void
create_polyhedrons
(
  PolyhedronList &    polyhedrons,
  C3t3 const * const  c3t3
)
{
  typedef std::map<SurfaceKeys,BuilderInfo*>  Mapping;
  Mapping                                     mapping;

  // GO thru all 2D faces: collect face and vertice counts
  std::cerr << "*** 2D face count and index vertices ***"  << std::endl;
  std::cerr << "* faces - " << c3t3->number_of_facets() << std::endl;
  std::cerr << "* vertices - " << c3t3->number_of_vertices() << std::endl;
  for( C3t3::Facet_iterator fi = c3t3->facets_begin();
        fi != c3t3->facets_end(); ++fi )
  {
    C3t3::Facet const &                         f = *fi;
    C3t3::Surface_index const &         si = c3t3->surface_index( f );

    C3t3::Cell_handle const &           ch = f.first;
    C3t3::Subdomain_index const &       sdi = c3t3->subdomain_index( ch );
    if( ( sdi == CGAL::NEVER_CLASSIFIED ) ||
        ( si.first == CGAL::NEVER_CLASSIFIED ) ||
        ( si.second == CGAL::NEVER_CLASSIFIED ) )
    {
      // Ignore domains and surfaces related to non-classified volume
      continue;
    }

    // Select Polyhedron incremental builder
    SurfaceKeys    k( sdi, si );
    if( mapping.find( k ) == mapping.end() )
    {
      mapping.insert( std::make_pair( k, new BuilderInfo() ) );
    }
    BuilderInfo *    bi = mapping[ k ];

    // Count faces
    bi->number_of_faces++;

    // Index points and count vertices
    C3t3::Triangulation::Cell const &   c = *ch;
    for( int vi = 0; vi < 4; ++vi )
    {
      if( vi == f.second )
      {
        continue;
      }
      C3t3::Vertex_handle const &            vh = c.vertex( vi );
      C3t3::Triangulation::Vertex const &    v = *vh;
      Point const &                          p = v.point();

      if( bi->index.find( p ) == bi->index.end() )
      {
        ///std::cerr << "* p - " << p.x() << " " << p.y() << " " << p.z() <<
std::endl;
        size_t const    counter = bi->index.size();

        bi->index[ p ] = counter;
      }
    }
  }


  // GO thru all 2D surfaces: add vertices for surfaces
  std::cerr << "*** Add vertices to surface ***"  << std::endl;
  for( Mapping::iterator ki = mapping.begin(); ki != mapping.end(); ++ki )
  {
    SurfaceKeys const &            k = ki->first;
    C3t3::Subdomain_index const &  sdi = k.first;
    C3t3::Surface_index const &    si = k.second;
    std::cerr << "* subdomain index - " << sdi << std::endl;
    std::cerr << "* surface index - pair<" << si.first << "," << si.second <<
">" << std::endl;

    BuilderInfo *    bi = ki->second;
    std::cerr << "* faces - " << bi->number_of_faces << std::endl;
    std::cerr << "* vertices - " << bi->index.size() << std::endl;

    bi->hds = new HalfedgeDS( bi->index.size(), 3 * bi->number_of_faces,
bi->number_of_faces );
    bi->builder = new Pib( *( bi->hds ) /*, true */ );
    bi->builder->begin_surface( bi->index.size(), bi->number_of_faces,
                                3 * bi->number_of_faces );

    for( PointIndex::iterator pi = bi->index.begin();
         pi != bi->index.end(); ++pi )
    {
      Point const &     p = pi->first;

      bi->builder->add_vertex( p );
    }
  }

  // GO thru all 2D faces: add faces to surfaces
  std::cerr << "*** Add 2D faces to surfaces ***"  << std::endl;
  std::cerr << "* faces - " << c3t3->number_of_facets() << std::endl;
  for( C3t3::Facet_iterator fi = c3t3->facets_begin();
       fi != c3t3->facets_end(); ++fi )
  {
    C3t3::Facet const &          f = *fi;
    C3t3::Cell_handle const &    ch = f.first;
    C3t3::Surface_index const    si = c3t3->surface_index( f );
    C3t3::Subdomain_index const  sdi = c3t3->subdomain_index( ch );
    if( ( sdi == CGAL::NEVER_CLASSIFIED ) ||
        ( si.first == CGAL::NEVER_CLASSIFIED ) ||
        ( si.second == CGAL::NEVER_CLASSIFIED ) )
    {
      // Ignore domains and surfaces related to non-classified volume
      continue;
    }
    ///std::cerr << "* face - " << sdi << ", pair<"<< si.first << "," <<
si.second << ">" << std::endl;

    // Select Polyhedron incremental builder
    SurfaceKeys    k( sdi, si );
    BuilderInfo *    bi = mapping[ k ];
    bi->builder->begin_facet();

    C3t3::Triangulation::Triangle const &  t = 
c3t3->triangulation().triangle( f
);
    int         v_index = 0;
    for( size_t i = 0; i < 3;
         ++i, v_index = c3t3->triangulation().ccw( v_index ) )
    {
      Point const &     p = t.vertex( v_index );
      size_t const      idx = bi->index[ p ];

      ///std::cerr << "* vertex - " << idx << std::endl;

      bi->builder->add_vertex_to_facet( idx );
    }

    bi->builder->end_facet();
  }

  // Close all surfaces
  std::cerr << "*** Close all surfaces ***"  << std::endl;
  for( Mapping::iterator ki = mapping.begin();
       ki != mapping.end(); ++ki )
  {
    std::cerr << "Surface:" << std::endl;
    SurfaceKeys const &    k = ki->first;
    C3t3::Subdomain_index const &    sdi = k.first;
    C3t3::Surface_index const &    si = k.second;
    std::cerr << "* surface - " << sdi << ", pair<"<< si.first << "," <<
si.second << ">" << std::endl;
    BuilderInfo *    bi = ki->second;

    bi->builder->end_surface();
    std::cerr << "* halfedges - " << bi->hds->size_of_halfedges() << 
std::endl;
    std::cerr << "* vertices - " << bi->hds->size_of_vertices() << std::endl;


    std::cerr << "Remove unconnected vertices:"  << std::endl;
    (void) bi->builder->remove_unconnected_vertices();
    std::cerr << "* halfedges - " << bi->hds->size_of_halfedges() << 
std::endl;
    std::cerr << "* vertices - " << bi->hds->size_of_vertices() << std::endl;


    // Fill Polyhedron
    std::cerr << "Polyhedron:" << std::endl;
    SurfacePolyhedron   SP;
    SP.first = k;
    Polyhedron &    P = SP.second;
    Build_3D_Mesh_triangulation<HalfedgeDS>    triangulation( *( bi->hds ) );
    P.delegate( triangulation );
    std::cerr << "* facets - " << P.size_of_facets() << std::endl;
    std::cerr << "* halfedges - " << P.size_of_halfedges() << std::endl;
    std::cerr << "* vertices - " << P.size_of_vertices() << std::endl;
    std::cerr << "* is_valid - " << ( P.is_valid( false ) ? "YES" : "NO" ) <<
std::endl;


    std::cerr << "Border of polyhedron:" << std::endl;
    P.normalize_border();
    std::cerr << "* border halfedges - " <<  P.size_of_border_halfedges() <<
std::endl;
    std::cerr << "* border edges - " <<  P.size_of_border_edges() << 
std::endl;
    std::cerr << "* border is_valid - " << ( P.normalized_border_is_valid() ?
"YES" : "NO" ) << std::endl;


    // Return Polyhedron
    if( !P.empty() )
    {
      polyhedrons.push_back( SP );
    }
  }

  std::cerr << "*** Polyhedron list ***" << std::endl;
  std::cerr << "* polyhedrons - " << polyhedrons.size() << std::endl;
}


Thx for any help or suggestions
Lehtonen, Matti

Researcher, head programmer - Helsinki Institute for Information Technology 
HIIT
http://www.hiit.fi/
-- 
Life is complex. It has real and imaginary parts.