filters_Primal.h 29.7 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29
/*******************************************************************************
* CGoGN: Combinatorial and Geometric modeling with Generic N-dimensional Maps  *
* version 0.1                                                                  *
* Copyright (C) 2009-2012, IGG Team, LSIIT, University of Strasbourg           *
*                                                                              *
* 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.           *
*                                                                              *
* Web site: http://cgogn.unistra.fr/                                           *
* Contact information: cgogn@unistra.fr                                        *
*                                                                              *
*******************************************************************************/

#ifndef __3MR_FILTERS_PRIMAL__
#define __3MR_FILTERS_PRIMAL__

#include <cmath>
#include "Algo/Geometry/centroid.h"
30
#include "Algo/Modelisation/tetrahedralization.h"
31 32 33 34

namespace CGoGN
{

35 36 37
namespace Algo
{

38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53
namespace Multiresolution
{

class MRFilter
{
public:
	MRFilter() {}
	virtual ~MRFilter() {}
	virtual void operator() () = 0 ;
} ;


/*********************************************************************************
 *                           LOOP BASIC FUNCTIONS
 *********************************************************************************/
template <typename PFP>
54
typename PFP::VEC3 loopOddVertex(typename PFP::MAP& map, const VertexAttribute<typename PFP::VEC3>& position, Dart d1)
55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73
{
	Dart d2 = map.phi2(d1) ;
	Dart d3 = map.phi_1(d1) ;
	Dart d4 = map.phi_1(d2) ;

	typename PFP::VEC3 p1 = position[d1] ;
	typename PFP::VEC3 p2 = position[d2] ;
	typename PFP::VEC3 p3 = position[d3] ;
	typename PFP::VEC3 p4 = position[d4] ;

	p1 *= 3.0 / 8.0 ;
	p2 *= 3.0 / 8.0 ;
	p3 *= 1.0 / 8.0 ;
	p4 *= 1.0 / 8.0 ;

	return p1 + p2 + p3 + p4 ;
}

template <typename PFP>
74
typename PFP::VEC3 loopEvenVertex(typename PFP::MAP& map, const VertexAttribute<typename PFP::VEC3>& position, Dart d)
75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99
{
	map.incCurrentLevel() ;

	typename PFP::VEC3 np(0) ;
	unsigned int degree = 0 ;
	Traversor2VVaE<typename PFP::MAP> trav(map, d) ;
	for(Dart it = trav.begin(); it != trav.end(); it = trav.next())
	{
		++degree ;
		np += position[it] ;
	}

	map.decCurrentLevel() ;

	float mu = 3.0/8.0 + 1.0/4.0 * cos(2.0 * M_PI / degree) ;
	mu = (5.0/8.0 - (mu * mu)) / degree ;
	np *= 8.0/5.0 * mu ;

	return np ;
}

/*********************************************************************************
 *          SHW04 BASIC FUNCTIONS : tetrahedral/octahedral meshes
 *********************************************************************************/
template <typename PFP>
100
typename PFP::VEC3 SHW04Vertex(typename PFP::MAP& map, const VertexAttribute<typename PFP::VEC3>& position, Dart d)
101 102 103
{
	typename PFP::VEC3 res(0);

104
	if(Algo::Modelisation::Tetrahedralization::isTetrahedron<PFP>(map, d))
105 106 107 108 109 110 111 112 113 114 115 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 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168
	{
		Dart d1 = map.phi1(d) ;
		Dart d2 = map.phi_1(d);
		Dart d3 = map.phi_1(map.phi2(d));

		typename PFP::VEC3 p = position[d];
		typename PFP::VEC3 p1 = position[d1] ;
		typename PFP::VEC3 p2 = position[d2] ;
		typename PFP::VEC3 p3 = position[d3] ;

		p *= -1;
		p1 *= 17.0 / 3.0;
		p2 *= 17.0 / 3.0;
		p3 *= 17.0 / 3.0;

		res += p + p1 + p2 + p3;
		res *= 1.0 / 16.0;
	}
	else
	{
		Dart d1 = map.phi1(d);
		Dart d2 = map.phi_1(d);
		Dart d3 = map.phi_1(map.phi2(d));
		Dart d4 = map.phi_1(map.phi2(d3));
		Dart d5 = map.phi_1(map.phi2(map.phi_1(d)));

		typename PFP::VEC3 p = position[d];
		typename PFP::VEC3 p1 = position[d1] ;
		typename PFP::VEC3 p2 = position[d2] ;
		typename PFP::VEC3 p3 = position[d3] ;
		typename PFP::VEC3 p4 = position[d4] ;
		typename PFP::VEC3 p5 = position[d5] ;

		p *= 3.0 / 4.0;
		p1 *= 1.0 / 6.0;
		p2 *= 1.0 / 6.0;
		p3 *= 1.0 / 6.0;
		p4 *= 7.0 / 12.0;
		p5 *= 1.0 / 6.0;

		res += p + p1 + p2 + p3 + p4 + p5;
		res *= 1.0 / 2.0;
	}

	return res;
}

/*********************************************************************************
 *          BSXW02 BASIC FUNCTIONS : polyhedral meshes
 *********************************************************************************/



/*********************************************************************************
 *                           ANALYSIS FILTERS
 *********************************************************************************/

/* Loop
 *********************************************************************************/
template <typename PFP>
class LoopEvenAnalysisFilter : public MRFilter
{
protected:
	typename PFP::MAP& m_map ;
169
	VertexAttribute<typename PFP::VEC3>& m_position ;
170 171

public:
172
	LoopEvenAnalysisFilter(typename PFP::MAP& m, VertexAttribute<typename PFP::VEC3>& p) : m_map(m), m_position(p)
173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194
	{}

	void operator() ()
	{
		TraversorV<typename PFP::MAP> trav(m_map) ;
		for (Dart d = trav.begin(); d != trav.end(); d = trav.next())
		{
			if(m_map.isBoundaryVertex(d))
			{
				Dart db = m_map.findBoundaryFaceOfVertex(d);
				typename PFP::VEC3 p = loopEvenVertex<PFP>(m_map, m_position, db) ;
				m_position[db] -= p ;
			}
		}
	}
} ;

template <typename PFP>
class LoopNormalisationAnalysisFilter : public MRFilter
{
protected:
	typename PFP::MAP& m_map ;
195
	VertexAttribute<typename PFP::VEC3>& m_position ;
196 197

public:
198
	LoopNormalisationAnalysisFilter(typename PFP::MAP& m, VertexAttribute<typename PFP::VEC3>& p) : m_map(m), m_position(p)
199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224
	{}

	void operator() ()
	{
		TraversorV<typename PFP::MAP> trav(m_map) ;
		for (Dart d = trav.begin(); d != trav.end(); d = trav.next())
		{
			if(m_map.isBoundaryVertex(d))
			{
				Dart db = m_map.findBoundaryFaceOfVertex(d);

				unsigned int degree = m_map.vertexDegreeOnBoundary(db) ;
				float n = 3.0/8.0 + 1.0/4.0 * cos(2.0 * M_PI / degree) ;
				n = 8.0/5.0 * (n * n) ;

				m_position[db] /= n ;
			}
		}
	}
} ;

template <typename PFP>
class LoopOddAnalysisFilter : public MRFilter
{
protected:
	typename PFP::MAP& m_map ;
225
	VertexAttribute<typename PFP::VEC3>& m_position ;
226 227

public:
228
	LoopOddAnalysisFilter(typename PFP::MAP& m, VertexAttribute<typename PFP::VEC3>& p) : m_map(m), m_position(p)
229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273
	{}

	void operator() ()
	{
		TraversorE<typename PFP::MAP> trav(m_map) ;
		for (Dart d = trav.begin(); d != trav.end(); d = trav.next())
		{
			if(m_map.isBoundaryEdge(d))
			{
				Dart db = m_map.findBoundaryFaceOfEdge(d);
				typename PFP::VEC3 p = loopOddVertex<PFP>(m_map, m_position, db) ;

				m_map.incCurrentLevel() ;

				Dart oddV = m_map.phi2(db) ;
				m_position[oddV] -= p ;

				m_map.decCurrentLevel() ;
			}
			else
			{
				typename PFP::VEC3 p = (m_position[d] + m_position[m_map.phi2(d)]) * typename PFP::REAL(0.5);

				m_map.incCurrentLevel() ;

				Dart midV = m_map.phi2(d) ;
				m_position[midV] -= p ;

				m_map.decCurrentLevel() ;
			}
		}
	}
} ;

/*********************************************************************************
 *                           SYNTHESIS FILTERS
 *********************************************************************************/

/* Linear Interpolation
 *********************************************************************************/
template <typename PFP>
class LerpEdgeSynthesisFilter : public MRFilter
{
protected:
	typename PFP::MAP& m_map ;
274
	VertexAttribute<typename PFP::VEC3>& m_position ;
275 276

public:
277
	LerpEdgeSynthesisFilter(typename PFP::MAP& m, VertexAttribute<typename PFP::VEC3>& p) : m_map(m), m_position(p)
278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301
	{}

	void operator() ()
	{
		TraversorE<typename PFP::MAP> trav(m_map) ;
		for (Dart d = trav.begin(); d != trav.end(); d = trav.next())
		{
			typename PFP::VEC3 p = (m_position[d] + m_position[m_map.phi2(d)]) * typename PFP::REAL(0.5);

			m_map.incCurrentLevel() ;

			Dart midV = m_map.phi2(d) ;
			m_position[midV] = p ;

			m_map.decCurrentLevel() ;
		}
	}
} ;

template <typename PFP>
class LerpFaceSynthesisFilter : public MRFilter
{
protected:
	typename PFP::MAP& m_map ;
302
	VertexAttribute<typename PFP::VEC3>& m_position ;
303 304

public:
305
	LerpFaceSynthesisFilter(typename PFP::MAP& m, VertexAttribute<typename PFP::VEC3>& p) : m_map(m), m_position(p)
306 307 308 309 310 311 312 313 314 315
	{}

	void operator() ()
	{
		TraversorF<typename PFP::MAP> trav(m_map) ;
		for (Dart d = trav.begin(); d != trav.end(); d = trav.next())
		{
			typename PFP::VEC3 p = Algo::Geometry::faceCentroid<PFP>(m_map, d, m_position);

			m_map.incCurrentLevel() ;
316 317 318 319 320
			if(m_map.faceDegree(d) != 3)
			{
				Dart midF = m_map.phi2(m_map.phi1(d));
				m_position[midF] = p ;
			}
321
			m_map.decCurrentLevel() ;
322

323 324 325 326 327 328 329 330 331
		}
	}
} ;

template <typename PFP>
class LerpVolumeSynthesisFilter : public MRFilter
{
protected:
	typename PFP::MAP& m_map ;
332
	VertexAttribute<typename PFP::VEC3>& m_position ;
333 334

public:
335
	LerpVolumeSynthesisFilter(typename PFP::MAP& m, VertexAttribute<typename PFP::VEC3>& p) : m_map(m), m_position(p)
336 337 338 339 340 341 342 343 344 345 346
	{}

	void operator() ()
	{
		TraversorW<typename PFP::MAP> trav(m_map) ;
		for (Dart d = trav.begin(); d != trav.end(); d = trav.next())
		{
			typename PFP::VEC3 p = Algo::Geometry::volumeCentroid<PFP>(m_map, d, m_position);

			m_map.incCurrentLevel() ;

347
			if(!Algo::Modelisation::Tetrahedralization::isTetrahedron<PFP>(m_map,d))
348
			{
349

350 351 352
				Dart midV = m_map.phi_1(m_map.phi2(m_map.phi1(d)));
				m_position[midV] = p ;
			}
353
			m_map.decCurrentLevel() ;
354

355 356 357 358
		}
	}
} ;

untereiner's avatar
untereiner committed
359
/* Ber02
360
 *********************************************************************************/
untereiner's avatar
untereiner committed
361 362

//w-lift(a)
363
template <typename PFP>
untereiner's avatar
untereiner committed
364
class Ber02OddSynthesisFilter : public MRFilter
365 366 367
{
protected:
	typename PFP::MAP& m_map ;
368
	VertexAttribute<typename PFP::VEC3>& m_position ;
369 370

public:
371
	Ber02OddSynthesisFilter(typename PFP::MAP& m, VertexAttribute<typename PFP::VEC3>& p) : m_map(m), m_position(p)
372 373 374 375
	{}

	void operator() ()
	{
untereiner's avatar
untereiner committed
376 377 378 379
		float a = 0.5;

		TraversorW<typename PFP::MAP> travW(m_map) ;
		for (Dart d = travW.begin(); d != travW.end(); d = travW.next())
380
		{
untereiner's avatar
untereiner committed
381 382 383 384 385 386
			typename PFP::VEC3 vc = Algo::Geometry::volumeCentroid<PFP>(m_map, d, m_position);

			unsigned int count = 0;
			typename PFP::VEC3 ec(0);
			Traversor3WE<typename PFP::MAP> travWE(m_map, d);
			for (Dart dit = travWE.begin(); dit != travWE.end(); dit = travWE.next())
387
			{
untereiner's avatar
untereiner committed
388 389 390 391
				m_map.incCurrentLevel();
				ec += m_position[m_map.phi2(dit)];
				m_map.decCurrentLevel();
				++count;
392
			}
untereiner's avatar
untereiner committed
393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411
			ec /= count;

			count = 0;
			typename PFP::VEC3 fc(0);
			Traversor3WF<typename PFP::MAP> travWF(m_map, d);
			for (Dart dit = travWF.begin(); dit != travWF.end(); dit = travWF.next())
			{
				m_map.incCurrentLevel();
				fc += m_position[m_map.phi2(m_map.phi1(dit))];
				m_map.decCurrentLevel();
				++count;
			}

			fc /= count;

			m_map.incCurrentLevel() ;
			Dart midV = m_map.phi_1(m_map.phi2(m_map.phi1(d)));
			m_position[midV] += 8 * a * a * a * vc + 12 * a * a * ec + 6 * a * fc;
			m_map.decCurrentLevel() ;
412
		}
untereiner's avatar
untereiner committed
413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447

		TraversorF<typename PFP::MAP> travF(m_map) ;
		for (Dart d = travF.begin(); d != travF.end(); d = travF.next())
		{
			typename PFP::VEC3 vf = Algo::Geometry::faceCentroid<PFP>(m_map, d, m_position);

			typename PFP::VEC3 ef(0);
			unsigned int count = 0;
			Traversor3FE<typename PFP::MAP> travFE(m_map, d);
			for (Dart dit = travFE.begin(); dit != travFE.end(); dit = travFE.next())
			{
				m_map.incCurrentLevel();
				ef += m_position[m_map.phi2(dit)];
				m_map.decCurrentLevel();
				++count;
			}
			ef /= count;

			m_map.incCurrentLevel() ;
			Dart midF = m_map.phi2(m_map.phi1(d));
			m_position[midF] += vf * 4.0 * a * a + ef * 4.0 * a;
			m_map.decCurrentLevel() ;
		}

		TraversorE<typename PFP::MAP> travE(m_map) ;
		for (Dart d = travE.begin(); d != travE.end(); d = travE.next())
		{
			typename PFP::VEC3 ve = (m_position[d] + m_position[m_map.phi2(d)]) * typename PFP::REAL(0.5);

			m_map.incCurrentLevel() ;
			Dart midV = m_map.phi2(d) ;
			m_position[midV] += ve * a * 2.0;
			m_map.decCurrentLevel() ;
		}

448 449 450
	}
} ;

untereiner's avatar
untereiner committed
451
// s-lift(a)
452
template <typename PFP>
untereiner's avatar
untereiner committed
453
class Ber02EvenSynthesisFilter : public MRFilter
454 455 456
{
protected:
	typename PFP::MAP& m_map ;
457
	VertexAttribute<typename PFP::VEC3>& m_position ;
458 459

public:
460
	Ber02EvenSynthesisFilter(typename PFP::MAP& m, VertexAttribute<typename PFP::VEC3>& p) : m_map(m), m_position(p)
461 462 463 464
	{}

	void operator() ()
	{
untereiner's avatar
untereiner committed
465 466 467 468
		float a = 0.5;

		TraversorV<typename PFP::MAP> travV(m_map);
		for(Dart d = travV.begin() ; d != travV.end() ; d = travV.next())
469
		{
untereiner's avatar
untereiner committed
470
			if(!m_map.isBoundaryVertex(d))
471
			{
untereiner's avatar
untereiner committed
472 473 474 475 476 477 478 479 480 481 482 483
				typename PFP::VEC3 cv(0);
				unsigned int count = 0;
				Traversor3VW<typename PFP::MAP> travVW(m_map,d);
				for(Dart dit = travVW.begin(); dit != travVW.end() ; dit = travVW.next())
				{
					m_map.incCurrentLevel() ;
					Dart midV = m_map.phi_1(m_map.phi2(m_map.phi1(dit)));
					cv += m_position[midV];
					m_map.decCurrentLevel() ;
					++count;
				}
				cv /= count;
484

untereiner's avatar
untereiner committed
485 486 487 488 489 490 491 492 493 494 495 496
				typename PFP::VEC3 fv(0);
				count = 0;
				Traversor3VF<typename PFP::MAP> travVF(m_map,d);
				for(Dart dit = travVF.begin(); dit != travVF.end() ; dit = travVF.next())
				{
					m_map.incCurrentLevel() ;
					Dart midV = m_map.phi2(m_map.phi1(dit));
					fv += m_position[midV];
					m_map.decCurrentLevel() ;
					++count;
				}
				fv /= count;
497

untereiner's avatar
untereiner committed
498 499 500 501 502 503 504 505 506 507 508 509 510 511
				typename PFP::VEC3 ev(0);
				count = 0;
				Traversor3VE<typename PFP::MAP> travVE(m_map,d);
				for(Dart dit = travVE.begin(); dit != travVE.end() ; dit = travVE.next())
				{
					m_map.incCurrentLevel() ;
					Dart midV = m_map.phi2(dit);
					ev += m_position[midV];
					m_map.decCurrentLevel() ;
					++count;
				}
				ev /= count;

				m_position[d] += cv * 8 * a * a * a + fv * 12 * a * a + ev * 6 * a;
512
			}
untereiner's avatar
untereiner committed
513
			else
514 515 516
			{
				Dart db = m_map.findBoundaryFaceOfVertex(d);

untereiner's avatar
untereiner committed
517 518 519 520 521 522 523 524 525 526 527 528
				typename PFP::VEC3 fv(0);
				unsigned int count = 0;
				Traversor2VF<typename PFP::MAP> travVF(m_map,db);
				for(Dart dit = travVF.begin(); dit != travVF.end() ; dit = travVF.next())
				{
					m_map.incCurrentLevel() ;
					Dart midV = m_map.phi2(m_map.phi1(dit));
					fv += m_position[midV];
					m_map.decCurrentLevel() ;
					++count;
				}
				fv /= count;
529

untereiner's avatar
untereiner committed
530 531 532 533 534 535 536 537 538 539 540 541
				typename PFP::VEC3 ev(0);
				count = 0;
				Traversor2VE<typename PFP::MAP> travVE(m_map,db);
				for(Dart dit = travVE.begin(); dit != travVE.end() ; dit = travVE.next())
				{
					m_map.incCurrentLevel() ;
					Dart midV = m_map.phi2(dit);
					ev += m_position[midV];
					m_map.decCurrentLevel() ;
					++count;
				}
				ev /= count;
542

untereiner's avatar
untereiner committed
543 544 545 546 547 548
				m_position[db] += fv * 4 * a * a + ev * 4 * a;
			}
		}

		TraversorE<typename PFP::MAP> travE(m_map);
		for(Dart d = travE.begin() ; d != travE.end() ; d = travE.next())
549 550 551 552
		{
			if(m_map.isBoundaryEdge(d))
			{
				Dart db = m_map.findBoundaryFaceOfEdge(d);
untereiner's avatar
untereiner committed
553 554

				typename PFP::VEC3 fe(0);
555 556

				m_map.incCurrentLevel() ;
untereiner's avatar
untereiner committed
557 558 559
				Dart midV = m_map.phi2(m_map.phi1(db));
				fe += m_position[midV];
				m_map.decCurrentLevel() ;
560

untereiner's avatar
untereiner committed
561 562 563 564 565 566
				m_map.incCurrentLevel() ;
				midV = m_map.phi2(m_map.phi1(m_map.phi2(db)));
				fe += m_position[midV];
				m_map.decCurrentLevel() ;

				fe /= 2;
567

untereiner's avatar
untereiner committed
568 569 570
				m_map.incCurrentLevel() ;
				Dart midF = m_map.phi2(db);
				m_position[midF] += fe * 2 * a;
571 572 573 574
				m_map.decCurrentLevel() ;
			}
			else
			{
untereiner's avatar
untereiner committed
575 576 577 578 579 580 581 582 583 584 585
				typename PFP::VEC3 ce(0);
				unsigned int count = 0;
				Traversor3EW<typename PFP::MAP> travEW(m_map, d);
				for(Dart dit = travEW.begin() ; dit != travEW.end() ; dit = travEW.next())
				{
					m_map.incCurrentLevel() ;
					Dart midV = m_map.phi_1(m_map.phi2(m_map.phi1(dit)));
					ce += m_position[midV];
					m_map.decCurrentLevel() ;
					++count;
				}
586

untereiner's avatar
untereiner committed
587
				ce /= count;
588

untereiner's avatar
untereiner committed
589 590 591 592 593 594 595 596 597 598 599 600 601
				typename PFP::VEC3 fe(0);
				count = 0;
				Traversor3FW<typename PFP::MAP> travFW(m_map, d);
				for(Dart dit = travFW.begin() ; dit != travFW.end() ; dit = travFW.next())
				{
					m_map.incCurrentLevel() ;
					Dart midV = m_map.phi2(m_map.phi1(dit));
					fe += m_position[midV];
					m_map.decCurrentLevel() ;
					++count;
				}

				fe /= count;
602

untereiner's avatar
untereiner committed
603 604 605
				m_map.incCurrentLevel() ;
				Dart midF = m_map.phi2(d);
				m_position[midF] += ce * 4 * a * a + fe * 4 * a;
606 607 608 609
				m_map.decCurrentLevel() ;
			}
		}

untereiner's avatar
untereiner committed
610 611
		TraversorF<typename PFP::MAP> travF(m_map) ;
		for (Dart d = travF.begin(); d != travF.end(); d = travF.next())
612
		{
untereiner's avatar
untereiner committed
613
			typename PFP::VEC3 cf(0);
614

untereiner's avatar
untereiner committed
615 616 617 618
			m_map.incCurrentLevel();
			Dart midV = m_map.phi_1(m_map.phi2(m_map.phi1(d)));
			cf += m_position[midV];
			m_map.decCurrentLevel();
619

untereiner's avatar
untereiner committed
620 621 622 623 624 625 626
			if(!m_map.isBoundaryFace(d))
			{
				Dart d3 = m_map.phi3(d);
				m_map.incCurrentLevel();
				Dart midV = m_map.phi_1(m_map.phi2(m_map.phi1(d3)));
				cf += m_position[midV];
				m_map.decCurrentLevel();
627

untereiner's avatar
untereiner committed
628
				cf /= 2;
629
			}
untereiner's avatar
untereiner committed
630 631 632 633 634

			m_map.incCurrentLevel() ;
			Dart midF = m_map.phi2(m_map.phi1(d));
			m_position[midF] += cf * 2 * a;
			m_map.decCurrentLevel() ;
635
		}
untereiner's avatar
untereiner committed
636

637 638 639
	}
} ;

untereiner's avatar
untereiner committed
640
// s-scale(a)
641
template <typename PFP>
untereiner's avatar
untereiner committed
642
class Ber02ScaleSynthesisFilter : public MRFilter
643 644 645
{
protected:
	typename PFP::MAP& m_map ;
646
	VertexAttribute<typename PFP::VEC3>& m_position ;
647 648

public:
649
	Ber02ScaleSynthesisFilter(typename PFP::MAP& m, VertexAttribute<typename PFP::VEC3>& p) : m_map(m), m_position(p)
650 651 652 653
	{}

	void operator() ()
	{
untereiner's avatar
untereiner committed
654 655 656 657
		float a = 0.5;

		TraversorV<typename PFP::MAP> travV(m_map) ;
		for (Dart d = travV.begin(); d != travV.end(); d = travV.next())
658 659 660 661
		{
			if(m_map.isBoundaryVertex(d))
			{
				Dart db = m_map.findBoundaryFaceOfVertex(d);
untereiner's avatar
untereiner committed
662 663 664 665 666 667 668
				m_position[db] *= a * a;
			}
			else
			{
				m_position[d] *= a * a * a;
			}
		}
669

untereiner's avatar
untereiner committed
670 671 672 673 674 675
		TraversorE<typename PFP::MAP> travE(m_map) ;
		for (Dart d = travE.begin(); d != travE.end(); d = travE.next())
		{
			if(m_map.isBoundaryEdge(d))
			{
				Dart db = m_map.findBoundaryFaceOfEdge(d);
676

untereiner's avatar
untereiner committed
677 678 679 680 681 682 683 684 685 686 687 688 689
				m_map.incCurrentLevel() ;
				Dart midE = m_map.phi2(db);
				m_position[midE] *= a ;
				m_map.decCurrentLevel() ;
			}
			else
			{
				m_map.incCurrentLevel() ;
				Dart midE = m_map.phi2(d);
				m_position[midE] *= a * a;
				m_map.decCurrentLevel() ;
			}
		}
690

untereiner's avatar
untereiner committed
691 692 693 694 695 696 697 698 699
		TraversorF<typename PFP::MAP> travF(m_map) ;
		for (Dart d = travF.begin(); d != travF.end(); d = travF.next())
		{
			if(!m_map.isBoundaryFace(d))
			{
				m_map.incCurrentLevel() ;
				Dart midF = m_map.phi2(m_map.phi1(d));
				m_position[midF] *= a ;
				m_map.decCurrentLevel() ;
700 701 702 703 704
			}
		}
	}
} ;

untereiner's avatar
untereiner committed
705
/* Loop on Boundary Vertices and SHW04 on Insides Vertices
706 707 708 709 710 711
 *********************************************************************************/
template <typename PFP>
class LoopOddSynthesisFilter : public MRFilter
{
protected:
	typename PFP::MAP& m_map ;
712
	VertexAttribute<typename PFP::VEC3>& m_position ;
713 714

public:
715
	LoopOddSynthesisFilter(typename PFP::MAP& m, VertexAttribute<typename PFP::VEC3>& p) : m_map(m), m_position(p)
716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750
	{}

	void operator() ()
	{
		TraversorE<typename PFP::MAP> trav(m_map) ;
		for (Dart d = trav.begin(); d != trav.end(); d = trav.next())
		{
			if(m_map.isBoundaryEdge(d))
			{
				Dart db = m_map.findBoundaryFaceOfEdge(d);
				typename PFP::VEC3 p = loopOddVertex<PFP>(m_map, m_position, db) ;

				m_map.incCurrentLevel() ;

				Dart oddV = m_map.phi2(db) ;
				m_position[oddV] += p ;

				m_map.decCurrentLevel() ;
			}
			else
			{
				typename PFP::VEC3 p = (m_position[d] + m_position[m_map.phi2(d)]) * typename PFP::REAL(0.5);

				m_map.incCurrentLevel() ;

				Dart midV = m_map.phi2(d) ;
				m_position[midV] += p ;

				m_map.decCurrentLevel() ;
			}
		}
	}
} ;

template <typename PFP>
untereiner's avatar
untereiner committed
751
class LoopNormalisationSynthesisFilter : public MRFilter
752 753 754
{
protected:
	typename PFP::MAP& m_map ;
755
	VertexAttribute<typename PFP::VEC3>& m_position ;
756 757

public:
758
	LoopNormalisationSynthesisFilter(typename PFP::MAP& m, VertexAttribute<typename PFP::VEC3>& p) : m_map(m), m_position(p)
759 760
	{}

untereiner's avatar
untereiner committed
761
	void operator() ()
762
	{
untereiner's avatar
untereiner committed
763 764 765 766 767 768
		TraversorV<typename PFP::MAP> trav(m_map) ;
		for (Dart d = trav.begin(); d != trav.end(); d = trav.next())
		{
			if(m_map.isBoundaryVertex(d))
			{
				Dart db = m_map.findBoundaryFaceOfVertex(d);
769

untereiner's avatar
untereiner committed
770 771 772
				unsigned int degree = m_map.vertexDegreeOnBoundary(db) ;
				float n = 3.0/8.0 + 1.0/4.0 * cos(2.0 * M_PI / degree) ;
				n = 8.0/5.0 * (n * n) ;
773

untereiner's avatar
untereiner committed
774 775 776
				m_position[db] *= n ;
			}
		}
777 778 779 780
	}
} ;

template <typename PFP>
untereiner's avatar
untereiner committed
781
class LoopEvenSynthesisFilter : public MRFilter
782 783 784
{
protected:
	typename PFP::MAP& m_map ;
785
	VertexAttribute<typename PFP::VEC3>& m_position ;
786 787

public:
788
	LoopEvenSynthesisFilter(typename PFP::MAP& m, VertexAttribute<typename PFP::VEC3>& p) : m_map(m), m_position(p)
789 790
	{}

untereiner's avatar
untereiner committed
791
	void operator() ()
792
	{
untereiner's avatar
untereiner committed
793 794 795 796 797 798 799 800 801 802
		TraversorV<typename PFP::MAP> trav(m_map) ;
		for (Dart d = trav.begin(); d != trav.end(); d = trav.next())
		{
			if(m_map.isBoundaryVertex(d))
			{
				Dart db = m_map.findBoundaryFaceOfVertex(d);
				typename PFP::VEC3 p = loopEvenVertex<PFP>(m_map, m_position, db) ;
				m_position[db] += p ;
			}
		}
803 804 805 806 807 808 809 810
	}
} ;

template <typename PFP>
class LoopVolumeSynthesisFilter : public MRFilter
{
protected:
	typename PFP::MAP& m_map ;
811
	VertexAttribute<typename PFP::VEC3>& m_position ;
812 813

public:
814
	LoopVolumeSynthesisFilter(typename PFP::MAP& m, VertexAttribute<typename PFP::VEC3>& p) : m_map(m), m_position(p)
815 816 817 818 819 820 821
	{}

	void operator() ()
	{
		TraversorW<typename PFP::MAP> trav(m_map) ;
		for (Dart d = trav.begin(); d != trav.end(); d = trav.next())
		{
822
			if(!Algo::Modelisation::Tetrahedralization::isTetrahedron<PFP>(m_map,d))
823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841
			{
				typename PFP::VEC3 p = Algo::Geometry::volumeCentroid<PFP>(m_map, d, m_position);

				m_map.incCurrentLevel() ;

				Dart midV = m_map.phi_1(m_map.phi2(m_map.phi1(d)));
				m_position[midV] = p ;

				m_map.decCurrentLevel() ;
			}
		}
	}
} ;

template <typename PFP>
class SHW04VolumeNormalisationSynthesisFilter : public MRFilter
{
protected:
	typename PFP::MAP& m_map ;
842
	VertexAttribute<typename PFP::VEC3>& m_position ;
843 844

public:
845
	SHW04VolumeNormalisationSynthesisFilter(typename PFP::MAP& m, VertexAttribute<typename PFP::VEC3>& p) : m_map(m), m_position(p)
846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874
	{}

	void operator() ()
	{
		m_map.incCurrentLevel() ;
		TraversorV<typename PFP::MAP> trav(m_map) ;
		for (Dart d = trav.begin(); d != trav.end(); d = trav.next())
		{
			if(!m_map.isBoundaryVertex(d))
			{
				typename PFP::VEC3 p = typename PFP::VEC3(0);
				unsigned int degree = 0;

				Traversor3VW<typename PFP::MAP> travVW(m_map, d);
				for(Dart dit = travVW.begin() ; dit != travVW.end() ; dit = travVW.next())
				{
					p += SHW04Vertex<PFP>(m_map, m_position, dit);
					++degree;
				}

				p /= degree;

				m_position[d] = p ;
			}
		}
		m_map.decCurrentLevel() ;
	}
} ;

untereiner's avatar
untereiner committed
875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137 1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 1164 1165 1166
/* Catmull-clark on Boundary Vertices and MJ96 on Insides Vertices
 *********************************************************************************/
template <typename PFP>
class MJ96VertexSubdivision : public MRFilter
{
protected:
	typename PFP::MAP& m_map ;
	VertexAttribute<typename PFP::VEC3>& m_position ;

public:
	MJ96VertexSubdivision(typename PFP::MAP& m, VertexAttribute<typename PFP::VEC3>& p) : m_map(m), m_position(p)
	{}

	void operator() ()
	{
		TraversorV<typename PFP::MAP> trav(m_map) ;
		for (Dart d = trav.begin(); d != trav.end(); d = trav.next())
		{
			if(m_map.isBoundaryVertex(d))
			{
				Dart db = m_map.findBoundaryFaceOfVertex(d);

				typename PFP::VEC3 np1(0) ;
				typename PFP::VEC3 np2(0) ;
				unsigned int degree1 = 0 ;
				unsigned int degree2 = 0 ;
				Dart it = db ;
				do
				{
					++degree1 ;
					Dart dd = m_map.phi1(it) ;
					np1 += m_position[dd] ;
					Dart end = m_map.phi_1(it) ;
					dd = m_map.phi1(dd) ;
					do
					{
						++degree2 ;
						np2 += m_position[dd] ;
						dd = m_map.phi1(dd) ;
					} while(dd != end) ;
					it = m_map.phi2(m_map.phi_1(it)) ;
				} while(it != db) ;

				float beta = 3.0 / (2.0 * degree1) ;
				float gamma = 1.0 / (4.0 * degree2) ;
				np1 *= beta / degree1 ;
				np2 *= gamma / degree2 ;

				typename PFP::VEC3 vp = m_position[db] ;
				vp *= 1.0 - beta - gamma ;

				m_map.incCurrentLevel() ;

				m_position[d] = np1 + np2 + vp ;

				m_map.decCurrentLevel() ;

			}
			else
			{
				typename PFP::VEC3 P = m_position[d];

				//vertex points
				typename PFP::VEC3 Cavg = typename PFP::VEC3(0);
				unsigned int degree = 0;
				Traversor3VW<typename PFP::MAP> travVW(m_map, d);
				for(Dart dit = travVW.begin() ; dit != travVW.end() ; dit = travVW.next())
				{
					Cavg += Algo::Geometry::volumeCentroid<PFP>(m_map, dit, m_position);
					++degree;
				}
				Cavg /= degree;

				typename PFP::VEC3 Aavg = typename PFP::VEC3(0);
				degree = 0;
				Traversor3VF<typename PFP::MAP> travVF(m_map, d);
				for(Dart dit = travVF.begin() ; dit != travVF.end() ; dit = travVF.next())
				{
					Aavg += Algo::Geometry::faceCentroid<PFP>(m_map, dit, m_position);
					++degree;
				}
				Aavg /= degree;

				typename PFP::VEC3 Mavg = typename PFP::VEC3(0);
				degree = 0;
				Traversor3VE<typename PFP::MAP> travVE(m_map, d);
				for(Dart dit = travVE.begin() ; dit != travVE.end() ; dit = travVE.next())
				{
					Dart d2 = m_map.phi2(dit);
					Aavg += (m_position[dit] + m_position[d2]) * typename PFP::REAL(0.5);
					++degree;
				}
				Aavg /= degree;

				typename PFP::VEC3 vp = Cavg + Aavg * 3 + Mavg * 3 + P;
				vp /= 8;

				m_map.incCurrentLevel() ;

				m_position[d] = P;//vp;

				m_map.decCurrentLevel() ;
			}
		}
	}
};

template <typename PFP>
class MJ96EdgeSubdivision : public MRFilter
{
protected:
	typename PFP::MAP& m_map ;
	VertexAttribute<typename PFP::VEC3>& m_position ;

public:
	MJ96EdgeSubdivision(typename PFP::MAP& m, VertexAttribute<typename PFP::VEC3>& p) : m_map(m), m_position(p)
	{}

	void operator() ()
	{
		TraversorE<typename PFP::MAP> trav(m_map) ;
		for (Dart d = trav.begin(); d != trav.end(); d = trav.next())
		{
			if(m_map.isBoundaryEdge(d))
			{
				Dart db = m_map.findBoundaryFaceOfEdge(d);

				Dart d1 = m_map.phi2(db) ;
				Dart d2 = m_map.phi2(d1) ;
				Dart d3 = m_map.phi_1(d1) ;
				Dart d4 = m_map.phi_1(d2) ;
				Dart d5 = m_map.phi1(m_map.phi1(d1)) ;
				Dart d6 = m_map.phi1(m_map.phi1(d2)) ;

				typename PFP::VEC3 p1 = m_position[d1] ;
				typename PFP::VEC3 p2 = m_position[d2] ;
				typename PFP::VEC3 p3 = m_position[d3] ;
				typename PFP::VEC3 p4 = m_position[d4] ;
				typename PFP::VEC3 p5 = m_position[d5] ;
				typename PFP::VEC3 p6 = m_position[d6] ;

				p1 *= 3.0 / 8.0 ;
				p2 *= 3.0 / 8.0 ;
				p3 *= 1.0 / 16.0 ;
				p4 *= 1.0 / 16.0 ;
				p5 *= 1.0 / 16.0 ;
				p6 *= 1.0 / 16.0 ;

				m_map.incCurrentLevel() ;

				Dart midV = m_map.phi2(d);

				m_position[midV] = p1 + p2 + p3 + p4 + p5 + p6 ;

				m_map.decCurrentLevel() ;
			}
			else
			{
				//edge points
				typename PFP::VEC3 Cavg = typename PFP::VEC3(0);
				unsigned int degree = 0;
				Traversor3EW<typename PFP::MAP> travEW(m_map, d);
				for(Dart dit = travEW.begin() ; dit != travEW.end() ; dit = travEW.next())
				{
					Cavg += Algo::Geometry::volumeCentroid<PFP>(m_map, dit, m_position);
					++degree;
				}
				Cavg /= degree;

				typename PFP::VEC3 Aavg = typename PFP::VEC3(0);
				degree = 0;
				Traversor3EF<typename PFP::MAP> travEF(m_map, d);
				for(Dart dit = travEF.begin() ; dit != travEF.end() ; dit = travEF.next())
				{
					Aavg += Algo::Geometry::faceCentroid<PFP>(m_map, dit, m_position);
					++degree;
				}
				Aavg /= degree;

				Dart d2 = m_map.phi2(d);
				typename PFP::VEC3 M = (m_position[d] + m_position[d2]) * typename PFP::REAL(0.5);

				typename PFP::VEC3 ep = Cavg + Aavg * 2 + M * (degree - 3);
				ep /= degree;

				m_map.incCurrentLevel() ;

				Dart midV = m_map.phi2(d);

				m_position[midV] = ep;

				m_map.decCurrentLevel() ;
			}
		}
	}
};

template <typename PFP>
class MJ96FaceSubdivision : public MRFilter
{
protected:
	typename PFP::MAP& m_map ;
	VertexAttribute<typename PFP::VEC3>& m_position ;

public:
	MJ96FaceSubdivision(typename PFP::MAP& m, VertexAttribute<typename PFP::VEC3>& p) : m_map(m), m_position(p)
	{}

	void operator() ()
	{
		TraversorF<typename PFP::MAP> trav(m_map) ;
		for (Dart d = trav.begin(); d != trav.end(); d = trav.next())
		{
			if(m_map.isBoundaryFace(d))
			{
				Dart db = m_map.phi3(d);

				typename PFP::VEC3 p(0) ;
				unsigned int degree = 0 ;
				Traversor2FV<typename PFP::MAP> trav(m_map, db) ;
				for(Dart it = trav.begin(); it != trav.end(); it = trav.next())
				{
					++degree ;
					p += m_position[it] ;
				}
				p /= degree ;

				m_map.incCurrentLevel() ;

				Dart df = m_map.phi1(m_map.phi1(d)) ;

				m_position[df] = p ;

				m_map.decCurrentLevel() ;
			}
			else
			{
				//face points
				typename PFP::VEC3 C0 = Algo::Geometry::volumeCentroid<PFP>(m_map, d, m_position);
				typename PFP::VEC3 C1 = Algo::Geometry::volumeCentroid<PFP>(m_map, m_map.phi3(d), m_position);

				typename PFP::VEC3 A = Algo::Geometry::faceCentroid<PFP>(m_map, m_map.phi3(d), m_position);

				typename PFP::VEC3 fp = C0 + A * 2 + C1;
				fp /= 4;

				m_map.incCurrentLevel() ;

				Dart df = m_map.phi1(m_map.phi1(d)) ;
				m_position[df] = fp;

				m_map.decCurrentLevel() ;
			}
		}
	}
};

template <typename PFP>
class MJ96VolumeSubdivision : public MRFilter
{
protected:
	typename PFP::MAP& m_map ;
	VertexAttribute<typename PFP::VEC3>& m_position ;

public:
	MJ96VolumeSubdivision(typename PFP::MAP& m, VertexAttribute<typename PFP::VEC3>& p) : m_map(m), m_position(p)
	{}

	void operator() ()
	{
		TraversorW<typename PFP::MAP> trav(m_map) ;
		for (Dart d = trav.begin(); d != trav.end(); d = trav.next())
		{
			//cell points : these points are the average of the
			//vertices of the lattice that bound the cell
			typename PFP::VEC3 p = Algo::Geometry::volumeCentroid<PFP>(m_map, d, m_position);

			m_map.incCurrentLevel() ;

			if(!Algo::Modelisation::Tetrahedralization::isTetrahedron<PFP>(m_map,d))
			{
				Dart midV = m_map.phi_1(m_map.phi2(m_map.phi1(d)));
				m_position[midV] = p ;
			}

			m_map.decCurrentLevel() ;

		}
	}
};


1167 1168 1169

} // namespace Multiresolution

1170 1171
} // namespace Algo

1172 1173 1174 1175
} // namespace CGoGN


#endif /* __3MR_FILTERS_PRIMAL__ */