normal.hpp 11.7 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 "Algo/Geometry/basic.h"
#include "Algo/Geometry/area.h"
Pierre Kraemer's avatar
Pierre Kraemer committed
27

28
#include "Topology/generic/traversorCell.h"
Pierre Kraemer's avatar
Pierre Kraemer committed
29 30
#include "Topology/generic/traversor2.h"

31 32
#include "Algo/Parallel/parallel_foreach.h"

Pierre Kraemer's avatar
Pierre Kraemer committed
33 34 35 36 37 38 39 40
#include <cmath>

namespace CGoGN
{

namespace Algo
{

41 42 43
namespace Surface
{

Pierre Kraemer's avatar
Pierre Kraemer committed
44 45 46
namespace Geometry
{

47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68
//template <typename PFP>
//typename PFP::VEC3 triangleNormal(typename PFP::MAP& map, Dart d, const VertexAttribute<typename PFP::VEC3>& position)
//{
//	typename PFP::VEC3 N = Geom::triangleNormal(position[d], position[map.phi1(d)], position[map.phi_1(d)]) ;
//	N.normalize() ;
//	return N ;
//}

//template<typename PFP>
//typename PFP::VEC3 newellNormal(typename PFP::MAP& map, Dart d, const VertexAttribute<typename PFP::VEC3>& position)
//{
//	typename PFP::VEC3 N(0);

//	Traversor2FV<typename PFP::MAP> t(map, d) ;
//	for(Dart it = t.begin(); it != t.end(); it = t.next())
//	{
//		const typename PFP::VEC3& P = position[it];
//		const typename PFP::VEC3& Q = position[map.phi1(it)];
//		N[0] += (P[1] - Q[1]) * (P[2] + Q[2]);
//		N[1] += (P[2] - Q[2]) * (P[0] + Q[0]);
//		N[2] += (P[0] - Q[0]) * (P[1] + Q[1]);
//	}
Pierre Kraemer's avatar
Pierre Kraemer committed
69

70 71 72
//	N.normalize();
//	return N;
//}
Pierre Kraemer's avatar
Pierre Kraemer committed
73

74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118
//template <typename PFP>
//typename PFP::VEC3 faceNormal(typename PFP::MAP& map, Dart d, const VertexAttribute<typename PFP::VEC3>& position)
//{
//	if(map.faceDegree(d) == 3)
//		return triangleNormal<PFP>(map, d, position) ;
//	else
//		return newellNormal<PFP>(map, d, position) ;
//}


//template <typename PFP>
//typename PFP::VEC3 vertexNormal(typename PFP::MAP& map, Dart d, const VertexAttribute<typename PFP::VEC3>& position)
//{
//	typedef typename PFP::VEC3 VEC3 ;

//	VEC3 N(0) ;

//	Traversor2VF<typename PFP::MAP> t(map, d) ;
//	for(Dart it = t.begin(); it != t.end(); it = t.next())
//	{
//		VEC3 n = faceNormal<PFP>(map, it, position) ;
//		if(!n.hasNan())
//		{
//			VEC3 v1 = vectorOutOfDart<PFP>(map, it, position) ;
//			VEC3 v2 = vectorOutOfDart<PFP>(map, map.phi_1(it), position) ;
//			n *= convexFaceArea<PFP>(map, it, position) / (v1.norm2() * v2.norm2()) ;
//			N += n ;
//		}
//	}

//	N.normalize() ;
//	return N ;
//}


//template <typename PFP>
//typename PFP::VEC3 vertexBorderNormal(typename PFP::MAP& map, Dart d, const VertexAttribute<typename PFP::VEC3>& position)
//{
//	assert(map.dimension() == 3);

//	typedef typename PFP::VEC3 VEC3 ;

//	VEC3 N(0) ;
//	std::vector<Dart> faces;
//	CellMarker<FACE> f(map);
119

120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145
//	FunctorStore fs(faces);
//	map.foreach_dart_of_vertex(d,fs);

//	for(std::vector<Dart>::iterator it = faces.begin() ; it != faces.end() ; ++it)
//	{
//		if(!f.isMarked(*it) && map.isBoundaryFace(*it))
//		{
//			f.mark(*it);
//			VEC3 n = faceNormal<PFP>(map, *it, position);
//			if(!n.hasNan())
//			{
//				VEC3 v1 = vectorOutOfDart<PFP>(map, *it, position);
//				VEC3 v2 = vectorOutOfDart<PFP>(map, map.phi_1(*it), position);
//				n *= convexFaceArea<PFP>(map, *it, position) / (v1.norm2() * v2.norm2());
//				N += n ;
//			}
//		}
//	}

//	N.normalize() ;
//	return N ;
//}


template<typename PFP, typename V_ATT>
typename V_ATT::DATA_TYPE triangleNormal(typename PFP::MAP& map, Dart d, const V_ATT& position)
Pierre Kraemer's avatar
Pierre Kraemer committed
146
{
147
	typename V_ATT::DATA_TYPE N = Geom::triangleNormal(position[d], position[map.phi1(d)], position[map.phi_1(d)]) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
148 149 150 151
	N.normalize() ;
	return N ;
}

152 153
template<typename PFP, typename V_ATT>
typename V_ATT::DATA_TYPE newellNormal(typename PFP::MAP& map, Dart d, const V_ATT& position)
154
{
155
	typedef typename V_ATT::DATA_TYPE VEC3;
156 157 158 159 160 161 162 163 164 165 166 167 168 169
    VEC3 N(0);

    Traversor2FV<typename PFP::MAP> t(map, d) ;
    for(Dart it = t.begin(); it != t.end(); it = t.next())
    {
        const VEC3& P = position[it];
        const VEC3& Q = position[map.phi1(it)];
        N[0] += (P[1] - Q[1]) * (P[2] + Q[2]);
        N[1] += (P[2] - Q[2]) * (P[0] + Q[0]);
        N[2] += (P[0] - Q[0]) * (P[1] + Q[1]);
    }

    N.normalize();
    return N;
170 171
}

172 173
template<typename PFP, typename V_ATT>
typename V_ATT::DATA_TYPE faceNormal(typename PFP::MAP& map, Dart d, const V_ATT& position)
Pierre Kraemer's avatar
Pierre Kraemer committed
174
{
175
	if(map.faceDegree(d) == 3)
Pierre Kraemer's avatar
Pierre Kraemer committed
176 177
		return triangleNormal<PFP>(map, d, position) ;
	else
178
		return newellNormal<PFP>(map, d, position) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
179 180
}

181 182 183

template<typename PFP, typename V_ATT>
typename V_ATT::DATA_TYPE vertexNormal(typename PFP::MAP& map, Dart d, const V_ATT& position)
Pierre Kraemer's avatar
Pierre Kraemer committed
184
{
185
	typedef typename V_ATT::DATA_TYPE VEC3 ;
Pierre Kraemer's avatar
Pierre Kraemer committed
186 187

	VEC3 N(0) ;
188 189 190

	Traversor2VF<typename PFP::MAP> t(map, d) ;
	for(Dart it = t.begin(); it != t.end(); it = t.next())
Pierre Kraemer's avatar
Pierre Kraemer committed
191 192
	{
		VEC3 n = faceNormal<PFP>(map, it, position) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
193 194 195 196 197 198 199
		if(!n.hasNan())
		{
			VEC3 v1 = vectorOutOfDart<PFP>(map, it, position) ;
			VEC3 v2 = vectorOutOfDart<PFP>(map, map.phi_1(it), position) ;
			n *= convexFaceArea<PFP>(map, it, position) / (v1.norm2() * v2.norm2()) ;
			N += n ;
		}
200 201
	}

Pierre Kraemer's avatar
Pierre Kraemer committed
202 203 204 205
	N.normalize() ;
	return N ;
}

206 207
template<typename PFP, typename V_ATT>
typename V_ATT::DATA_TYPE vertexBorderNormal(typename PFP::MAP& map, Dart d, const V_ATT& position)
208 209 210
{
	assert(map.dimension() == 3);

211
	typedef typename V_ATT::DATA_TYPE VEC3 ;
212 213 214

	VEC3 N(0) ;
	std::vector<Dart> faces;
Pierre Kraemer's avatar
Pierre Kraemer committed
215
	CellMarker<typename PFP::MAP, FACE> f(map);
216 217

	FunctorStore fs(faces);
218
	map.foreach_dart_of_vertex(d,fs);
219 220 221

	for(std::vector<Dart>::iterator it = faces.begin() ; it != faces.end() ; ++it)
	{
222
		if(!f.isMarked(*it) && map.isBoundaryFace(*it))
223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239
		{
			f.mark(*it);
			VEC3 n = faceNormal<PFP>(map, *it, position);
			if(!n.hasNan())
			{
				VEC3 v1 = vectorOutOfDart<PFP>(map, *it, position);
				VEC3 v2 = vectorOutOfDart<PFP>(map, map.phi_1(*it), position);
				n *= convexFaceArea<PFP>(map, *it, position) / (v1.norm2() * v2.norm2());
				N += n ;
			}
		}
	}

	N.normalize() ;
	return N ;
}

240 241
template <typename PFP, typename V_ATT, typename F_ATT>
void computeNormalFaces(typename PFP::MAP& map, const V_ATT& position, F_ATT& face_normal, unsigned int thread)
Pierre Kraemer's avatar
Pierre Kraemer committed
242
{
243
	TraversorF<typename PFP::MAP> trav(map, thread);
244
	for (Dart d = trav.begin(); d != trav.end(); d = trav.next())
Pierre Kraemer's avatar
Pierre Kraemer committed
245
		face_normal[d] = faceNormal<PFP>(map, d, position) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
246 247
}

248 249
template <typename PFP, typename V_ATT>
void computeNormalVertices(typename PFP::MAP& map, const V_ATT& position, V_ATT& normal, unsigned int thread)
Pierre Kraemer's avatar
Pierre Kraemer committed
250
{
251
	TraversorV<typename PFP::MAP> trav(map, thread);
252
	for (Dart d = trav.begin(); d != trav.end(); d = trav.next())
Pierre Kraemer's avatar
Pierre Kraemer committed
253
		normal[d] = vertexNormal<PFP>(map, d, position) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
254 255
}

256 257 258 259 260


namespace Parallel
{

261
template <typename PFP, typename V_ATT>
262 263
class FunctorComputeNormalVertices: public FunctorMapThreaded<typename PFP::MAP >
{
264 265
	 const V_ATT& m_position;
	 V_ATT& m_normal;
266
public:
267
	 FunctorComputeNormalVertices<PFP,V_ATT>(	typename PFP::MAP& map, const V_ATT& position, V_ATT& normal):
268 269 270
	 	 FunctorMapThreaded<typename PFP::MAP>(map), m_position(position), m_normal(normal)
	 { }

271
	void run(Dart d, unsigned int threadID)
272 273 274 275 276
	{
		m_normal[d] = vertexNormal<PFP>(this->m_map, d, m_position) ;
	}
};

277 278
template <typename PFP, typename V_ATT>
void computeNormalVertices(typename PFP::MAP& map, const V_ATT& position, V_ATT& normal, unsigned int nbth)
279
{
280
	FunctorComputeNormalVertices<PFP,V_ATT> funct(map,position,normal);
281
	Algo::Parallel::foreach_cell<typename PFP::MAP,VERTEX>(map, funct, nbth, false);
282 283 284
}


285
template <typename PFP, typename V_ATT, typename F_ATT>
286 287
class FunctorComputeNormalFaces: public FunctorMapThreaded<typename PFP::MAP >
{
288 289
	 const V_ATT& m_position;
	 F_ATT& m_normal;
290
public:
291
	 FunctorComputeNormalFaces<PFP,V_ATT,F_ATT>( typename PFP::MAP& map, const V_ATT& position, F_ATT& normal):
292 293 294
	 	 FunctorMapThreaded<typename PFP::MAP>(map), m_position(position), m_normal(normal)
	 { }

295
	void run(Dart d, unsigned int /*threadID*/)
296 297 298 299 300
	{
		m_normal[d] = faceNormal<PFP>(this->m_map, d, m_position) ;
	}
};

301 302
template <typename PFP, typename V_ATT, typename F_ATT>
void computeNormalFaces(typename PFP::MAP& map, const V_ATT& position, F_ATT& normal, unsigned int nbth)
303
{
304
	FunctorComputeNormalFaces<PFP,V_ATT,F_ATT> funct(map,position,normal);
305
	Algo::Parallel::foreach_cell<typename PFP::MAP,FACE>(map, funct, nbth, false);
306 307
}

308

309
template <typename PFP, typename V_ATT, typename E_ATT>
310 311
class FunctorComputeAngleBetweenNormalsOnEdge: public FunctorMapThreaded<typename PFP::MAP >
{
312 313
	 const V_ATT& m_position;
	 E_ATT& m_angles;
314
public:
315
	 FunctorComputeAngleBetweenNormalsOnEdge<PFP,V_ATT,E_ATT>( typename PFP::MAP& map, const V_ATT& position, E_ATT& angles):
316 317 318 319 320 321 322 323 324 325
	 	 FunctorMapThreaded<typename PFP::MAP>(map), m_position(position), m_angles(angles)
	 { }

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


326 327
template <typename PFP, typename V_ATT, typename E_ATT>
void computeAnglesBetweenNormalsOnEdges(typename PFP::MAP& map, const V_ATT& position, E_ATT& angles, unsigned int nbth)
328
{
329
	FunctorComputeAngleBetweenNormalsOnEdge<PFP,V_ATT,E_ATT> funct(map,position,angles);
330
	Algo::Parallel::foreach_cell<typename PFP::MAP,EDGE>(map, funct, nbth, false);
331 332
}

333 334
} // endnamespace Parallel

335 336


337 338
template <typename PFP, typename V_ATT>
typename PFP::REAL computeAngleBetweenNormalsOnEdge(typename PFP::MAP& map, Dart d, const V_ATT& position)
339
{
340
	typedef typename V_ATT::DATA_TYPE VEC3 ;
341

342
	if(map.isBoundaryEdge(d))
343 344
		return 0 ;

345
	Dart dd = map.phi2(d) ;
346 347 348 349 350 351 352
	const VEC3 n1 = faceNormal<PFP>(map, d, position) ;
	const VEC3 n2 = faceNormal<PFP>(map, dd, position) ;
	VEC3 e = position[dd] - position[d] ;
	e.normalize() ;
	typename PFP::REAL s = e * (n1 ^ n2) ;
	typename PFP::REAL c = n1 * n2 ;
	typename PFP::REAL a(0) ;
353

354
	// the following trick is useful for avoiding NaNs (due to floating point errors)
355
	if (c > 0.5) a = asin(s) ;
356 357
	else
	{
358 359 360
		if(c < -1) c = -1 ;
		if (s >= 0) a = acos(c) ;
		else a = -acos(c) ;
361
	}
Pierre Kraemer's avatar
Pierre Kraemer committed
362 363
//	if (isnan(a))
	if(a != a)
364
		std::cerr<< "Warning : computeAngleBetweenNormalsOnEdge returns NaN on edge " << d << "-" << dd << std::endl ;
365 366 367
	return a ;
}

368 369
template <typename PFP, typename V_ATT, typename E_ATT>
void computeAnglesBetweenNormalsOnEdges(typename PFP::MAP& map, const V_ATT& position, E_ATT& angles, unsigned int thread)
370
{
371
	TraversorE<typename PFP::MAP> trav(map, thread);
372
	for (Dart d = trav.begin(); d != trav.end(); d = trav.next())
Pierre Kraemer's avatar
Pierre Kraemer committed
373
		angles[d] = computeAngleBetweenNormalsOnEdge<PFP>(map, d, position) ;
374 375
}

Pierre Kraemer's avatar
Pierre Kraemer committed
376 377
} // namespace Geometry

378 379
}

Pierre Kraemer's avatar
Pierre Kraemer committed
380 381 382
} // namespace Algo

} // namespace CGoGN