Création d'un compte pour un collaborateur extérieur au laboratoire depuis l'intranet ICube : https://intranet.icube.unistra.fr/fr/labs/member/profile

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                                        *
*                                                                              *
*******************************************************************************/

Sylvain Thery's avatar
Sylvain Thery committed
25
#include "Topology/generic/traversor/traversorCell.h"
26
#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