normal.hpp 6.97 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"

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

namespace CGoGN
{

namespace Algo
{

namespace Geometry
{

template <typename PFP>
43
typename PFP::VEC3 triangleNormal(typename PFP::MAP& map, Dart d, const VertexAttribute<typename PFP::VEC3>& position)
Pierre Kraemer's avatar
Pierre Kraemer committed
44
{
45
	typename PFP::VEC3 N = Geom::triangleNormal(position[d], position[map.phi1(d)], position[map.phi_1(d)]) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
46 47 48 49
	N.normalize() ;
	return N ;
}

50
template<typename PFP>
51
typename PFP::VEC3 newellNormal(typename PFP::MAP& map, Dart d, const VertexAttribute<typename PFP::VEC3>& position)
52
{
Pierre Kraemer's avatar
Pierre Kraemer committed
53 54 55
	typename PFP::VEC3 N(0);

	Traversor2FV<typename PFP::MAP> t(map, d) ;
56
	for(Dart it = t.begin(); it != t.end(); it = t.next())
57
	{
Pierre Kraemer's avatar
Pierre Kraemer committed
58
		const typename PFP::VEC3& P = position[it];
59
		const typename PFP::VEC3& Q = position[map.phi1(it)];
Pierre Kraemer's avatar
Pierre Kraemer committed
60 61 62
		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]);
63 64
	}

Pierre Kraemer's avatar
Pierre Kraemer committed
65 66
	N.normalize();
	return N;
67 68
}

Pierre Kraemer's avatar
Pierre Kraemer committed
69
template <typename PFP>
70
typename PFP::VEC3 faceNormal(typename PFP::MAP& map, Dart d, const VertexAttribute<typename PFP::VEC3>& position)
Pierre Kraemer's avatar
Pierre Kraemer committed
71
{
72
	if(map.faceDegree(d) == 3)
Pierre Kraemer's avatar
Pierre Kraemer committed
73 74
		return triangleNormal<PFP>(map, d, position) ;
	else
75 76 77 78 79 80 81 82 83 84 85 86 87 88 89
		return newellNormal<PFP>(map, d, position) ;
//	{
//		VEC3 N(0) ;
//		Dart it = d ;
//		do
//		{
//			VEC3 n = triangleNormal<PFP>(map, it, position) ;
//			//if(!std::isnan(n[0]) && !std::isnan(n[1]) && !std::isnan(n[2]))
//			if(!n.hasNan())
//				N += n ;
//			it = map.phi1(it) ;
//		} while (it != d) ;
//		N.normalize() ;
//		return N ;
//	}
Pierre Kraemer's avatar
Pierre Kraemer committed
90 91 92
}

template <typename PFP>
93
typename PFP::VEC3 vertexNormal(typename PFP::MAP& map, Dart d, const VertexAttribute<typename PFP::VEC3>& position)
Pierre Kraemer's avatar
Pierre Kraemer committed
94 95 96 97
{
	typedef typename PFP::VEC3 VEC3 ;

	VEC3 N(0) ;
98 99 100

	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
101 102
	{
		VEC3 n = faceNormal<PFP>(map, it, position) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
103 104 105 106 107 108 109
		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 ;
		}
110 111
	}

Pierre Kraemer's avatar
Pierre Kraemer committed
112 113 114 115
	N.normalize() ;
	return N ;
}

116
template <typename PFP>
117
typename PFP::VEC3 vertexBorderNormal(typename PFP::MAP& map, Dart d, const VertexAttribute<typename PFP::VEC3>& position)
118 119 120 121 122 123 124
{
	assert(map.dimension() == 3);

	typedef typename PFP::VEC3 VEC3 ;

	VEC3 N(0) ;
	std::vector<Dart> faces;
Pierre Kraemer's avatar
Pierre Kraemer committed
125
	CellMarker<FACE> f(map);
126 127

	FunctorStore fs(faces);
128
	map.foreach_dart_of_vertex(d,fs);
129 130 131

	for(std::vector<Dart>::iterator it = faces.begin() ; it != faces.end() ; ++it)
	{
132
		if(!f.isMarked(*it) && map.isBoundaryFace(*it))
133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149
		{
			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 ;
}

Pierre Kraemer's avatar
Pierre Kraemer committed
150
template <typename PFP>
151
void computeNormalFaces(typename PFP::MAP& map, const VertexAttribute<typename PFP::VEC3>& position, FaceAttribute<typename PFP::VEC3>& face_normal, const FunctorSelect& select, unsigned int thread)
Pierre Kraemer's avatar
Pierre Kraemer committed
152
{
Pierre Kraemer's avatar
Pierre Kraemer committed
153
	TraversorF<typename PFP::MAP> trav(map, select, thread);
154
	for (Dart d = trav.begin(); d != trav.end(); d = trav.next())
Pierre Kraemer's avatar
Pierre Kraemer committed
155
		face_normal[d] = faceNormal<PFP>(map, d, position) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
156 157 158
}

template <typename PFP>
159
void computeNormalVertices(typename PFP::MAP& map, const VertexAttribute<typename PFP::VEC3>& position, VertexAttribute<typename PFP::VEC3>& normal, const FunctorSelect& select, unsigned int thread)
Pierre Kraemer's avatar
Pierre Kraemer committed
160
{
Pierre Kraemer's avatar
Pierre Kraemer committed
161
	TraversorV<typename PFP::MAP> trav(map, select, thread);
162
	for (Dart d = trav.begin(); d != trav.end(); d = trav.next())
Pierre Kraemer's avatar
Pierre Kraemer committed
163
		normal[d] = vertexNormal<PFP>(map, d, position) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
164 165
}

166
template <typename PFP>
167
typename PFP::REAL computeAngleBetweenNormalsOnEdge(typename PFP::MAP& map, Dart d, const VertexAttribute<typename PFP::VEC3>& position)
168 169 170
{
	typedef typename PFP::VEC3 VEC3 ;

171
	if(map.isBoundaryEdge(d))
172 173
		return 0 ;

174
	Dart dd = map.phi2(d) ;
175 176 177 178 179 180 181
	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) ;
182

183
	// the following trick is useful for avoiding NaNs (due to floating point errors)
184
	if (c > 0.5) a = asin(s) ;
185 186
	else
	{
187 188 189
		if(c < -1) c = -1 ;
		if (s >= 0) a = acos(c) ;
		else a = -acos(c) ;
190
	}
Pierre Kraemer's avatar
Pierre Kraemer committed
191 192
//	if (isnan(a))
	if(a != a)
193
		std::cerr<< "Warning : computeAngleBetweenNormalsOnEdge returns NaN on edge " << d << "-" << dd << std::endl ;
194 195 196 197
	return a ;
}

template <typename PFP>
198
void computeAnglesBetweenNormalsOnEdges(typename PFP::MAP& map, const VertexAttribute<typename PFP::VEC3>& position, EdgeAttribute<typename PFP::REAL>& angles, const FunctorSelect& select, unsigned int thread)
199
{
Pierre Kraemer's avatar
Pierre Kraemer committed
200
	TraversorE<typename PFP::MAP> trav(map, select, thread);
201
	for (Dart d = trav.begin(); d != trav.end(); d = trav.next())
Pierre Kraemer's avatar
Pierre Kraemer committed
202
		angles[d] = computeAngleBetweenNormalsOnEdge<PFP>(map, d, position) ;
203 204
}

Pierre Kraemer's avatar
Pierre Kraemer committed
205 206 207 208 209
} // namespace Geometry

} // namespace Algo

} // namespace CGoGN