average_normals.h 14.8 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
* Contact information: cgogn@unistra.fr                                        *
*                                                                              *
*******************************************************************************/

25 26
#include "Topology/generic/traversorCell.h"
#include "Topology/generic/traversor2.h"
Pierre Kraemer's avatar
Pierre Kraemer committed
27 28 29 30 31 32 33

namespace CGoGN
{

namespace Algo
{

34 35 36
namespace Surface
{

Pierre Kraemer's avatar
Pierre Kraemer committed
37
namespace Filtering
Pierre Kraemer's avatar
Pierre Kraemer committed
38 39 40 41 42 43 44 45 46
{

/**
 * compute new position of vertices from normals (normalAverage & MMSE filters)
 * @param map the map
 */
template <typename PFP>
void computeNewPositionsFromFaceNormals(
	typename PFP::MAP& map,
Pierre Kraemer's avatar
Pierre Kraemer committed
47 48 49 50 51 52
	const VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& position,
	VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& position2,
	const FaceAttribute<typename PFP::REAL, typename PFP::MAP::IMPL>& faceArea,
	const FaceAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& faceCentroid,
	const FaceAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& faceNormal,
	const FaceAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& faceNewNormal)
Pierre Kraemer's avatar
Pierre Kraemer committed
53 54 55 56
{
	typedef typename PFP::VEC3 VEC3 ;
	typedef typename PFP::REAL REAL ;

57
	TraversorV<typename PFP::MAP> t(map) ;
58
	for(Dart d = t.begin(); d != t.end(); d = t.next())
Pierre Kraemer's avatar
Pierre Kraemer committed
59
	{
Pierre Kraemer's avatar
Pierre Kraemer committed
60
		const VEC3& pos_d = position[d] ;
Pierre Kraemer's avatar
Pierre Kraemer committed
61

Pierre Kraemer's avatar
Pierre Kraemer committed
62 63
		VEC3 displ(0) ;
		REAL sumAreas = 0 ;
64

Pierre Kraemer's avatar
Pierre Kraemer committed
65 66 67 68 69 70 71
		Traversor2VF<typename PFP::MAP> tvf(map, d) ;
		for(Dart it = tvf.begin(); it != tvf.end(); it = tvf.next())
		{
			sumAreas += faceArea[it] ;
			VEC3 vT = faceCentroid[it] - pos_d ;
			vT = (vT * faceNewNormal[it]) * faceNormal[it] ;
			displ += faceArea[it] * vT ;
Pierre Kraemer's avatar
Pierre Kraemer committed
72
		}
Pierre Kraemer's avatar
Pierre Kraemer committed
73 74 75

		displ /= sumAreas ;
		position2[d] = pos_d + displ ;
Pierre Kraemer's avatar
Pierre Kraemer committed
76 77 78 79
	}
}

template <typename PFP>
Pierre Kraemer's avatar
Pierre Kraemer committed
80
void filterAverageNormals(typename PFP::MAP& map, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& position, VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& position2)
Pierre Kraemer's avatar
Pierre Kraemer committed
81
{
Pierre Kraemer's avatar
Pierre Kraemer committed
82
	typedef typename PFP::MAP::IMPL MAP_IMPL ;
Pierre Kraemer's avatar
Pierre Kraemer committed
83 84 85
	typedef typename PFP::VEC3 VEC3 ;
	typedef typename PFP::REAL REAL ;

Pierre Kraemer's avatar
Pierre Kraemer committed
86 87 88
	FaceAutoAttribute<REAL, MAP_IMPL> faceArea(map, "faceArea") ;
	FaceAutoAttribute<VEC3, MAP_IMPL> faceNormal(map, "faceNormal") ;
	FaceAutoAttribute<VEC3, MAP_IMPL> faceCentroid(map, "faceCentroid") ;
Pierre Kraemer's avatar
Pierre Kraemer committed
89

90 91 92
	Algo::Surface::Geometry::computeAreaFaces<PFP>(map, position, faceArea) ;
	Algo::Surface::Geometry::computeNormalFaces<PFP>(map, position, faceNormal) ;
	Algo::Surface::Geometry::computeCentroidFaces<PFP>(map, position, faceCentroid) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
93

Pierre Kraemer's avatar
Pierre Kraemer committed
94
	FaceAutoAttribute<VEC3, MAP_IMPL> faceNewNormal(map, "faceNewNormal") ;
Pierre Kraemer's avatar
Pierre Kraemer committed
95 96

	// Compute new normals
97
	TraversorF<typename PFP::MAP> tf(map) ;
98
	for(Dart d = tf.begin(); d != tf.end(); d = tf.next())
Pierre Kraemer's avatar
Pierre Kraemer committed
99
	{
Pierre Kraemer's avatar
Pierre Kraemer committed
100 101
		REAL sumArea = 0 ;
		VEC3 meanFilter(0) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
102

Pierre Kraemer's avatar
Pierre Kraemer committed
103 104 105 106 107 108
		// traversal of adjacent faces (by edges and vertices)
		Traversor2FFaV<typename PFP::MAP> taf(map, d) ;
		for(Dart it = taf.begin(); it != taf.end(); it = taf.next())
		{
			sumArea += faceArea[it] ;
			meanFilter += faceArea[it] * faceNormal[it] ;
Pierre Kraemer's avatar
Pierre Kraemer committed
109
		}
Pierre Kraemer's avatar
Pierre Kraemer committed
110 111 112 113 114 115

		// finalize the computation of meanFilter normal
		meanFilter /= sumArea ;
		meanFilter.normalize() ;
		// and store it
		faceNewNormal[d] = meanFilter ;
Pierre Kraemer's avatar
Pierre Kraemer committed
116 117 118 119
	}

	// Compute new vertices position
	computeNewPositionsFromFaceNormals<PFP>(
120
		map, position, position2, faceArea, faceCentroid, faceNormal, faceNewNormal) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
121 122 123
}

template <typename PFP>
Pierre Kraemer's avatar
Pierre Kraemer committed
124
void filterMMSE(typename PFP::MAP& map, float sigmaN2, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& position, VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& position2)
Pierre Kraemer's avatar
Pierre Kraemer committed
125
{
Pierre Kraemer's avatar
Pierre Kraemer committed
126
	typedef typename PFP::MAP::IMPL MAP_IMPL ;
Pierre Kraemer's avatar
Pierre Kraemer committed
127 128 129
	typedef typename PFP::VEC3 VEC3 ;
	typedef typename PFP::REAL REAL ;

Pierre Kraemer's avatar
Pierre Kraemer committed
130 131 132
	FaceAutoAttribute<REAL, MAP_IMPL> faceArea(map, "faceArea") ;
	FaceAutoAttribute<VEC3, MAP_IMPL> faceNormal(map, "faceNormal") ;
	FaceAutoAttribute<VEC3, MAP_IMPL> faceCentroid(map, "faceCentroid") ;
Pierre Kraemer's avatar
Pierre Kraemer committed
133

134 135 136
	Algo::Surface::Geometry::computeAreaFaces<PFP>(map, position, faceArea) ;
	Algo::Surface::Geometry::computeNormalFaces<PFP>(map, position, faceNormal) ;
	Algo::Surface::Geometry::computeCentroidFaces<PFP>(map, position, faceCentroid) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
137

Pierre Kraemer's avatar
Pierre Kraemer committed
138
	FaceAutoAttribute<VEC3, MAP_IMPL> faceNewNormal(map, "faceNewNormal") ;
Pierre Kraemer's avatar
Pierre Kraemer committed
139 140

	// Compute new normals
141
	TraversorF<typename PFP::MAP> tf(map) ;
142
	for(Dart d = tf.begin(); d != tf.end(); d = tf.next())
Pierre Kraemer's avatar
Pierre Kraemer committed
143
	{
Pierre Kraemer's avatar
Pierre Kraemer committed
144 145 146 147 148
		// traversal of neighbour vertices
		REAL sumArea = 0 ;
		REAL sigmaX2 = 0 ;
		REAL sigmaY2 = 0 ;
		REAL sigmaZ2 = 0 ;
Pierre Kraemer's avatar
Pierre Kraemer committed
149

Pierre Kraemer's avatar
Pierre Kraemer committed
150
		VEC3 meanFilter(0) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
151

Pierre Kraemer's avatar
Pierre Kraemer committed
152 153 154 155 156 157 158 159 160 161 162 163 164
		// traversal of adjacent faces (by edges and vertices)
		Traversor2FFaV<typename PFP::MAP> taf(map, d) ;
		for(Dart it = taf.begin(); it != taf.end(); it = taf.next())
		{
			// get info from face embedding and sum
			REAL area = faceArea[it] ;
			sumArea += area ;
			VEC3 normal = faceNormal[it] ;
			meanFilter += area * normal ;
			sigmaX2 += area * normal[0] * normal[0] ;
			sigmaY2 += area * normal[1] * normal[1] ;
			sigmaZ2 += area * normal[2] * normal[2] ;
		}
Pierre Kraemer's avatar
Pierre Kraemer committed
165

Pierre Kraemer's avatar
Pierre Kraemer committed
166 167 168 169 170 171 172
		meanFilter /= sumArea ;
		sigmaX2 /= sumArea ;
		sigmaX2 -= meanFilter[0] * meanFilter[0] ;
		sigmaY2 /= sumArea ;
		sigmaY2 -= meanFilter[1] * meanFilter[1] ;
		sigmaZ2 /= sumArea ;
		sigmaZ2 -= meanFilter[2] * meanFilter[2] ;
Pierre Kraemer's avatar
Pierre Kraemer committed
173

Pierre Kraemer's avatar
Pierre Kraemer committed
174 175
		VEC3& oldNormal = faceNormal[d] ;
		VEC3 newNormal ;
Pierre Kraemer's avatar
Pierre Kraemer committed
176

Pierre Kraemer's avatar
Pierre Kraemer committed
177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196
		if(sigmaX2 < sigmaN2)
			newNormal[0] = meanFilter[0] ;
		else
		{
			newNormal[0] = (1 - (sigmaN2 / sigmaX2)) * oldNormal[0] ;
			newNormal[0] += (sigmaN2 / sigmaX2) * meanFilter[0] ;
		}
		if(sigmaY2 < sigmaN2)
			newNormal[1] = meanFilter[1] ;
		else
		{
			newNormal[1] = (1 - (sigmaN2 / sigmaY2)) * oldNormal[1] ;
			newNormal[1] += (sigmaN2 / sigmaY2) * meanFilter[1] ;
		}
		if(sigmaZ2 < sigmaN2)
			newNormal[2] = meanFilter[2] ;
		else
		{
			newNormal[2] = (1 - (sigmaN2 / sigmaZ2)) * oldNormal[2] ;
			newNormal[2] += (sigmaN2 / sigmaZ2) * meanFilter[2] ;
Pierre Kraemer's avatar
Pierre Kraemer committed
197
		}
Pierre Kraemer's avatar
Pierre Kraemer committed
198 199 200

		newNormal.normalize() ;
		faceNewNormal[d] = newNormal ;
Pierre Kraemer's avatar
Pierre Kraemer committed
201 202 203 204
	}

	// Compute new vertices position
	computeNewPositionsFromFaceNormals<PFP>(
205
		map, position, position2, faceArea, faceCentroid, faceNormal, faceNewNormal) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
206 207 208
}

template <typename PFP>
Pierre Kraemer's avatar
Pierre Kraemer committed
209
void filterTNBA(typename PFP::MAP& map, float sigmaN2, float SUSANthreshold, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& position, VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& position2)
Pierre Kraemer's avatar
Pierre Kraemer committed
210
{
Pierre Kraemer's avatar
Pierre Kraemer committed
211
	typedef typename PFP::MAP::IMPL MAP_IMPL ;
Pierre Kraemer's avatar
Pierre Kraemer committed
212 213 214
	typedef typename PFP::VEC3 VEC3 ;
	typedef typename PFP::REAL REAL ;

Pierre Kraemer's avatar
Pierre Kraemer committed
215 216 217
	FaceAutoAttribute<REAL, MAP_IMPL> faceArea(map, "faceArea") ;
	FaceAutoAttribute<VEC3, MAP_IMPL> faceNormal(map, "faceNormal") ;
	FaceAutoAttribute<VEC3, MAP_IMPL> faceCentroid(map, "faceCentroid") ;
Pierre Kraemer's avatar
Pierre Kraemer committed
218

219 220 221
	Algo::Surface::Geometry::computeAreaFaces<PFP>(map, position, faceArea) ;
	Algo::Surface::Geometry::computeNormalFaces<PFP>(map, position, faceNormal) ;
	Algo::Surface::Geometry::computeCentroidFaces<PFP>(map, position, faceCentroid) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
222

Pierre Kraemer's avatar
Pierre Kraemer committed
223
	FaceAutoAttribute<VEC3, MAP_IMPL> faceNewNormal(map, "faceNewNormal") ;
Pierre Kraemer's avatar
Pierre Kraemer committed
224 225 226 227 228 229

	// Compute new normals
	long nbTot = 0 ;
	long nbAdapt = 0 ;
	long nbSusan = 0 ;

230
	TraversorF<typename PFP::MAP> tf(map) ;
231
	for(Dart d = tf.begin(); d != tf.end(); d = tf.next())
Pierre Kraemer's avatar
Pierre Kraemer committed
232
	{
Pierre Kraemer's avatar
Pierre Kraemer committed
233
		const VEC3& normF = faceNormal[d] ;
Pierre Kraemer's avatar
Pierre Kraemer committed
234

Pierre Kraemer's avatar
Pierre Kraemer committed
235 236 237 238 239
		// traversal of neighbour vertices
		REAL sumArea = 0 ;
		REAL sigmaX2 = 0 ;
		REAL sigmaY2 = 0 ;
		REAL sigmaZ2 = 0 ;
Pierre Kraemer's avatar
Pierre Kraemer committed
240

Pierre Kraemer's avatar
Pierre Kraemer committed
241 242 243 244 245 246 247 248 249
		VEC3 meanFilter(0) ;
		bool SUSANregion = false ;

		// traversal of adjacent faces (by edges and vertices)
		Traversor2FFaV<typename PFP::MAP> taf(map, d) ;
		for(Dart it = taf.begin(); it != taf.end(); it = taf.next())
		{
			// get info from face embedding and sum
			const VEC3& normal = faceNormal[it] ;
Pierre Kraemer's avatar
Pierre Kraemer committed
250

Pierre Kraemer's avatar
Pierre Kraemer committed
251 252
			float angle = Geom::angle(normF, normal) ;
			if(angle <= SUSANthreshold)
Pierre Kraemer's avatar
Pierre Kraemer committed
253
			{
Pierre Kraemer's avatar
Pierre Kraemer committed
254 255 256 257 258 259
				REAL area = faceArea[it] ;
				sumArea += area ;
				meanFilter += area * normal ;
				sigmaX2 += area * normal[0] * normal[0] ;
				sigmaY2 += area * normal[1] * normal[1] ;
				sigmaZ2 += area * normal[2] * normal[2] ;
260
			}
Pierre Kraemer's avatar
Pierre Kraemer committed
261 262
			else SUSANregion = true ;
		}
Pierre Kraemer's avatar
Pierre Kraemer committed
263

Pierre Kraemer's avatar
Pierre Kraemer committed
264 265
		if(SUSANregion)
			++nbSusan ;
Pierre Kraemer's avatar
Pierre Kraemer committed
266

Pierre Kraemer's avatar
Pierre Kraemer committed
267
		++nbTot ;
Pierre Kraemer's avatar
Pierre Kraemer committed
268

Pierre Kraemer's avatar
Pierre Kraemer committed
269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285
		if(sumArea > 0.0f)
		{
			meanFilter /= sumArea ;
			sigmaX2 /= sumArea ;
			sigmaX2 -= meanFilter[0] * meanFilter[0] ;
			sigmaY2 /= sumArea ;
			sigmaY2 -= meanFilter[1] * meanFilter[1] ;
			sigmaZ2 /= sumArea ;
			sigmaZ2 -= meanFilter[2] * meanFilter[2] ;

			VEC3& oldNormal = faceNormal[d] ;
			VEC3 newNormal ;

			bool adapt = false ;
			if(sigmaX2 < sigmaN2)
				newNormal[0] = meanFilter[0] ;
			else
Pierre Kraemer's avatar
Pierre Kraemer committed
286
			{
Pierre Kraemer's avatar
Pierre Kraemer committed
287 288 289
				adapt = true ;
				newNormal[0] = (1 - (sigmaN2 / sigmaX2)) * oldNormal[0] ;
				newNormal[0] += (sigmaN2 / sigmaX2) * meanFilter[0] ;
Pierre Kraemer's avatar
Pierre Kraemer committed
290
			}
Pierre Kraemer's avatar
Pierre Kraemer committed
291 292
			if(sigmaY2 < sigmaN2)
				newNormal[1] = meanFilter[1] ;
Pierre Kraemer's avatar
Pierre Kraemer committed
293 294
			else
			{
Pierre Kraemer's avatar
Pierre Kraemer committed
295 296 297
				adapt = true ;
				newNormal[1] = (1 - (sigmaN2 / sigmaY2)) * oldNormal[1] ;
				newNormal[1] += (sigmaN2 / sigmaY2) * meanFilter[1] ;
Pierre Kraemer's avatar
Pierre Kraemer committed
298
			}
Pierre Kraemer's avatar
Pierre Kraemer committed
299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315
			if(sigmaZ2 < sigmaN2)
				newNormal[2] = meanFilter[2] ;
			else
			{
				adapt = true ;
				newNormal[2] = (1 - (sigmaN2 / sigmaZ2)) * oldNormal[2] ;
				newNormal[2] += (sigmaN2 / sigmaZ2) * meanFilter[2] ;
			}
			if(adapt)
				++nbAdapt ;

			newNormal.normalize() ;
			faceNewNormal[d] = newNormal;
		}
		else
		{
			faceNewNormal[d] = normF ;
Pierre Kraemer's avatar
Pierre Kraemer committed
316 317 318 319 320
		}
	}

	// Compute new vertices position
	computeNewPositionsFromFaceNormals<PFP>(
321
		map, position, position2, faceArea, faceCentroid, faceNormal, faceNewNormal) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
322

323 324
//	CGoGNout <<" susan rate = "<< float(nbSusan)/float(nbTot)<<CGoGNendl;
//	CGoGNout <<" adaptive rate = "<< float(nbAdapt)/float(nbTot)<<CGoGNendl;
Pierre Kraemer's avatar
Pierre Kraemer committed
325 326 327
}

template <typename PFP>
Pierre Kraemer's avatar
Pierre Kraemer committed
328
void filterVNBA(typename PFP::MAP& map, float sigmaN2, float SUSANthreshold, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& position, VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& position2, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP::IMPL>& normal)
Pierre Kraemer's avatar
Pierre Kraemer committed
329
{
Pierre Kraemer's avatar
Pierre Kraemer committed
330
	typedef typename PFP::MAP::IMPL MAP_IMPL ;
Pierre Kraemer's avatar
Pierre Kraemer committed
331 332 333
	typedef typename PFP::VEC3 VEC3 ;
	typedef typename PFP::REAL REAL ;

Pierre Kraemer's avatar
Pierre Kraemer committed
334 335 336
	FaceAutoAttribute<REAL, MAP_IMPL> faceArea(map, "faceArea") ;
	FaceAutoAttribute<VEC3, MAP_IMPL> faceNormal(map, "faceNormal") ;
	FaceAutoAttribute<VEC3, MAP_IMPL> faceCentroid(map, "faceCentroid") ;
Pierre Kraemer's avatar
Pierre Kraemer committed
337

338 339 340
	Algo::Surface::Geometry::computeAreaFaces<PFP>(map, position, faceArea) ;
	Algo::Surface::Geometry::computeNormalFaces<PFP>(map, position, faceNormal) ;
	Algo::Surface::Geometry::computeCentroidFaces<PFP>(map, position, faceCentroid) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
341

Pierre Kraemer's avatar
Pierre Kraemer committed
342 343 344
	VertexAutoAttribute<REAL, MAP_IMPL> vertexArea(map, "vertexArea") ;
	FaceAutoAttribute<VEC3, MAP_IMPL> faceNewNormal(map, "faceNewNormal") ;
	VertexAutoAttribute<VEC3, MAP_IMPL> vertexNewNormal(map, "vertexNewNormal") ;
Pierre Kraemer's avatar
Pierre Kraemer committed
345 346 347 348 349

	long nbTot = 0 ;
	long nbAdapt = 0 ;
	long nbSusan = 0 ;

350
	TraversorV<typename PFP::MAP> tv(map) ;
351
	for(Dart d = tv.begin(); d != tv.end(); d = tv.next())
Pierre Kraemer's avatar
Pierre Kraemer committed
352
	{
Pierre Kraemer's avatar
Pierre Kraemer committed
353
		const VEC3& normV = normal[d] ;
Pierre Kraemer's avatar
Pierre Kraemer committed
354

Pierre Kraemer's avatar
Pierre Kraemer committed
355 356 357 358
		REAL sumArea = 0 ;
		REAL sigmaX2 = 0 ;
		REAL sigmaY2 = 0 ;
		REAL sigmaZ2 = 0 ;
Pierre Kraemer's avatar
Pierre Kraemer committed
359

Pierre Kraemer's avatar
Pierre Kraemer committed
360
		VEC3 meanFilter(0) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
361

Pierre Kraemer's avatar
Pierre Kraemer committed
362
		bool SUSANregion = false ;
Pierre Kraemer's avatar
Pierre Kraemer committed
363

Pierre Kraemer's avatar
Pierre Kraemer committed
364 365 366 367 368 369 370
		// traversal of neighbour vertices
		Traversor2VVaE<typename PFP::MAP> tav(map, d) ;
		for(Dart it = tav.begin(); it != tav.end(); it = tav.next())
		{
			const VEC3& neighborNormal = normal[it] ;
			float angle = Geom::angle(normV, neighborNormal) ;
			if( angle <= SUSANthreshold )
Pierre Kraemer's avatar
Pierre Kraemer committed
371
			{
Sylvain Thery's avatar
Sylvain Thery committed
372
				REAL umbArea = Algo::Surface::Geometry::vertexOneRingArea<PFP>(map, it, position) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
373 374 375 376 377 378 379
				vertexArea[it] = umbArea ;

				sumArea += umbArea ;
				sigmaX2 += umbArea * neighborNormal[0] * neighborNormal[0] ;
				sigmaY2 += umbArea * neighborNormal[1] * neighborNormal[1] ;
				sigmaZ2 += umbArea * neighborNormal[2] * neighborNormal[2] ;
				meanFilter += neighborNormal * umbArea ;
380
			}
Pierre Kraemer's avatar
Pierre Kraemer committed
381 382
			else SUSANregion = true ;
		}
Pierre Kraemer's avatar
Pierre Kraemer committed
383

Pierre Kraemer's avatar
Pierre Kraemer committed
384 385
		if(SUSANregion)
			++nbSusan ;
Pierre Kraemer's avatar
Pierre Kraemer committed
386

Pierre Kraemer's avatar
Pierre Kraemer committed
387
		++nbTot ;
Pierre Kraemer's avatar
Pierre Kraemer committed
388

Pierre Kraemer's avatar
Pierre Kraemer committed
389 390 391 392 393 394 395 396 397 398 399 400 401 402 403
		if(sumArea > 0.0f)
		{
			meanFilter /= sumArea ;
			sigmaX2 /= sumArea ;
			sigmaX2 -= meanFilter[0] * meanFilter[0] ;
			sigmaY2 /= sumArea ;
			sigmaY2 -= meanFilter[1] * meanFilter[1] ;
			sigmaZ2 /= sumArea ;
			sigmaZ2 -= meanFilter[2] * meanFilter[2] ;

			VEC3 newNormal ;
			bool adapt = false ;
			if(sigmaX2 < sigmaN2)
				newNormal[0] = meanFilter[0] ;
			else
Pierre Kraemer's avatar
Pierre Kraemer committed
404
			{
Pierre Kraemer's avatar
Pierre Kraemer committed
405 406 407
				adapt = true ;
				newNormal[0] = (1 - (sigmaN2 / sigmaX2)) * normV[0] ;
				newNormal[0] += (sigmaN2 / sigmaX2) * meanFilter[0] ;
Pierre Kraemer's avatar
Pierre Kraemer committed
408
			}
Pierre Kraemer's avatar
Pierre Kraemer committed
409 410
			if(sigmaY2 < sigmaN2)
				newNormal[1] = meanFilter[1] ;
Pierre Kraemer's avatar
Pierre Kraemer committed
411 412
			else
			{
Pierre Kraemer's avatar
Pierre Kraemer committed
413 414 415
				adapt = true ;
				newNormal[1] = (1 - (sigmaN2 / sigmaY2)) * normV[1] ;
				newNormal[1] += (sigmaN2 / sigmaY2) * meanFilter[1] ;
Pierre Kraemer's avatar
Pierre Kraemer committed
416
			}
Pierre Kraemer's avatar
Pierre Kraemer committed
417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434
			if(sigmaZ2 < sigmaN2)
				newNormal[2] = meanFilter[2] ;
			else
			{
				adapt = true ;
				newNormal[2] = (1 - (sigmaN2 / sigmaZ2)) * normV[2] ;
				newNormal[2] += (sigmaN2 / sigmaZ2) * meanFilter[2] ;
			}

			if(adapt)
				++nbAdapt ;

			newNormal.normalize() ;
			vertexNewNormal[d] = newNormal ;
		}
		else
		{
			vertexNewNormal[d] = normV ;
Pierre Kraemer's avatar
Pierre Kraemer committed
435 436 437 438
		}
	}

	// Compute face normals from vertex normals
439
	TraversorF<typename PFP::MAP> tf(map) ;
440
	for(Dart d = tf.begin(); d != tf.end(); d = tf.next())
Pierre Kraemer's avatar
Pierre Kraemer committed
441
	{
Pierre Kraemer's avatar
Pierre Kraemer committed
442 443 444 445
		VEC3 newNormal(0) ;
		REAL totArea = 0 ;
		Traversor2FV<typename PFP::MAP> tav(map, d) ;
		for(Dart it = tav.begin(); it != tav.end(); it = tav.next())
Pierre Kraemer's avatar
Pierre Kraemer committed
446
		{
Pierre Kraemer's avatar
Pierre Kraemer committed
447 448 449 450 451
			VEC3 vNorm = vertexNewNormal[it] ;
			REAL area = vertexArea[it] ;
			totArea += area ;
			vNorm *= area ;
			newNormal += vNorm ;
Pierre Kraemer's avatar
Pierre Kraemer committed
452
		}
Pierre Kraemer's avatar
Pierre Kraemer committed
453 454
		newNormal.normalize() ;
		faceNewNormal[d] = newNormal ;
Pierre Kraemer's avatar
Pierre Kraemer committed
455 456 457 458
	}

	// Compute new vertices position
	computeNewPositionsFromFaceNormals<PFP>(
459
		map, position, position2, faceArea, faceCentroid, faceNormal, faceNewNormal	) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
460

461 462
//	CGoGNout <<" susan rate = "<< float(nbSusan)/float(nbTot)<<CGoGNendl;
//	CGoGNout <<" adaptive rate = "<< float(nbAdapt)/float(nbTot)<<CGoGNendl;
Pierre Kraemer's avatar
Pierre Kraemer committed
463 464
}

Pierre Kraemer's avatar
Pierre Kraemer committed
465
} // namespace Filtering
Pierre Kraemer's avatar
Pierre Kraemer committed
466

Pierre Kraemer's avatar
Pierre Kraemer committed
467
} // namespace Surface
468

Pierre Kraemer's avatar
Pierre Kraemer committed
469
} // namespace Algo
Pierre Kraemer's avatar
Pierre Kraemer committed
470

Pierre Kraemer's avatar
Pierre Kraemer committed
471
} // namespace CGoGN