CGAL users discussion list

["Lehtonen, Matti/HIIT" <>]

Quoting Stephane Tayeb 
<>:

>
> Hi,
>
> As APIs of C2T3 and C3T3 are not so different, it is probably not too
> hard to convert the C2T3 builder.
>
> Regards,
> Stéphane.

I managed to write something that almost works (see attachment). The problem 
is
that some facets are ignored and not included to surface. Any ideas, why it
could fail? Some facets are systematically ignored when queried by edge?


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_2_in_triangulation_3_polyhedron_builder.h $
// $Id: Complex_3_in_triangulation_3_polyhedron_builder.h,v 1.1.2.5 2010-07-26 13:14:16 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;
	}


	void
	construct_surface_from_faces
	(
		CGAL::Polyhedron_incremental_builder_3<HDS> &	builder,
		Tr const &										tr,
		std::set<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::set<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() )
				{
					size_t const	count = V.size();

					V[ vh ] = count;
				}
			}
		}

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

		// Add points
		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();

			builder.add_vertex( p );
		}

		for( typename std::set<Facet>::const_iterator fit = oriented_set.begin(); 
			 fit != oriented_set.end(); ++fit )
		{
			size_t		indices[ 3 ];
			int		index = 0;
			for( int i = 0; i < 3; i++ )
			{
				indices[ index++ ] = V[ fit->first->vertex( tr.vertex_triple_index( fit->second, i ) ) ];
			}
			builder.begin_facet();
			builder.add_vertex_to_facet( indices[ 0 ] );
			builder.add_vertex_to_facet( regular_orientation ? indices[ 1 ] : indices[ 2 ] );
			builder.add_vertex_to_facet( regular_orientation ? indices[ 2 ] : indices[ 1 ] );
			builder.end_facet();
		}
		builder.end_surface();
	}


	bool
	select_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;
	}


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::set<Facet> oriented_set;
		std::stack<Facet> stack;
		for( Facet_iterator fit = c3t3.facets_begin(); fit != c3t3.facets_end(); ++fit )
		{
			// Already processed?
			if( ( processed_set.find( *fit ) != processed_set.end() ) ||
				( processed_set.find( tr.mirror_facet( *fit ) ) != processed_set.end() ) )
			{
				continue;
			}

			// Surface selection
			if( !select_face( subdomain, surface, *fit ) )
			{
				continue;
			}

			// Locate all faces that are some way connected to this face
			processed_set.insert( *fit );
			stack.push( *fit );
			while( !stack.empty() ) 
			{
				Facet f = stack.top();
				stack.pop();
				oriented_set.insert( 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;

						// Already processed?
						if( ( processed_set.find( fn ) != processed_set.end() ) ||
							( processed_set.find( tr.mirror_facet( fn ) ) != processed_set.end() ) )
						{
							continue;
						}

						// Surface selection
						if( !select_face( subdomain, surface, fn ) ||
							( fn.first->is_facet_on_surface( fn.second ) != true ) )
						{
							continue;
						}

						processed_set.insert( fn );
						stack.push( fn );
						oriented_set.insert( fn );
					}
				}
			}

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

		// Most likely not needed, but be on safe side
		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