area.hpp 9.35 KB
Newer Older
Pierre Kraemer's avatar
Pierre Kraemer committed
1 2 3
/*******************************************************************************
 * CGoGN: Combinatorial and Geometric modeling with Generic N-dimensional Maps  *
 * version 0.1                                                                  *
4
 * Copyright (C) 2009-2012, IGG Team, LSIIT, University of Strasbourg           *
Pierre Kraemer's avatar
Pierre Kraemer committed
5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
 *                                                                              *
 * This library is free software; you can redistribute it and/or modify it      *
 * under the terms of the GNU Lesser General Public License as published by the *
 * Free Software Foundation; either version 2.1 of the License, or (at your     *
 * option) any later version.                                                   *
 *                                                                              *
 * This library is distributed in the hope that it will be useful, but WITHOUT  *
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or        *
 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License  *
 * for more details.                                                            *
 *                                                                              *
 * You should have received a copy of the GNU Lesser General Public License     *
 * along with this library; if not, write to the Free Software Foundation,      *
 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301 USA.           *
 *                                                                              *
20
 * Web site: http://cgogn.unistra.fr/                                           *
Pierre Kraemer's avatar
Pierre Kraemer committed
21 22 23 24 25 26
 * Contact information: cgogn@unistra.fr                                        *
 *                                                                              *
 *******************************************************************************/

#include "Geometry/basic.h"
#include "Algo/Geometry/centroid.h"
27 28
#include "Topology/generic/traversorCell.h"
#include "Topology/generic/traversor2.h"
Pierre Kraemer's avatar
Pierre Kraemer committed
29 30 31 32 33 34 35 36 37 38 39

namespace CGoGN
{

namespace Algo
{

namespace Geometry
{

template <typename PFP>
40
typename PFP::REAL triangleArea(typename PFP::MAP& map, Dart d, const VertexAttribute<typename PFP::VEC3>& position)
Pierre Kraemer's avatar
Pierre Kraemer committed
41
{
Pierre Kraemer's avatar
Pierre Kraemer committed
42 43 44
	typename PFP::VEC3 p1 = position[d] ;
	typename PFP::VEC3 p2 = position[map.phi1(d)] ;
	typename PFP::VEC3 p3 = position[map.phi_1(d)] ;
Pierre Kraemer's avatar
Pierre Kraemer committed
45 46 47 48 49

	return Geom::triangleArea(p1, p2, p3) ;
}

template <typename PFP>
50
typename PFP::REAL convexFaceArea(typename PFP::MAP& map, Dart d, const VertexAttribute<typename PFP::VEC3>& position)
Pierre Kraemer's avatar
Pierre Kraemer committed
51
{
Pierre Kraemer's avatar
Pierre Kraemer committed
52
	typedef typename PFP::VEC3 VEC3 ;
Pierre Kraemer's avatar
Pierre Kraemer committed
53

54
	if(map.faceDegree(d) == 3)
Pierre Kraemer's avatar
Pierre Kraemer committed
55 56 57 58 59
		return triangleArea<PFP>(map, d, position) ;
	else
	{
		float area = 0.0f ;
		VEC3 centroid = Algo::Geometry::faceCentroid<PFP>(map, d, position) ;
60 61
		Traversor2FE<typename PFP::MAP> t(map, d) ;
		for(Dart it = t.begin(); it != t.end(); it = t.next())
Pierre Kraemer's avatar
Pierre Kraemer committed
62
		{
Pierre Kraemer's avatar
Pierre Kraemer committed
63 64
			VEC3 p1 = position[it] ;
			VEC3 p2 = position[map.phi1(it)] ;
Pierre Kraemer's avatar
Pierre Kraemer committed
65
			area += Geom::triangleArea(p1, p2, centroid) ;
66
		}
Pierre Kraemer's avatar
Pierre Kraemer committed
67 68 69 70 71
		return area ;
	}
}

template <typename PFP>
72
typename PFP::REAL totalArea(typename PFP::MAP& map, const VertexAttribute<typename PFP::VEC3>& position, const FunctorSelect& select, unsigned int thread)
Pierre Kraemer's avatar
Pierre Kraemer committed
73
{
74
	typename PFP::REAL area(0) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
75
	TraversorF<typename PFP::MAP> t(map, select) ;
76
	for(Dart d = t.begin(); d != t.end(); d = t.next())
Pierre Kraemer's avatar
Pierre Kraemer committed
77
		area += convexFaceArea<PFP>(map, d, position) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
78 79 80
	return area ;
}

81
template <typename PFP>
82
typename PFP::REAL vertexOneRingArea(typename PFP::MAP& map, Dart d, const VertexAttribute<typename PFP::VEC3>& position)
83 84
{
	typename PFP::REAL area(0) ;
85 86
	Traversor2VF<typename PFP::MAP> t(map, d) ;
	for(Dart it = t.begin(); it != t.end(); it = t.next())
87 88 89 90 91
		area += convexFaceArea<PFP>(map, it, position) ;
	return area ;
}

template <typename PFP>
92
typename PFP::REAL vertexBarycentricArea(typename PFP::MAP& map, Dart d, const VertexAttribute<typename PFP::VEC3>& position)
93 94
{
	typename PFP::REAL area(0) ;
95 96
	Traversor2VF<typename PFP::MAP> t(map, d) ;
	for(Dart it = t.begin(); it != t.end(); it = t.next())
97
		area += convexFaceArea<PFP>(map, it, position) / 3 ;
98 99 100
	return area ;
}

101
template <typename PFP>
102
typename PFP::REAL vertexVoronoiArea(typename PFP::MAP& map, Dart d, const VertexAttribute<typename PFP::VEC3>& position)
103 104
{
	typename PFP::REAL area(0) ;
105 106
	Traversor2VF<typename PFP::MAP> t(map, d) ;
	for(Dart it = t.begin(); it != t.end(); it = t.next())
107
	{
108 109 110 111
		const typename PFP::VEC3& p1 = position[it] ;
		const typename PFP::VEC3& p2 = position[map.phi1(it)] ;
		const typename PFP::VEC3& p3 = position[map.phi_1(it)] ;
		if(!Geom::isTriangleObtuse(p1, p2, p3))
112
		{
113 114 115
			typename PFP::REAL a = Geom::angle(p3 - p2, p1 - p2) ;
			typename PFP::REAL b = Geom::angle(p1 - p3, p2 - p3) ;
			area += ( (p2 - p1).norm2() / tan(b) + (p3 - p1).norm2() / tan(a) ) / 8 ;
116 117 118
		}
		else
		{
119 120
			typename PFP::REAL tArea = Geom::triangleArea(p1, p2, p3) ;
			if(Geom::angle(p2 - p1, p3 - p1) > M_PI / 2)
121 122 123 124
				area += tArea / 2 ;
			else
				area += tArea / 4 ;
		}
125
	}
126 127 128
	return area ;
}

Pierre Kraemer's avatar
Pierre Kraemer committed
129
template <typename PFP>
130
void computeAreaFaces(typename PFP::MAP& map, const VertexAttribute<typename PFP::VEC3>& position, FaceAttribute<typename PFP::REAL>& face_area, const FunctorSelect& select)
Pierre Kraemer's avatar
Pierre Kraemer committed
131
{
Pierre Kraemer's avatar
Pierre Kraemer committed
132
	TraversorF<typename PFP::MAP> t(map, select) ;
133
	for(Dart d = t.begin(); d != t.end(); d = t.next())
Pierre Kraemer's avatar
Pierre Kraemer committed
134
		face_area[d] = convexFaceArea<PFP>(map, d, position) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
135 136
}

137
template <typename PFP>
138
void computeOneRingAreaVertices(typename PFP::MAP& map, const VertexAttribute<typename PFP::VEC3>& position, VertexAttribute<typename PFP::REAL>& vertex_area, const FunctorSelect& select)
139
{
Pierre Kraemer's avatar
Pierre Kraemer committed
140
	TraversorV<typename PFP::MAP> t(map, select) ;
141
	for(Dart d = t.begin(); d != t.end(); d = t.next())
Pierre Kraemer's avatar
Pierre Kraemer committed
142
		vertex_area[d] = vertexOneRingArea<PFP>(map, d, position) ;
143 144
}

145
template <typename PFP>
146
void computeBarycentricAreaVertices(typename PFP::MAP& map, const VertexAttribute<typename PFP::VEC3>& position, VertexAttribute<typename PFP::REAL>& vertex_area, const FunctorSelect& select)
147
{
Pierre Kraemer's avatar
Pierre Kraemer committed
148
	TraversorV<typename PFP::MAP> t(map, select) ;
149
	for(Dart d = t.begin(); d != t.end(); d = t.next())
Pierre Kraemer's avatar
Pierre Kraemer committed
150
		vertex_area[d] = vertexBarycentricArea<PFP>(map, d, position) ;
151 152
}

153
template <typename PFP>
154
void computeVoronoiAreaVertices(typename PFP::MAP& map, const VertexAttribute<typename PFP::VEC3>& position, VertexAttribute<typename PFP::REAL>& vertex_area, const FunctorSelect& select)
155
{
Pierre Kraemer's avatar
Pierre Kraemer committed
156
	TraversorV<typename PFP::MAP> t(map, select) ;
157
	for(Dart d = t.begin(); d != t.end(); d = t.next())
Pierre Kraemer's avatar
Pierre Kraemer committed
158
		vertex_area[d] = vertexVoronoiArea<PFP>(map, d, position) ;
159 160
}

161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244

namespace Parallel
{
template <typename PFP>
class FunctorConvexFaceArea: public FunctorMapThreaded<typename PFP::MAP >
{
	 const VertexAttribute<typename PFP::VEC3>& m_position;
	 FaceAttribute<typename PFP::REAL>& m_area;
public:
	 FunctorConvexFaceArea<PFP>( typename PFP::MAP& map, const VertexAttribute<typename PFP::VEC3>& position, FaceAttribute<typename PFP::REAL>& area):
	 	 FunctorMapThreaded<typename PFP::MAP>(map), m_position(position), m_area(area)
	 { }

	void run(Dart d, unsigned int threadID)
	{
		m_area[d] = convexFaceArea<PFP>(this->m_map, d, m_position) ;
	}
};

template <typename PFP>
void computeAreaFaces(typename PFP::MAP& map, const VertexAttribute<typename PFP::VEC3>& position, FaceAttribute<typename PFP::REAL>& area, const FunctorSelect& select, unsigned int nbth, unsigned int current_thread)
{
	FunctorConvexFaceArea<PFP> funct(map,position,area);
	Algo::Parallel::foreach_cell<typename PFP::MAP,FACE>(map, funct, nbth, false, select, current_thread);
}


template <typename PFP>
class FunctorVertexOneRingArea: public FunctorMapThreaded<typename PFP::MAP >
{
	 const VertexAttribute<typename PFP::VEC3>& m_position;
	 VertexAttribute<typename PFP::REAL>& m_area;
public:
	 FunctorVertexOneRingArea<PFP>( typename PFP::MAP& map, const VertexAttribute<typename PFP::VEC3>& position, VertexAttribute<typename PFP::REAL>& area):
	 	 FunctorMapThreaded<typename PFP::MAP>(map), m_position(position), m_area(area)
	 { }

	void run(Dart d, unsigned int threadID)
	{
		m_area[d] = vertexOneRingArea<PFP>(this->m_map, d, m_position) ;
	}
};

template <typename PFP>
void computeOneRingAreaVertices(typename PFP::MAP& map, const VertexAttribute<typename PFP::VEC3>& position, VertexAttribute<typename PFP::REAL>& area, const FunctorSelect& select, unsigned int nbth, unsigned int current_thread)
{
	FunctorConvexFaceArea<PFP> funct(map,position,area);
	Algo::Parallel::foreach_cell<typename PFP::MAP,VERTEX>(map, funct, nbth, false, select, current_thread);
}



template <typename PFP>
class FunctorVertexVoronoiArea: public FunctorMapThreaded<typename PFP::MAP >
{
	 const VertexAttribute<typename PFP::VEC3>& m_position;
	 VertexAttribute<typename PFP::REAL>& m_area;
public:
	 FunctorVertexVoronoiArea<PFP>( typename PFP::MAP& map, const VertexAttribute<typename PFP::VEC3>& position, VertexAttribute<typename PFP::REAL>& area):
	 	 FunctorMapThreaded<typename PFP::MAP>(map), m_position(position), m_area(area)
	 { }

	void run(Dart d, unsigned int threadID)
	{
		m_area[d] = vertexVoronoiArea<PFP>(this->m_map, d, m_position) ;
	}
};

template <typename PFP>
void computeVoronoiAreaVertices(typename PFP::MAP& map, const VertexAttribute<typename PFP::VEC3>& position, VertexAttribute<typename PFP::REAL>& area, const FunctorSelect& select, unsigned int nbth, unsigned int current_thread)
{
	FunctorConvexFaceArea<PFP> funct(map,position,area);
	Algo::Parallel::foreach_cell<typename PFP::MAP,VERTEX>(map, funct, nbth, false, select, current_thread);
}


}







Pierre Kraemer's avatar
Pierre Kraemer committed
245 246 247 248 249
} // namespace Geometry

} // namespace Algo

} // namespace CGoGN