normal.hpp 6.48 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 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42
 * Contact information: cgogn@unistra.fr                                        *
 *                                                                              *
 *******************************************************************************/

#include "Algo/Geometry/basic.h"
#include "Algo/Geometry/area.h"
#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)
{
	typedef typename PFP::VEC3 VEC3 ;

Pierre Kraemer's avatar
Pierre Kraemer committed
43 44 45
	const VEC3& p1 = position[d];
	const VEC3& p2 = position[map.phi1(d)];
	const VEC3& p3 = position[map.phi_1(d)];
Pierre Kraemer's avatar
Pierre Kraemer committed
46 47 48 49 50 51 52

	VEC3 N = Geom::triangleNormal(p1, p2, p3) ;
	N.normalize() ;

	return N ;
}

53 54 55
template<typename PFP>
typename PFP::VEC3 newellNormal(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position)
{
56
	Dart e = d;
57 58 59 60 61 62
	typename PFP::VEC3 normal(0);
	do
	{
		const typename PFP::VEC3& P = position[e];
		e = map.phi1(e);
		const typename PFP::VEC3& Q = position[e];
63 64 65 66
		normal[0] += (P[1] - Q[1]) * (P[2] + Q[2]);
		normal[1] += (P[2] - Q[2]) * (P[0] + Q[0]);
		normal[2] += (P[0] - Q[0]) * (P[1] + Q[1]);
	} while (e != d);
67 68 69 70 71

	normal.normalize();
	return normal;
}

Pierre Kraemer's avatar
Pierre Kraemer committed
72 73 74 75 76 77 78 79 80 81 82 83 84 85
template <typename PFP>
typename PFP::VEC3 faceNormal(typename PFP::MAP& map, Dart d, const typename PFP::TVEC3& position)
{
	typedef typename PFP::VEC3 VEC3 ;

	if(map.isFaceTriangle(d))
		return triangleNormal<PFP>(map, d, position) ;
	else
	{
		VEC3 N(0) ;
		Dart it = d ;
		do
		{
			VEC3 n = triangleNormal<PFP>(map, it, position) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
86 87
			//if(!std::isnan(n[0]) && !std::isnan(n[1]) && !std::isnan(n[2]))
			if (n[0] == n[0] && n[1] == n[1] && n[2] == n[2])
Pierre Kraemer's avatar
Pierre Kraemer committed
88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105
				N += n ;
			it = map.phi1(it) ;
		} while (it != d) ;
		N.normalize() ;
		return N ;
	}
}

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) ;
	Dart it = d ;
	do
	{
		VEC3 n = faceNormal<PFP>(map, it, position) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
106 107 108 109 110 111 112
		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 ;
		}
Pierre Kraemer's avatar
Pierre Kraemer committed
113 114 115 116 117 118 119
		it = map.phi1(map.phi2(it)) ;
	} while (it != d) ;
	N.normalize() ;
	return N ;
}

template <typename PFP>
Sylvain Thery's avatar
Sylvain Thery committed
120
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
121
{
122
	CellMarker marker(map, FACE,thread);
Pierre Kraemer's avatar
Pierre Kraemer committed
123 124 125 126 127 128 129 130 131 132
	for(Dart d = map.begin(); d != map.end(); map.next(d))
	{
		if(select(d) && !marker.isMarked(d))
		{
			marker.mark(d);
			face_normal[d] = faceNormal<PFP>(map, d, position) ;
		}
	}
}

Sylvain Thery's avatar
Sylvain Thery committed
133
template <typename PFP>
134
class computeNormalVerticesFunctor : public FunctorMap<typename PFP::MAP>
Sylvain Thery's avatar
Sylvain Thery committed
135 136 137 138 139 140 141
{
protected:
	typename PFP::MAP& m_map;
	const typename PFP::TVEC3& m_position;
	typename PFP::TVEC3& m_normal;
public:
	computeNormalVerticesFunctor(typename PFP::MAP& map, const typename PFP::TVEC3& position, typename PFP::TVEC3& normal):
142 143
		m_map(map), m_position(position), m_normal(normal)
	{}
Sylvain Thery's avatar
Sylvain Thery committed
144 145 146 147 148 149 150
	bool operator()(Dart d)
	{
		m_normal[d] = vertexNormal<PFP>(m_map, d, m_position) ;
		return false;
	}
};

Pierre Kraemer's avatar
Pierre Kraemer committed
151
template <typename PFP>
Sylvain Thery's avatar
Sylvain Thery committed
152
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
153
{
154
	CellMarker marker(map, VERTEX, thread);
Pierre Kraemer's avatar
Pierre Kraemer committed
155 156 157 158 159 160 161 162 163 164
	for(Dart d = map.begin(); d != map.end(); map.next(d))
	{
		if(select(d) && !marker.isMarked(d))
		{
			marker.mark(d);
			normal[d] = vertexNormal<PFP>(map, d, position) ;
		}
	}
}

165 166 167 168 169 170 171 172 173 174 175 176 177 178
template <typename PFP>
typename PFP::REAL computeAngleBetweenNormalsOnEdge(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position)
{
	typedef typename PFP::VEC3 VEC3 ;

	Dart dd = map.phi2(d) ;
	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) ;
	// the following trick is useful for avoiding NaNs (due to floating point errors)
179
	if (c > 0.5) a = asin(s) ;
180 181
	else
	{
182 183 184
		if(c < -1) c = -1 ;
		if (s >= 0) a = acos(c) ;
		else a = -acos(c) ;
185 186
	}
	if (isnan(a))
187
		std::cerr<< "Warning : computeAngleBetweenNormalsOnEdge returns NaN on edge " << d << "-" << dd << std::endl ;
188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204
	return a ;
}

template <typename PFP>
void computeAnglesBetweenNormalsOnEdges(typename PFP::MAP& map, typename PFP::TVEC3& position, typename PFP::TREAL& angles, const FunctorSelect& select, unsigned int thread)
{
	CellMarker me(map, EDGE, thread) ;
	for(Dart d = map.begin(); d != map.end(); map.next(d))
	{
		if(select(d) && !me.isMarked(d))
		{
			me.mark(d) ;
			angles[d] = computeAngleBetweenNormalsOnEdge<PFP>(map, d, position) ;
		}
	}
}

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

} // namespace Algo

} // namespace CGoGN