Computing area of all facets in CGAL::Polyhedron_3

In this post I will show how to compute area of a facet in CGAL::Polyhedron_3.

ComputeFacetArea() Functor

ComputeFacetArea() is a thread-safe functor for computing the area of a given facet in a CGAL::Polyhedron_3. The facet must be a planar polygon with arbitrary number of sides. We need facet’s normal vector to compute it’s area. The facet normals must be be initialized using the ComputeFacetArea()’s constructor. The code for computing the facet normals is presented in this post: Computing normal of all facets in CGAL::Polyhedron_3.

#ifndef _SMESHLIB_OPERATIONS_COMPUTEFACETAREA_H_
#define _SMESHLIB_OPERATIONS_COMPUTEFACETAREA_H_

#include "PropertyMap.h"

namespace SMeshLib   {
namespace Operations {
;

// The ComputeFacetArea is a functor (delegate) to compute the area of a facet in CGAL::Polyhdeon_3.
// operator() is thread-safe.
// TPolyhedron is a type of CGAL::Polyhdeon_3.
// TFacetNormals is a property map associating a normal vector for each facet in the CGAL::Polyhdeon_3.
template<class TPolyhedron, class TFacetNormals>
struct ComputeFacetArea
{
public:
	
	// Redefine types from TPolyhedron for convenience.
	typedef typename TPolyhedron::Facet                                  Facet;
	typedef typename TPolyhedron::Traits::Vector_3                       Vector3;
	typedef typename TPolyhedron::Halfedge_around_facet_const_circulator HalfEdgeConstCirculator;
	
	// Return type of operator() required by QtConcurrent.
	typedef double result_type;
	
public:
	
	ComputeFacetArea(const TFacetNormals& facetNormals_)
		: facetNormals(facetNormals_)
	{}
	
	// Compute the area of a given facet.
	// The formula for computing facet area, which in general is a planar polygon, is from
	// http://softsurfer.com/Archive/algorithm_0101/algorithm_0101.htm
	inline double operator() (const Facet& f) const
	{
		std::vector<Vector3> _vectors;
		HalfEdgeConstCirculator s = f.facet_begin();
		HalfEdgeConstCirculator e = s;
		CGAL_For_all(s, e)
		{
			_vectors.push_back(s->vertex()->point() - CGAL::ORIGIN);
		}
		
		Vector3 _sumCrossProducts(CGAL::NULL_VECTOR);
		const size_t N = _vectors.size();
		for(size_t i=0 ; i<N ; ++i)
		{
			_sumCrossProducts = _sumCrossProducts + CGAL::cross_product(_vectors[i], _vectors[(i+1)%N]);
		}
		
		return fabs((facetNormals.value(&f) * _sumCrossProducts) / 2.0);
	}
	
public:
	
	// Facet normals are required for computing the area of a facet.
	const TFacetNormals& facetNormals;
};

};	// End namespace Operations.
};	// End namespace SMeshLib.

#endif // _SMESHLIB_OPERATIONS_COMPUTEFACETAREA_H_

Using ComputeFacetArea() Functor

Area of a facet f can be computed as double area = ComputeFacetArea(h);.

For most purposes, it is better to compute area of all facets once and cache them for later use. It is best to store the results in an associative container which associates the facet handle with the area. In the following example, I use PropertyMap which is a wrapper for std::set.

#include "ImportOBJ.h"
#include "ComputeFacetArea.h"
#include "ComputeFacetNormal.h"
#include "PropertyMap.h"
#include "CGAL/Simple_cartesian.h"
#include "CGAL/Polyhedron_items_3.h"
#include "CGAL/HalfedgeDS_list.h"
#include "CGAL/Polyhedron_3.h"

typedef CGAL::Simple_cartesian<double> Kernel;
typedef CGAL::Polyhedron_3<Kernel, 
	                       CGAL::Polyhedron_items_3, 
						   CGAL::HalfedgeDS_list> CgalPolyhedron;

// Compute the area of all facets in a CGAL::Polyhedron_3 and stores it in 
// facetAreas.
// TPolyhedron is a type of CGAL::Polyhedron_3.
// TFacetAreas is a property map which associates the facet area to facet in 
// the CGAL::Polyhdeon_3.
// TFacetNormals is a property map which associates a normal vector to facet in 
// the CGAL::Polyhdeon_3.
template<class TPolyhedron, class TFacetAreas, class TFacetNormals>
void computeFacetAreas(const TPolyhedron& polyhedron, TFacetAreas* facetAreas, const TFacetNormals& facetNormals)
{
	if(facetAreas == 0)
	{
		return;
	}
	
	typename TPolyhedron::Facet_const_iterator _begin = polyhedron.facets_begin();
	typename TPolyhedron::Facet_const_iterator _end   = polyhedron.facets_end();
	
	SMeshLib::Operations::ComputeFacetArea<TPolyhedron, TFacetNormals> _computeFacetArea(facetNormals);
	CGAL_For_all(_begin, _end)
	{
		facetAreas->setValue(&*_begin, _computeFacetArea(*_begin));
	}
}

// Compute the normal of all facets in a CGAL::Polyhedron_3 and stores it in facetNormals.
// TPolyhedron is a type of CGAL::Polyhedron_3.
// TFacetNormals is a property map which associates a normal vector to facet in the CGAL::Polyhdeon_3. 
template<class TPolyhedron, class TFacetNormals>
void computeFacetNormals(const TPolyhedron& polyhedron, TFacetNormals* facetNormals)
{
	if(facetNormals == 0)
	{
		return;
	}
	
	typename TPolyhedron::Facet_const_iterator _begin = polyhedron.facets_begin();
	typename TPolyhedron::Facet_const_iterator _end   = polyhedron.facets_end();

	SMeshLib::Operations::ComputeFacetNormal<TPolyhedron> _computeFacetNormal;
	CGAL_For_all(_begin, _end)
	{
		facetNormals->setValue(&*_begin, _computeFacetNormal(*_begin));
	}
}

void testComputeFacetArea()
{
	CgalPolyhedron _poly;
	SMeshLib::IO::importOBJ("Venus.obj", &_poly);
	
	typedef SMeshLib::PropertyMap<const CgalPolyhedron::Facet*, CgalPolyhedron::Traits::Vector_3> FacetNormalPM;
	FacetNormalPM _facetNormals;
	computeFacetNormals(_poly, &_facetNormals);
	
	typedef SMeshLib::PropertyMap<const CgalPolyhedron::Facet*, double> FacetAreaPM;
	FacetAreaPM _facetAreas;
	computeFacetAreas(_poly, &_facetAreas, _facetNormals);
}

Downloads

ImportOBJ.h
PropertyMap.h
ComputeFacetNormal.h
ComputeFacetArea.h
TestComputeFacetArea.cpp
Venus.obj
ComputeFacetArea.zip