normal.hpp 6.79 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-2011, 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.u-strasbg.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 43 44
#include <cmath>

namespace CGoGN
{

namespace Algo
{

namespace Geometry
{

template <typename PFP>
typename PFP::VEC3 triangleNormal(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position)
{
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 51 52
template<typename PFP>
typename PFP::VEC3 newellNormal(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position)
{
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 70 71
template <typename PFP>
typename PFP::VEC3 faceNormal(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position)
{
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 93 94 95 96 97
}

template <typename PFP>
typename PFP::VEC3 vertexNormal(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position)
{
	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 117 118 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 146 147 148 149
template <typename PFP>
typename PFP::VEC3 vertexBorderNormal(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position)
{
	assert(map.dimension() == 3);

	typedef typename PFP::VEC3 VEC3 ;

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

	FunctorStore fs(faces);
	map.foreach_dart_of_oriented_vertex(d,fs);

	for(std::vector<Dart>::iterator it = faces.begin() ; it != faces.end() ; ++it)
	{
		if(!f.isMarked(*it) && map.phi3(*it)==*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 ;
}

Pierre Kraemer's avatar
Pierre Kraemer committed
150
template <typename PFP>
Sylvain Thery's avatar
Sylvain Thery committed
151
void computeNormalFaces(typename PFP::MAP& map, const typename PFP::TVEC3& position, typename PFP::TVEC3& 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>
Sylvain Thery's avatar
Sylvain Thery committed
159
void computeNormalVertices(typename PFP::MAP& map, const typename PFP::TVEC3& position, typename PFP::TVEC3& 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>
Pierre Kraemer's avatar
Pierre Kraemer committed
167
typename PFP::REAL computeAngleBetweenNormalsOnEdge(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& 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>
Pierre Kraemer's avatar
Pierre Kraemer committed
198
void computeAnglesBetweenNormalsOnEdges(typename PFP::MAP& map, const typename PFP::TVEC3& position, typename PFP::TREAL& 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