CGAL users discussion list

["Lehtonen, Matti/HIIT" <>]

Hi!

This seems to be now fully working, without complaining any skipped triangles.

Feel free to improve and import into CGAL.

Lehtonen, Matti

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

// Copyright (c) 2003-2007  INRIA Sophia-Antipolis (France).
// Copyright (c) 2008       GeometryFactory (France)
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org); you may redistribute it under
// the terms of the Q Public License version 1.0.
// See the file LICENSE.QPL distributed with CGAL.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL: svn+ssh://scm.gforge.inria.fr/svn/cgal/branches/CGAL-3.6-branch/Surface_mesher/include/CGAL/IO/Complex_3_in_triangulation_3_polyhedron_builder.h $
// $Id: Complex_3_in_triangulation_3_polyhedron_builder.h,v 1.1.2.10 2010-07-29 23:39:52 lehtomat Exp $
//
//
// Author(s)     : Matti Lehtonen, based on work of Laurent Rineau

#ifndef CGAL_IO_COMPLEX_3_IN_TRIANGULATION_3_POLYHEDRON_BUILDER_H
#define CGAL_IO_COMPLEX_3_IN_TRIANGULATION_3_POLYHEDRON_BUILDER_H

#include <CGAL/Modifier_base.h>
#include <CGAL/Polyhedron_incremental_builder_3.h>

namespace CGAL {

template <class C3T3, class Polyhedron_>
class Complex_3_in_triangulation_3_polyhedron_builder
  : public CGAL::Modifier_base<typename Polyhedron_::HalfedgeDS>
{
public:
	typedef C3T3 C3t3;
	typedef typename C3t3::Subdomain_index Subdomain_index;
	typedef Polyhedron_ Polyhedron;
	typedef typename Polyhedron::HalfedgeDS HDS;

private:
	Subdomain_index const	subdomain;
	Subdomain_index const	surface;
	C3t3 const &			c3t3;

	typedef CGAL::Modifier_base<typename Polyhedron::HalfedgeDS> Base;

	typedef typename C3t3::Triangulation Tr;
	typedef typename C3t3::Facet_iterator Facet_iterator;

	typedef typename Tr::Facet Facet;
	typedef typename Tr::Vertex_handle Vertex_handle;
	typedef typename Tr::Geom_traits::Vector_3 Vector;


	bool
	define_orientation
	(
		Tr const &		tr,
		Facet_iterator	first,
		Facet_iterator	end
	)
	{
		// Orients the whole mesh towards outside:
		// - find the facet with max z	
		Facet	top_facet = *first;
		double	top_z = 
		  ( top_facet.first->vertex( tr.vertex_triple_index( top_facet.second, 0 ) )->point().z() +
			top_facet.first->vertex( tr.vertex_triple_index( top_facet.second, 1 ) )->point().z() +
			top_facet.first->vertex( tr.vertex_triple_index( top_facet.second, 2 ) )->point().z() ) / 3.;
		for( Facet_iterator fit = first; fit != end; ++fit )
		{
			double	z = 
			  ( fit->first->vertex( tr.vertex_triple_index( fit->second, 0 ) )->point().z() +
				fit->first->vertex( tr.vertex_triple_index( fit->second, 1 ) )->point().z() +
				fit->first->vertex( tr.vertex_triple_index( fit->second, 2 ) )->point().z() ) / 3.;
			if( top_z < z )
			{
				top_facet = *fit;
				top_z = z;
			}
		}

		//	- orient the facet with max z towards +Z axis
		static Vector const		Z( 0, 0, 1 );
		Vertex_handle const		v0 = top_facet.first->vertex( tr.vertex_triple_index( top_facet.second, 0 ) );
		Vertex_handle const		v1 = top_facet.first->vertex( tr.vertex_triple_index( top_facet.second, 1 ) );
		Vertex_handle const		v2 = top_facet.first->vertex( tr.vertex_triple_index( top_facet.second, 2 ) );
		Vector const			normal = cross_product( v1->point() - v0->point(), 
														v2->point() - v1->point() );
		bool const				regular_orientation = ( Z * normal >= 0 );

		return regular_orientation;
	}


	size_t
	construct_surface_from_faces
	(
		CGAL::Polyhedron_incremental_builder_3<HDS> &	builder,
		Tr const &										tr,
		std::vector<Facet> const &						oriented_set,
		bool const										regular_orientation
	)
	{
		typedef typename Tr::Vertex Vertex;
		typedef typename Tr::Point Point;
	
		// Map finite vertices of faces
		std::map<Vertex_handle, size_t>		V;
		for( typename std::vector<Facet>::const_iterator fi = oriented_set.begin();
			 fi != oriented_set.end(); ++fi )
		{
			for( unsigned i = 0; i < 3; ++i )
			{
				Vertex_handle const &	vh = fi->first->vertex( tr.vertex_triple_index( fi->second, i ) );

				if( V.find( vh ) == V.end() )
				{
					// Update later
					V[ vh ] = 0;
				}
			}
		}

		// Begin surface
		builder.begin_surface( V.size(), oriented_set.size() );

		// Add points
		size_t	counter = 0;
		for( typename std::map<Vertex_handle, size_t>::const_iterator vi = V.begin();
			 vi != V.end(); ++vi )
		{
			Vertex_handle const &	vh = vi->first;
			Point const &			p  = vh->point();

			// Update now, in correct order
			V[ vh ] = counter++;

			builder.add_vertex( p );
		}

		size_t	skipped = 0;
		for( typename std::vector<Facet>::const_iterator fit = oriented_set.begin(); 
			 fit != oriented_set.end(); ++fit )
		{
			size_t		indices[ 3 ];
			for( int i = 0; i < 3; i++ )
			{
				indices[ i ] = V[ fit->first->vertex( tr.vertex_triple_index( fit->second, i ) ) ];
			}
			
			// Test
			std::vector<size_t>		list;
			list.push_back( indices[ 0 ] );
			list.push_back( regular_orientation ? indices[ 1 ] : indices[ 2 ] );
			list.push_back( regular_orientation ? indices[ 2 ] : indices[ 1 ] );
			list.push_back( indices[ 0 ] );

			if( builder.test_facet_indices( list ) )
			{
				builder.begin_facet();
				builder.add_vertex_to_facet( list[ 0 ] );
				builder.add_vertex_to_facet( list[ 1 ] );
				builder.add_vertex_to_facet( list[ 2 ] );
				builder.end_facet();
			}
			else
			{
#if 0
				builder.begin_facet();
				builder.add_vertex_to_facet( list[ 2 ] );
				builder.add_vertex_to_facet( list[ 1 ] );
				builder.add_vertex_to_facet( list[ 0 ] );
				builder.end_facet();
#else
				// Test again, counterwise
				list.clear();
				list.push_back( indices[ 0 ] );
				list.push_back( regular_orientation ? indices[ 2 ] : indices[ 1 ] );
				list.push_back( regular_orientation ? indices[ 1 ] : indices[ 2 ] );
				list.push_back( indices[ 0 ] );
				if( builder.test_facet_indices( list ) )
				{
					builder.begin_facet();
					builder.add_vertex_to_facet( list[ 0 ] );
					builder.add_vertex_to_facet( list[ 1 ] );
					builder.add_vertex_to_facet( list[ 2 ] );
					builder.end_facet();
				}
				else
				{
					// Crap!
					++skipped;
				}
#endif
			}
		}
		builder.end_surface();
		
		return skipped;
	}


	// Maybe this function should be a template parameter or virtual function
	inline
	bool
	accept_face
	(
		Subdomain_index const	subdomain,
		Subdomain_index const	surface,
		Facet const &			f
	)
	{
		if( ( c3t3.subdomain_index( f.first ) == subdomain ) &&
			( ( c3t3.surface_index( f ).first == surface ) ||
			  ( c3t3.surface_index( f ).second == surface ) ) )
		{
			return true;
		}

		return false;
	}

	inline
	bool
	on_selection
	(
		std::set<Facet> const &		processed_set,
		Facet const &				f
	)
	{
		// Surface selection
		// Already processed?
		return	accept_face( subdomain, surface, f ) && 
			///( f.first->is_facet_on_surface( f.second ) ) &&
				( processed_set.find( f ) == processed_set.end() );
	}

	inline
	void
	do_selection
	(
		std::stack<Facet> &		stack,
		std::set<Facet> &		processed_set,
		Facet const &			f,
		Facet const &			fm
	)
	{
		if( on_selection( processed_set, f ) )
		{
			processed_set.insert( f );
			stack.push( f );
		}
		else if( on_selection( processed_set, fm ) ) 
		{
			processed_set.insert( fm );
			stack.push( fm );
		}	
	}

public:
	Complex_3_in_triangulation_3_polyhedron_builder
	(
		Subdomain_index const	subdomain,
		Subdomain_index const	surface,
		C3t3 const &			c3t3
	)
	: Base(), subdomain( subdomain ), surface( surface ), c3t3( c3t3 )
	{}

	void 
	operator()
	(
		HDS&	hds
	) 
	{
		typedef typename Tr::Edge Edge;
		typedef typename Tr::Facet_circulator Facet_circulator;

		CGAL::Polyhedron_incremental_builder_3<HDS> builder( hds, true );
		Tr const &	tr = c3t3.triangulation();

		bool const	regular_orientation = define_orientation( tr, 
															  c3t3.facets_begin(), 
															  c3t3.facets_end() );

		std::set<Facet>		processed_set;
		std::vector<Facet>	oriented_set; // Ordered set
		std::stack<Facet>	stack;
		size_t				rejected = 0;
		for( Facet_iterator fit = c3t3.facets_begin(); fit != c3t3.facets_end(); ++fit )
		{
			do_selection( stack, processed_set, 
						  *fit, tr.mirror_facet( *fit ) );

			// Locate all faces that are some way connected to faces in stack
			while( !stack.empty() ) 
			{
				Facet f = stack.top();
				stack.pop();
				oriented_set.push_back( f );

				// Process neighbor faces of f
				for( int ih = 0 ; ih < 3 ; ++ih ) 
				{
					int const			i1  = tr.vertex_triple_index( f.second, tr. cw( ih ) );
					int const			i2  = tr.vertex_triple_index( f.second, tr.ccw( ih ) );
					Facet_circulator	fc = tr.incident_facets( Edge( f.first, i1, i2 ), f );
					Facet_circulator	done( fc );

					CGAL_For_all( fc, done )
					{
						Facet const & fn = *fc;

						do_selection( stack, processed_set, 
									  fn, tr.mirror_facet( fn ) );
					}
				}
			}

			// Process the oriented set
			rejected += construct_surface_from_faces( builder, tr, oriented_set, regular_orientation );
			oriented_set.clear();

		}

		if( rejected )
		{
			std::cerr << "Warning: skipped " << rejected 
					  << " triangle" <<  ( rejected == 1 ? "" : "s"  ) << std::endl;
			builder.remove_unconnected_vertices();
		}
	} // end operator()
}; // end Complex_3_in_triangulation_3_polyhedron_builder

} // end namespace CGAL

#endif  // CGAL_IO_COMPLEX_3_IN_TRIANGULATION_3_POLYHEDRON_BUILDER_H