Coupure prévue mardi 3 Août au matin pour maintenance du serveur. Nous faisons au mieux pour que celle-ci soit la plus brève possible.

subdivision.hpp 22.9 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
* Contact information: cgogn@unistra.fr                                        *
*                                                                              *
*******************************************************************************/

#include "Algo/Geometry/centroid.h"
#include "Topology/generic/autoAttributeHandler.h"

thery's avatar
thery committed
28
29
30
#define _USE_MATH_DEFINES
#include <math.h>

Pierre Kraemer's avatar
Pierre Kraemer committed
31
32
33
34
35
36
37
38
39
40
41
42
43
44
namespace CGoGN
{

namespace Algo
{

namespace Modelisation
{

template <typename PFP>
Dart trianguleFace(typename PFP::MAP& map, Dart d)
{
	Dart d1 = map.phi1(d);
	if (d1 == d)
45
		CGoGNout << "Warning: triangulation of a face with only one edge" << CGoGNendl;
Pierre Kraemer's avatar
Pierre Kraemer committed
46
	if (map.phi1(d1) == d)
47
		CGoGNout << "Warning: triangulation of a face with only two edges" << CGoGNendl;
Pierre Kraemer's avatar
Pierre Kraemer committed
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76

	map.splitFace(d, d1) ;
	map.cutEdge(map.phi_1(d)) ;
	Dart x = map.phi2(map.phi_1(d)) ;
	Dart dd = map.template phi<111>(x) ;
	while(dd != x)
	{
		Dart next = map.phi1(dd) ;
		map.splitFace(dd, map.phi1(x)) ;
		dd = next ;
	}

	return map.phi2(x);	// Return a dart of the central vertex
}

template <typename PFP, typename EMBV, typename EMB>
void trianguleFaces(typename PFP::MAP& map, EMBV& attributs, const FunctorSelect& selected)
{
	DartMarker m(map) ;
	for (Dart d = map.begin(); d != map.end(); map.next(d))
	{
		if (selected(d) && !m.isMarked(d))
		{
			EMB center = Algo::Geometry::faceCentroidGen<PFP,EMBV,EMB>(map, d, attributs);	// compute center
			Dart cd = trianguleFace<PFP>(map, d);	// triangule the face
			attributs[cd] = center;					// affect the data to the central vertex
			Dart fit = cd ;
			do
			{
77
				m.markOrbit(FACE, fit);
Pierre Kraemer's avatar
Pierre Kraemer committed
78
79
80
81
82
83
				fit = map.alpha1(fit);
			} while(fit != cd);
		}
	}
}

84
85
86
87
88
89
template <typename PFP>
void trianguleFaces(typename PFP::MAP& map, typename PFP::TVEC3& position, const FunctorSelect& selected)
{
	trianguleFaces<PFP, typename PFP::TVEC3, typename PFP::VEC3>(map, position, selected) ;
}

Pierre Kraemer's avatar
Pierre Kraemer committed
90
91
92
93
94
95
template <typename PFP>
void trianguleFaces(
		typename PFP::MAP& map,
		typename PFP::TVEC3& position, typename PFP::TVEC3& positionF,
		const FunctorSelect& selected)
{
96
97
	assert(position.getOrbit() == VERTEX) ;
	assert(positionF.getOrbit() == FACE) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
98
99
100
101
102
103
104
105
106
107
108
109

	DartMarker m(map) ;
	for (Dart d = map.begin(); d != map.end(); map.next(d))
	{
		if (selected(d) && !m.isMarked(d))
		{
			typename PFP::VEC3 p = positionF[d] ;
			Dart cd = trianguleFace<PFP>(map, d) ;	// triangule the face
			position[cd] = p ;						// affect the data to the central vertex
			Dart fit = cd ;
			do
			{
110
				m.markOrbit(FACE, fit);
Pierre Kraemer's avatar
Pierre Kraemer committed
111
112
113
114
115
116
				fit = map.alpha1(fit);
			} while(fit != cd);
		}
	}
}

Pierre Kraemer's avatar
Pierre Kraemer committed
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
template <typename PFP>
Dart quadranguleFace(typename PFP::MAP& map, Dart d)
{
	d = map.phi1(d) ;
	map.splitFace(d, map.template phi<11>(d)) ;
	map.cutEdge(map.phi_1(d)) ;
	Dart x = map.phi2(map.phi_1(d)) ;
	Dart dd = map.template phi<1111>(x) ;
	while(dd != x)
	{
		Dart next = map.template phi<11>(dd) ;
		map.splitFace(dd, map.phi1(x)) ;
		dd = next ;
	}

	return map.phi2(x);	// Return a dart of the central vertex
}

template <typename PFP, typename EMBV, typename EMB>
void quadranguleFaces(typename PFP::MAP& map, EMBV& attributs, const FunctorSelect& selected)
{
	DartMarker me(map) ;
	DartMarker mf(map) ;

	// first pass: cut the edges
	for (Dart d = map.begin(); d != map.end(); map.next(d))
	{
		if (selected(d) && !me.isMarked(d))
		{
			Dart f = map.phi1(d);
			map.cutEdge(d);
			Dart e = map.phi1(d);
//			TODO trouver pourquoi lerp bug avec ECell
//			attributs[m] = AttribOps::lerp<EMB,PFP>(attributs[d],attributs[f], 0.5);
			attributs[e] = attributs[d];
			attributs[e] += attributs[f];
			attributs[e] *= 0.5;

155
156
157
			me.markOrbit(EDGE, d);
			me.markOrbit(EDGE, e);
			mf.markOrbit(VERTEX, e);
Pierre Kraemer's avatar
Pierre Kraemer committed
158
159
160
161
162
163
164
165
166
167
168
169
170
171
		}
	}

	// second pass: quandrangule faces
	for (Dart d = map.begin(); d != map.end(); map.next(d))
	{
		if (selected(d) && !mf.isMarked(d))
		{
			EMB center = Algo::Geometry::faceCentroidGen<PFP,EMBV,EMB>(map, d, attributs);	// compute center
			Dart cf = quadranguleFace<PFP>(map, d);	// quadrangule the face
			attributs[cf] = center;					// affect the data to the central vertex
			Dart e = cf;
			do
			{
172
				mf.markOrbit(FACE, e);
Pierre Kraemer's avatar
Pierre Kraemer committed
173
174
175
176
177
178
				e = map.alpha1(e);
			} while (e != cf);
		}
	}
}

179
180
181
182
183
184
template <typename PFP>
void quadranguleFaces(typename PFP::MAP& map, typename PFP::TVEC3& position, const FunctorSelect& selected)
{
	quadranguleFaces<PFP, typename PFP::TVEC3, typename PFP::VEC3>(map, position, selected) ;
}

Pierre Kraemer's avatar
Pierre Kraemer committed
185
template <typename PFP, typename EMBV, typename EMB>
186
void CatmullClarkSubdivision(typename PFP::MAP& map, EMBV& attributs, const FunctorSelect& selected)
Pierre Kraemer's avatar
Pierre Kraemer committed
187
188
189
190
{
	std::vector<Dart> l_middles;
	std::vector<Dart> l_verts;

191
	CellMarkerNoUnmark m0(map, VERTEX);
Pierre Kraemer's avatar
Pierre Kraemer committed
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
	DartMarkerNoUnmark mf(map);
	DartMarkerNoUnmark me(map);

	// first pass: cut edges
	for (Dart d = map.begin(); d != map.end(); map.next(d))
	{
		if (selected(d) && !me.isMarked(d))
		{
			if (!m0.isMarked(d))
			{
				m0.mark(d);
				l_verts.push_back(d);
			}
			if (!m0.isMarked(map.phi2(d)))
			{
				m0.mark(map.phi2(d));
				l_verts.push_back(map.phi2(d));
			}

			Dart f = map.phi1(d);
			map.cutEdge(d);
			Dart e = map.phi1(d) ;

			attributs[e] =  attributs[d];
			attributs[e] += attributs[f];
			attributs[e] *= 0.5;

219
220
			me.markOrbit(EDGE, d);
			me.markOrbit(EDGE, e);
Pierre Kraemer's avatar
Pierre Kraemer committed
221
222
223
224
225
226
227
228
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
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305

			mf.mark(d) ;
			mf.mark(map.phi2(e)) ;

			l_middles.push_back(e);
		}
	}

	// second pass: quandrangule faces
	for (Dart d = map.begin(); d != map.end(); map.next(d))
	{
		if (selected(d) && mf.isMarked(d)) // for each face not subdivided
		{
			// compute center skip darts of new vertices non embedded
			EMB center = AttribOps::zero<EMB,PFP>();
			unsigned int count = 0 ;
			Dart it = d;
			do
			{
				mf.unmark(it);
				me.unmark(it);
				me.unmark(map.phi1(it));
				center += attributs[it];
				++count ;
				it = map.template phi<11>(it) ;
			} while(it != d) ;
			center /= double(count);
			Dart cf = quadranguleFace<PFP>(map, d);	// quadrangule the face
			attributs[cf] = center;					// affect the data to the central vertex
		}
	}

	// Compute edge points
	for(typename std::vector<Dart>::iterator mid = l_middles.begin(); mid != l_middles.end(); ++mid)
	{
		Dart x = *mid;
		// other side of the edge
		Dart y = map.phi2(x);
		if (y != x)
		{
			Dart f1 = map.phi_1(x);
			Dart f2 = map.phi2(map.phi1(y));
			EMB temp = AttribOps::zero<EMB,PFP>();
			temp =  attributs[f1];
			temp += attributs[f2];			// E' = (V0+V1+F1+F2)/4
			temp *= 0.25;
			attributs[x] *= 0.5;
			attributs[x] += temp;
		}
		// else nothing to do point already in the middle of segment
	}

	// Compute vertex points
	for(typename std::vector<Dart>::iterator vert = l_verts.begin(); vert != l_verts.end(); ++vert)
	{
		m0.unmark(*vert);

		EMB temp = AttribOps::zero<EMB,PFP>();
		EMB temp2 = AttribOps::zero<EMB,PFP>();

		unsigned int n = 0;
		Dart x = *vert;
		do
		{
			Dart m = map.phi1(x);
			Dart f = map.phi2(m);
			Dart v = map.template phi<11>(f);

			temp += attributs[f];
			temp2 += attributs[v];

			++n;
			x = map.alpha1(x);
		} while (x != *vert);

		EMB emcp = attributs[*vert];
		emcp *= double((n-2)*n);		// V' = (n-2)/n*V + 1/n2 *(F+E)
		emcp += temp;
		emcp += temp2;
		emcp /= double(n*n);

		attributs[*vert] = emcp ;
	}
}

306
307
308
309
310
311
template <typename PFP>
void CatmullClarkSubdivision(typename PFP::MAP& map, typename PFP::TVEC3& position, const FunctorSelect& selected)
{
	CatmullClarkSubdivision<PFP, typename PFP::TVEC3, typename PFP::VEC3>(map, position, selected) ;
}

Pierre Kraemer's avatar
Pierre Kraemer committed
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
inline double betaF(unsigned int n)
{
	switch(n)
	{
		case 1: return 0.234375 ;
		case 2: return 0.609375 ;
		case 3: return 0.5625 ;
		case 4: return 0.484375 ;
		case 5: return 0.420466 ;
		case 6: return 0.375 ;
		case 7: return 0.343174 ;
		case 8: return 0.320542 ;
		case 9: return 0.304065 ;
		case 10: return 0.291778 ;
		case 11: return 0.282408 ;
		case 12: return 0.27512 ;
		case 13: return 0.26935 ;
		case 14: return 0.264709 ;
		default:
			double t = 3.0 + 2.0 * cos((2.0*M_PI)/double(n)) ;
			return 5.0/8.0 - (t * t) / 64.0 ;
	}
}

template <typename PFP, typename EMBV, typename EMB>
void LoopSubdivision(typename PFP::MAP& map, EMBV& attributs, const FunctorSelect& selected)
{
	std::vector<Dart> l_middles;
	std::vector<Dart> l_verts;

342
	CellMarkerNoUnmark m0(map, VERTEX);
Pierre Kraemer's avatar
Pierre Kraemer committed
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
	DartMarkerNoUnmark mv(map);
	DartMarkerNoUnmark me(map);

	// first pass cut edges
	for (Dart d = map.begin(); d != map.end(); map.next(d))
	{
		if (selected(d) && !me.isMarked(d))
		{
			if (!m0.isMarked(d))
			{
				m0.mark(d);
				l_verts.push_back(d);
			}
			if (!m0.isMarked(map.phi2(d)))
			{
				m0.mark(map.phi2(d));
				l_verts.push_back(map.phi2(d));
			}

			Dart f = map.phi1(d);

			map.cutEdge(d);
			Dart e = map.phi1(d) ;

			attributs[e] =  attributs[d];
			attributs[e] += attributs[f];
			attributs[e] *= 0.5;

371
372
			me.markOrbit(EDGE, d);
			me.markOrbit(EDGE, e);
Pierre Kraemer's avatar
Pierre Kraemer committed
373

374
			mv.markOrbit(VERTEX, e);
Pierre Kraemer's avatar
Pierre Kraemer committed
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
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
448
449
450
451
452
453
454
455
456
457
458
459
460

			l_middles.push_back(e);
		}
	}

	// Compute edge points
	for(typename std::vector<Dart>::iterator mid = l_middles.begin(); mid != l_middles.end(); ++mid)
	{
		Dart d = *mid;
		Dart dd = map.phi2(d);
		if (dd != d)
		{
			attributs[d] *= 0.75;
			Dart e1 = map.template phi<111>(d);
			EMB temp = AttribOps::zero<EMB,PFP>();
			temp = attributs[e1];
			e1 = map.phi_1(map.phi_1(dd));
			temp += attributs[e1];
			temp *= 1.0/8.0;
			attributs[d] += temp;
		}
		// else nothing to do point already in the middle of segment
	}

	// Compute vertex points
	for(typename std::vector<Dart>::iterator vert = l_verts.begin(); vert != l_verts.end(); ++vert)
	{
		m0.unmark(*vert);

		EMB temp = AttribOps::zero<EMB,PFP>();
		int n = 0;
		Dart x = *vert;
		do
		{
			Dart y = map.phi1(map.phi1(x));
			temp += attributs[y];
			++n;
			x = map.alpha1(x);
		} while ((x != *vert));
		EMB emcp = attributs[*vert];
		if (n == 6)
		{
			temp /= 16.0;
			emcp *= 10.0/16.0;
			emcp += temp;
		}
		else
		{
			double beta = betaF(n) ;
			temp *= (beta / double(n));
			emcp *= (1.0 - beta);
			emcp += temp;
		}
		attributs[*vert] = emcp;
	}

	// insert new edges
	for (Dart d = map.begin(); d != map.end(); map.next(d))
	{
		if (mv.isMarked(d))
		{
			// unmark the darts of the face
			Dart fit = d ;
			do
			{
				mv.unmark(fit) ;
				me.unmark(fit) ;
				me.unmark(map.phi_1(fit)) ;
				fit = map.template phi<11>(fit) ;
			} while(fit != d) ;

			Dart dd = d;
			Dart e = map.template phi<11>(dd) ;
			map.splitFace(dd, e);

			dd = e;
			e = map.template phi<11>(dd) ;
			map.splitFace(dd, e);

			dd = e;
			e = map.template phi<11>(dd) ;
			map.splitFace(dd, e);
		}
	}
}

461
462
463
template <typename PFP, typename EMBV, typename EMB>
void TwoNPlusOneSubdivision(typename PFP::MAP& map, EMBV& attributs, const FunctorSelect& selected)
{
464
465
	AutoAttributeHandler< Dart > tablePred(map,EDGE);
	CellMarker m0(map, EDGE);
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497

	std::vector<Dart> dOrig;
	std::vector<EMB> cOrig;

	//first pass cut edge
	for (Dart d = map.begin(); d != map.end(); map.next(d))
	{
		if(selected(d))
		{
			if(!m0.isMarked(d)) {
				dOrig.push_back(d);

				Dart dd = d;
				do {
					if(!m0.isMarked(dd)) {
						EMB e1 = attributs[dd];
						EMB e2 = attributs[map.phi1(dd)];
						map.cutEdge(dd);
						attributs[map.phi1(dd)] = e1*2.0f/3.0f+e2/3.0f;
						map.cutEdge(map.phi1(dd));
						attributs[map.phi1(map.phi1(dd))] = e2*2.0f/3.0f+e1/3.0f;
						m0.mark(dd);
						m0.mark(map.phi1(dd));
						m0.mark(map.template phi<11>(dd));
					}
					dd = map.template phi<111>(dd);
				} while(dd!=d);

			}
		}
	}

498
	CGoGNout << "nb orig : " << dOrig.size() << CGoGNendl;
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528

	DartMarkerNoUnmark mCorner(map);
//	//second pass create corner face
	for (std::vector<Dart>::iterator it = dOrig.begin(); it != dOrig.end(); ++it)
	{
		EMB c = Algo::Geometry::faceCentroid<PFP>(map,*it,attributs);
		Dart dd = *it;
		do
		{
			map.splitFace(map.phi1(dd),map.phi_1(dd));
			map.cutEdge(map.phi1(dd));
			mCorner.mark(map.phi2(map.phi1(dd)));
			attributs[map.template phi<11>(dd)] = c*2.0/3.0f + attributs[dd]/3.0f;
			dd = map.phi1(map.phi1(map.phi1(map.phi2(map.phi1(dd)))));
		} while(!mCorner.isMarked(dd));
	}

	//third pass create center face
	for (std::vector<Dart>::iterator it = dOrig.begin(); it != dOrig.end(); ++it)
	{
		Dart dd = map.phi2(map.phi1(*it));
		do {
			mCorner.unmark(dd);
			Dart dNext = map.phi1(map.phi1(map.phi1(dd)));
			map.splitFace(dd,dNext);
			dd = dNext;
		} while(mCorner.isMarked(dd));
	}
}

529
530
template <typename PFP>
void LoopSubdivision(typename PFP::MAP& map, typename PFP::TVEC3& position, const FunctorSelect& selected)
Pierre Kraemer's avatar
Pierre Kraemer committed
531
{
532
	LoopSubdivision<PFP, typename PFP::TVEC3, typename PFP::VEC3>(map, position, selected) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
533
534
}

535
536
537
template <typename PFP>
void reverseOrientation(typename PFP::MAP& map)
{
538
	AttributeHandler<unsigned int> emb0(&map, map.getEmbeddingAttributeVector(VERTEX)) ;
539
540
	if(emb0.isValid())
	{
541
		AttributeHandler<unsigned int> new_emb0 = map.template addAttribute<unsigned int>(DART, "new_EMB_0") ;
542
543
544
545
546
547
		for(Dart d = map.begin(); d != map.end(); map.next(d))
			new_emb0[d] = emb0[map.phi1(d)] ;
		map.template swapAttributes<unsigned int>(emb0, new_emb0) ;
		map.template removeAttribute<unsigned int>(new_emb0) ;
	}

548
549
	AttributeHandler<Dart> phi1 = map.template getAttribute<Dart>(DART, "phi1") ;
	AttributeHandler<Dart> phi_1 = map.template getAttribute<Dart>(DART, "phi_1") ;
550
551
552
	map.template swapAttributes<Dart>(phi1, phi_1) ;
}

553
template <typename PFP>
Pierre Kraemer's avatar
Pierre Kraemer committed
554
void computeDual(typename PFP::MAP& map, const FunctorSelect& selected)
Pierre Kraemer's avatar
Pierre Kraemer committed
555
{
556
557
558
559
	AttributeHandler<Dart> phi1 = map.template getAttribute<Dart>(DART, "phi1") ;
	AttributeHandler<Dart> phi_1 = map.template getAttribute<Dart>(DART, "phi_1") ;
	AttributeHandler<Dart> new_phi1 = map.template addAttribute<Dart>(DART, "new_phi1") ;
	AttributeHandler<Dart> new_phi_1 = map.template addAttribute<Dart>(DART, "new_phi_1") ;
560
561
562

	for(Dart d = map.begin(); d != map.end(); map.next(d))
	{
563
		Dart dd = map.alpha_1(d) ;
564
		new_phi1[d] = dd ;
Pierre Kraemer's avatar
Pierre Kraemer committed
565
		new_phi_1[dd] = d ;
566
567
568
	}

	map.template swapAttributes<Dart>(phi1, new_phi1) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
569
	map.template swapAttributes<Dart>(phi_1, new_phi_1) ;
570
571

	map.template removeAttribute<Dart>(new_phi1) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
572
573
	map.template removeAttribute<Dart>(new_phi_1) ;

574
	map.swapEmbeddingContainers(VERTEX, FACE) ;
575
576

	reverseOrientation<PFP>(map) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
577
578
}

Pierre Kraemer's avatar
Pierre Kraemer committed
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
inline double sqrt3_K(unsigned int n)
{
	switch(n)
	{
		case 1: return 0.333333 ;
		case 2: return 0.555556 ;
		case 3: return 0.5 ;
		case 4: return 0.444444 ;
		case 5: return 0.410109 ;
		case 6: return 0.388889 ;
		case 7: return 0.375168 ;
		case 8: return 0.365877 ;
		case 9: return 0.359328 ;
		case 10: return 0.354554 ;
		case 11: return 0.350972 ;
		case 12: return 0.348219 ;
		default:
			double t = cos((2.0*M_PI)/double(n)) ;
			return (4.0 - t) / 9.0 ;
	}
Pierre Kraemer's avatar
Pierre Kraemer committed
599
600
}

601
602
template <typename PFP>
void Sqrt3Subdivision(typename PFP::MAP& map, typename PFP::TVEC3& position, const FunctorSelect& selected)
Pierre Kraemer's avatar
Pierre Kraemer committed
603
{
Pierre Kraemer's avatar
Pierre Kraemer committed
604
605
606
	typedef typename PFP::VEC3 VEC3 ;
	typedef typename PFP::REAL REAL ;

607
	AttributeHandler<VEC3> positionF = map.template getAttribute<VEC3>(FACE, "position") ;
Pierre Kraemer's avatar
Pierre Kraemer committed
608
	if(!positionF.isValid())
609
		positionF = map.template addAttribute<VEC3>(FACE, "position") ;
Pierre Kraemer's avatar
Pierre Kraemer committed
610
611
	Algo::Geometry::computeCentroidFaces<PFP>(map, position, positionF) ;

612
	computeDual<PFP>(map, selected);
Pierre Kraemer's avatar
Pierre Kraemer committed
613
614
615
616
617

	AttributeHandler<VEC3> tmp = position ;
	position = positionF ;
	positionF = tmp ;

618
	CellMarker m(map, VERTEX) ;
Pierre Kraemer's avatar
Pierre Kraemer committed
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
	m.markAll() ;

	trianguleFaces<PFP>(map, position, positionF, selected);

	for(Dart d = map.begin(); d != map.end(); map.next(d))
	{
		if(!m.isMarked(d))
		{
			m.mark(d) ;
			VEC3 P = position[d] ;
			VEC3 newP(0) ;
			unsigned int val = 0 ;
			Dart vit = d ;
			do
			{
				newP += position[map.phi2(vit)] ;
				++val ;
				vit = map.alpha1(vit) ;
			} while(vit != d) ;
			REAL K = sqrt3_K(val) ;
			newP *= REAL(3) ;
			newP -= REAL(val) * P ;
			newP *= K / REAL(2 * val) ;
			newP += (REAL(1) - K) * P ;
			position[d] = newP ;
		}
	}
Pierre Kraemer's avatar
Pierre Kraemer committed
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
}


template <typename PFP, typename EMBV, typename EMB>
void hexaCutVolume(typename PFP::MAP& map, Dart d, EMBV& attributs)
{
	DartMarker mf(map) ; //mark face
	DartMarker me(map) ; //mark edge

	DartMarkerStore mark(map);		// Lock a marker

	// compute volume centroid
	EMB volCenter = Algo::Geometry::volumeCentroidGen<PFP,EMBV,EMB>(map, d, attributs);

	//Store faces that are traversed and start with the face of d
	std::vector<Dart> visitedFaces;
	visitedFaces.reserve(100);
	visitedFaces.push_back(d);
	std::vector<Dart>::iterator face;

	//Store the edges before the cutEdge
	std::vector<Dart> oldEdges;
	oldEdges.reserve(100);
	std::vector<Dart>::iterator edge;

	//Store a dart from a each face
	std::vector<Dart> faces;
	faces.reserve(100);

	//Store the darts from quadrangulated faces
	std::vector<Dart> quadfaces;
	quadfaces.reserve(100);

	// First pass : for every face added to the list save a dart
	for (face = visitedFaces.begin(); face != visitedFaces.end(); ++face)
	{
		if (!mark.isMarked(*face))		// Face has not been visited yet
		{

			faces.push_back(*face);

			Dart dNext = *face ;
			do
			{
				mark.mark(dNext);					// Mark
				oldEdges.push_back(dNext);
				Dart adj = map.phi2(dNext);				// Get adjacent face
				if (adj != dNext && !mark.isMarked(adj))
					visitedFaces.push_back(adj);	// Add it
				dNext = map.phi1(dNext);
			} while(dNext != *face);
		}
	}

	//Second pass : cut the edges && quadrangule the faces
	for (face = faces.begin(); face != faces.end(); ++face)
	{
		Dart dNext = *face;

		//parcours de la face pour couper les aretes
		do
		{
			Dart d = dNext;
			dNext = map.phi1(dNext);

			if(!me.isMarked(d))
			{
				Dart f = map.phi1(d);
				map.cutEdge(d);
				Dart e = map.phi1(d);
				attributs[e] = attributs[d];
				attributs[e] += attributs[f];
				attributs[e] *= 0.5;

720
721
				me.markOrbit(EDGE, d);
				me.markOrbit(EDGE, e);
Pierre Kraemer's avatar
Pierre Kraemer committed
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
			}

		}while (dNext != *face);

		//quadrangulation
		EMB center = Algo::Geometry::faceCentroidGen<PFP,EMBV,EMB>(map, *face, attributs);	// compute center
		Dart cf = Algo::Modelisation::quadranguleFace<PFP>(map, *face);	// quadrangule the face
		attributs[cf] = center;	// affect the data to the central vertex

		Dart e = cf;
		do
		{
			quadfaces.push_back(e);
			quadfaces.push_back(map.phi2(e));
			e = map.phi2(map.phi_1(e));
		}while (e != cf);
	}

	DartMarker moe(map) ; //mark edge

	//Third pass : deconnect the corners
	for (edge = oldEdges.begin(); edge != oldEdges.end(); ++edge)
	{

		map.unsewFaces(map.phi1(*edge));
747
		moe.markOrbit(DART,map.phi1(*edge));
Pierre Kraemer's avatar
Pierre Kraemer committed
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
	}

	//Thourth pass : close the hole
	for (edge = oldEdges.begin(); edge != oldEdges.end(); ++edge)
	{
		if(moe.isMarked(map.phi1(*edge)))
		{
			map.closeHole(map.phi1(*edge));
			moe.unmark(map.phi1(*edge));
			moe.unmark(map.phi1(map.phi2(map.phi_1(*edge))));
			moe.unmark(map.phi1(map.phi1(map.phi2(*edge))));

			Dart cf = Algo::Modelisation::quadranguleFace<PFP>(map,  map.phi1(map.phi2(map.phi1(*edge))));	// quadrangule the face
			attributs[cf] = volCenter;	// affect the data to the central vertex
		}
	}

765
766
	moe.unmarkAll();

Pierre Kraemer's avatar
Pierre Kraemer committed
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
	//Fifth pass : traversal to phi3 sewing
	std::vector<Dart>::iterator nvol;
	for (nvol = quadfaces.begin(); nvol != quadfaces.end(); nvol = nvol + 2)
	{
		map.sewVolumes(map.phi2(*nvol), map.phi2(*(nvol+1)));
	}

//	// Compute edge points
//	for (edge = oldEdges.begin(); edge != oldEdges.end(); ++edge)
//	{
//		Dart x = *edge;
//		// other side of the edge
//		Dart y = map.phi2(x);
//		if (y != x)
//		{
//			Dart f1 = map.phi_1(x);
//			Dart f2 = map.phi2(map.phi1(y));
//			EMB temp = AttribOps::zero<EMB,PFP>();
//			temp =  attributs[f1];
//			temp += attributs[f2];			// E' = (V0+V1+F1+F2)/4
//			temp *= 0.25;
//			attributs[x] *= 0.5;
//			attributs[x] += temp;
//		}
//		// else nothing to do point already in the middle of segment
//	}
//
//	// Compute vertex points
//	for (face = faces.begin(); face != faces.end(); ++face)
//	{
//		//m0.unmark(*face);
//
//		EMB temp = AttribOps::zero<EMB,PFP>();
//		EMB temp2 = AttribOps::zero<EMB,PFP>();
//
//		unsigned int n = 0;
//		Dart x = *face;
//		do
//		{
//			Dart m = map.phi1(x);
//			Dart f = map.phi2(m);
//			Dart v = map.template phi<11>(f);
//
//			temp += attributs[f];
//			temp2 += attributs[v];
//
//			++n;
//			x = map.alpha1(x);
//		} while (x != *face);
//
//		EMB emcp = attributs[*face];
//		emcp *= double((n-2)*n);		// V' = (n-2)/n*V + 1/n2 *(F+E)
//		emcp += temp;
//		emcp += temp2;
//		emcp /= double(n*n);
//
//		attributs[*face] = emcp ;
//	}
}

827

untereiner's avatar
untereiner committed
828
829
830
831
832
833
template <typename PFP>
void splitVolumes(typename PFP::MAP& map, Dart d, typename PFP::TVEC3& position)
{
	//Cut the edges
	DartMarker me(map, EDGE);
	for(Dart d = map.begin(); d != map.end(); map.next(d))
834
	{
untereiner's avatar
untereiner committed
835
		if(!me.isMarked(d))
836
		{
untereiner's avatar
untereiner committed
837
838
839
840
841
			// Cut the edge
			Dart dd = map.phi2(d) ;
			typename PFP::VEC3 p1 = position[d] ;
			typename PFP::VEC3 p2 = position[map.phi1(d)] ;
			map.cutEdge(d) ;
842

untereiner's avatar
untereiner committed
843
			position[map.phi1(d)] = (p1 + p2) * typename PFP::REAL(0.5) ;
844

untereiner's avatar
untereiner committed
845
846
			me.markOrbit(EDGE, d);
			me.markOrbit(EDGE, map.phi1(d));
847
848
849
		}
	}

untereiner's avatar
untereiner committed
850
851
	DartMarker mf(map, FACE);
	for(Dart d = map.begin(); d != map.end(); map.next(d))
852
	{
untereiner's avatar
untereiner committed
853
		if(!mf.isMarked(d))
854
855
856
		{


untereiner's avatar
untereiner committed
857
858
			mf.markOrbit(FACE, d);
		}
859
860
	}

untereiner's avatar
untereiner committed
861
862

//	//and split the faces
863
//
untereiner's avatar
untereiner committed
864
//	//Insert the middleFaces
865
//
untereiner's avatar
untereiner committed
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
//	// first cut the edges (if they are not already)
//	Dart t = d;
//	do
//	{
//		if(!map.edgeIsSubdivided(map.phi1(map.phi2(t))))
//			Algo::IHM::subdivideEdge<PFP>(map, map.phi1(map.phi2(t)), position) ;
//		t = map.phi1(t);
//	}
//	while(t != d);
//
//	Dart neighboordVolume = map.phi1(map.phi1(map.phi2(d)));
//
//	//Split the faces and open the midlle
//	do
//	{
//		Dart t2 = map.phi2(t);
//
//		Dart face2 = map.phi1(map.phi1(t2));
//		map.splitFace(map.phi_1(t2), face2);
//		map.unsewFaces(map.phi1(map.phi1(t2)));
//
//		t = map.phi1(t);
//	}
//	while(t != d);
//
//
//	map.closeHole(map.phi1(map.phi1(map.phi2(d))));
//	map.closeHole(map.phi_1(neighboordVolume));
//	map.sewVolumes(map.phi2(map.phi1(map.phi1(map.phi2(d)))), map.phi2(map.phi_1(neighboordVolume)));
895
896
897

}

untereiner's avatar
untereiner committed
898

Pierre Kraemer's avatar
Pierre Kraemer committed
899
900
901
902
903
} // namespace Modelisation

} // namespace Algo

} // namespace CGoGN