#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_