subdivision3.hpp 25.2 KB
Newer Older
untereiner's avatar
untereiner 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           *
untereiner's avatar
untereiner 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/                                           *
untereiner's avatar
untereiner committed
21 22 23 24 25 26 27
* Contact information: cgogn@unistra.fr                                        *
*                                                                              *
*******************************************************************************/

#include "Algo/Geometry/centroid.h"
#include "Algo/Modelisation/subdivision.h"
#include "Algo/Modelisation/extrusion.h"
28
#include "Geometry/intersection.h"
29 30 31
#include "NL/nl.h"
#include "Algo/LinearSolving/basic.h"
#include "Algo/Geometry/laplacian.h"
untereiner's avatar
untereiner committed
32 33 34 35 36 37 38

namespace CGoGN
{

namespace Algo
{

39 40 41
namespace Volume
{

untereiner's avatar
untereiner committed
42 43 44 45 46 47
namespace Modelisation
{

template <typename PFP>
Dart cut3Ear(typename PFP::MAP& map, Dart d)
{
Pierre Kraemer's avatar
Pierre Kraemer committed
48 49 50 51 52 53 54 55
	Dart e = d;
	int nb = 0;
	Dart dNew;

	Dart dRing;
	Dart dRing2;

	//count the valence of the vertex
untereiner's avatar
untereiner committed
56 57
	do
	{
Pierre Kraemer's avatar
Pierre Kraemer committed
58 59 60
		nb++;
		e = map.phi1(map.phi2(e));
	} while (e != d);
untereiner's avatar
untereiner committed
61

Pierre Kraemer's avatar
Pierre Kraemer committed
62 63 64 65 66 67 68
	if(nb < 3)
	{
		CGoGNout << "Warning : cannot cut 2 volumes without creating a degenerated face " << CGoGNendl;
		return d;
	}
	else
	{
69 70
		std::vector<Dart> vPath;

Pierre Kraemer's avatar
Pierre Kraemer committed
71 72 73 74 75 76
		//triangulate around the vertex
		do
		{
			Dart dN = map.phi1(map.phi2(e));
			if(map.template phi<111>(e) != e)
				map.splitFace(map.phi_1(e), map.phi1(e));
untereiner's avatar
untereiner committed
77

78
			dRing = map.phi2(map.phi1(e));
untereiner's avatar
untereiner committed
79

80
			vPath.push_back(dRing); //remember all darts from the ring
untereiner's avatar
untereiner committed
81

Pierre Kraemer's avatar
Pierre Kraemer committed
82 83
			e = dN;
		} while (e != d);
untereiner's avatar
untereiner committed
84

85
		map.splitVolume(vPath);
Pierre Kraemer's avatar
Pierre Kraemer committed
86
	}
untereiner's avatar
untereiner committed
87

Pierre Kraemer's avatar
Pierre Kraemer committed
88 89
	return map.phi2(dRing);
}
untereiner's avatar
untereiner committed
90

91
template <typename PFP>
Thomas's avatar
Thomas committed
92
Dart sliceConvexVolume(typename PFP::MAP& map, VertexAttribute<typename PFP::VEC3>& position, Dart d, Geom::Plane3D<typename PFP::REAL > pl)
93 94 95 96
{
	Dart dRes=NIL;
	unsigned int nbInter = 0;
	unsigned int nbVertices = 0;
Thomas's avatar
Thomas committed
97 98
	CellMarkerStore<VERTEX> vs(map);			//marker for new vertices from edge cut
	CellMarkerStore<FACE> cf(map);
99 100
	Dart dPath;

Thomas's avatar
Thomas committed
101 102
	MarkerForTraversor<typename PFP::MAP::ParentMap, EDGE > mte(map);
	MarkerForTraversor<typename PFP::MAP::ParentMap, FACE > mtf(map);
103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118

	//search edges and vertices crossing the plane
	Traversor3WE<typename PFP::MAP::ParentMap > te(map,d);
	for(Dart dd = te.begin() ;dd != te.end() ; dd = te.next())
	{
		if(!mte.isMarked(dd))
		{
			if(fabs(pl.distance(position[dd]))<0.000001f)
			{
				nbVertices++;
				vs.mark(dd); //mark vertex on slicing path
				mte.mark(dd);
			}
			else
			{
				typename PFP::VEC3 interP;
119
				typename PFP::VEC3 vec(Surface::Geometry::vectorOutOfDart<PFP>(map,dd,position));
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
				Geom::Intersection inter = Geom::intersectionLinePlane<typename PFP::VEC3, typename Geom::Plane3D<typename PFP::REAL > >(position[dd],vec,pl,interP);

				if(inter==Geom::FACE_INTERSECTION)
				{
					Dart dOp = map.phi1(dd);
					typename PFP::VEC3 v2(interP-position[dd]);
					typename PFP::VEC3 v3(interP-position[dOp]);
					if(vec.norm2()>v2.norm2() && vec.norm2()>v3.norm2())
					{
						nbInter++;

						cf.mark(dd);			//mark face and opposite face to split
						cf.mark(map.phi2(dd));

						map.cutEdge(dd);
						Dart dN = map.phi1(dd);

						mte.mark(dN);

						vs.mark(dN);			//mark vertex for split
						position[dN] = interP; 	//place
					}
				}
			}
		}
	}

//	std::cout << "edges cut: " << nbInter << std::endl;
	unsigned int nbSplit=0;

	//slice when at least two edges are concerned
	if(nbInter>1)
	{
		Traversor3WF<typename PFP::MAP::ParentMap > tf(map,d);
		for(Dart dd = tf.begin() ; dd != tf.end() ; dd = tf.next())
		{
			//for faces with a new vertex
			if(cf.isMarked(dd))
			{
				cf.unmark(dd);

				Dart dS = dd;
				bool split=false;

164 165
				do
				{
166 167 168 169 170
					//find the new vertex
					if(vs.isMarked(dS))
					{
						Dart dSS = map.phi1(dS);
						//search an other new vertex (or an existing vertex intersected with the plane) in order to split the face
171 172
						do
						{
173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214
							if(vs.isMarked(dSS))
							{
								nbSplit++;
								map.splitFace(dS,dSS);
								dPath=map.phi_1(dS);
								split=true;
							}
							dSS = map.phi1(dSS);
						} while(!split && dSS!=dS);
					}
					dS = map.phi1(dS);
				} while(!split && dS!=dd);
			}
		}

//		std::cout << "face split " << nbSplit << std::endl;

		//define the path to split
		std::vector<Dart> vPath;
		vPath.reserve((nbSplit+nbVertices)+1);
		vPath.push_back(dPath);
		for(std::vector<Dart>::iterator it = vPath.begin() ;it != vPath.end() ; ++it)
		{
			Dart dd = map.phi1(*it);

			Dart ddd = map.phi1(map.phi2(dd));

			while(!vs.isMarked(map.phi1(ddd)) && ddd!=dd)
				ddd = map.phi1(map.phi2(ddd));

			if(vs.isMarked(map.phi1(ddd)) && !map.sameVertex(ddd,*vPath.begin()))
				vPath.push_back(ddd);
		}

		assert(vPath.size()>2);
		map.splitVolume(vPath);
		dRes = map.phi2(*vPath.begin());
	}

	return dRes;
}

215
template <typename PFP>
Thomas's avatar
Thomas committed
216
Dart sliceConvexVolume(typename PFP::MAP& map, VertexAttribute<typename PFP::VEC3>& position, Dart d, CellMarker<EDGE>& edgesToCut, CellMarker<VERTEX>& verticesToSplit)
217 218 219 220 221 222
{
	typedef typename PFP::VEC3 VEC3;

	Dart dRes;
	unsigned int nbInter = 0;
	unsigned int nbVertices = 0;
Thomas's avatar
Thomas committed
223 224
	CellMarkerStore<VERTEX> vs(map);			//marker for new vertices from edge cut
	CellMarkerStore<FACE> cf(map);
225 226
	Dart dPath;

Thomas's avatar
Thomas committed
227 228
	MarkerForTraversor<typename PFP::MAP::ParentMap, EDGE > mte(map);
	MarkerForTraversor<typename PFP::MAP::ParentMap, FACE > mtf(map);
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 306 307 308 309 310 311 312 313 314

	//search edges and vertices crossing the plane
	Traversor3WE<typename PFP::MAP::ParentMap > te(map,d);
	for(Dart dd = te.begin() ;dd != te.end() ; dd = te.next())
	{
		if(!mte.isMarked(dd) && edgesToCut.isMarked(dd))
		{
			nbInter++;
			VEC3 p = (position[dd]+position[map.phi1(dd)])*0.5f;
			cf.mark(dd);			//mark face and opposite face to split
			cf.mark(map.phi2(dd));

			map.cutEdge(dd);
			Dart dN = map.phi1(dd);

			mte.mark(dN);

			vs.mark(dN);		//mark vertex for split
			position[dN] = p;
		}
	}

//	std::cout << "edges cut: " << nbInter << std::endl;
	unsigned int nbSplit=0;

	//at least two edges are concerned
	assert(nbInter>1);

	Traversor3WF<typename PFP::MAP::ParentMap > tf(map,d);
	for(Dart dd = tf.begin() ; dd != tf.end() ; dd = tf.next())
	{
		//for faces with a new vertex
		if(cf.isMarked(dd))
		{
			cf.unmark(dd);

			Dart dS = dd;
			bool split=false;

			do {
				//find the new vertex
				if(vs.isMarked(dS) || verticesToSplit.isMarked(dS))
				{
					Dart dSS = map.phi1(dS);
					//search an other new vertex (or an existing vertex intersected with the plane) in order to split the face
					do {
						if(vs.isMarked(dSS) || verticesToSplit.isMarked(dSS))
						{
							nbSplit++;
							map.splitFace(dS,dSS);
							dPath=map.phi_1(dS);
							split=true;
						}
						dSS = map.phi1(dSS);
					} while(!split && dSS!=dS);
				}
				dS = map.phi1(dS);
			} while(!split && dS!=dd);
		}

		//define the path to split
		std::vector<Dart> vPath;
		vPath.reserve((nbSplit+nbVertices)+1);
		vPath.push_back(dPath);
		for(std::vector<Dart>::iterator it = vPath.begin() ;it != vPath.end() ; ++it)
		{
			Dart dd = map.phi1(*it);

			Dart ddd = map.phi1(map.phi2(dd));

			while(!vs.isMarked(map.phi1(ddd)) && ddd!=dd)
				ddd = map.phi1(map.phi2(ddd));

			if(vs.isMarked(map.phi1(ddd)) && !map.sameVertex(ddd,*vPath.begin()))
				vPath.push_back(ddd);
		}

		assert(vPath.size()>2);
		map.splitVolume(vPath);
		dRes = map.phi2(*vPath.begin());
	}

	return dRes;
}

template <typename PFP>
Thomas's avatar
Thomas committed
315
std::vector<Dart> sliceConvexVolumes(typename PFP::MAP& map, VertexAttribute<typename PFP::VEC3>& position,CellMarker<VOLUME>& volumesToCut, CellMarker<EDGE>& edgesToCut, CellMarker<VERTEX>& verticesToSplit)
316
{
Thomas's avatar
Thomas committed
317
    std::vector<Dart> vRes;
318 319

    typedef typename PFP::VEC3 VEC3;
Thomas's avatar
Thomas committed
320
    CellMarker<VERTEX> localVerticesToSplit(map); //marker for new vertices from edge cut
321 322 323

    //Step 1: Cut the edges and mark the resulting vertices as vertices to be face-split
    TraversorE<typename PFP::MAP> te(map);
Thomas's avatar
Thomas committed
324
    CellMarkerStore<FACE> cf(map);
325 326 327 328 329 330

    for(Dart d = te.begin(); d != te.end(); d=te.next()) //cut all edges
    {
        if(edgesToCut.isMarked(d))
        {
            VEC3 p = (position[d]+position[map.phi1(d)])*0.5f;
331 332 333 334 335

            //turn around the edge and mark for future split face
            Traversor3EF<typename PFP::MAP> t3ef(map,d);
            for(Dart dd = t3ef.begin() ; dd != t3ef.end() ; dd = t3ef.next())
            	cf.mark(dd);			//mark face to split
336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356

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

            localVerticesToSplit.mark(dN);		//mark vertex for split
            position[dN] = p;
        }
    }

    //Step 2: Split faces with cut edges
    TraversorF<typename PFP::MAP> tf(map);
    for(Dart d = tf.begin(); d != tf.end(); d=tf.next())
    {
        if(cf.isMarked(d))
        {
            cf.unmark(d);
            Dart dS = d;
            bool split=false;
            do
            {
                //find the new vertex
357
                if(localVerticesToSplit.isMarked(dS) || verticesToSplit.isMarked(dS))
358
                {
359 360
                	//start from phi1(phi1()) to avoid the creation of faces of degree 2
                    Dart dSS = map.phi1(map.phi1(dS));
361
                    //search an other new vertex (or an existing vertex to split) in order to split the face
362

363 364
                    do
                    {
365 366
                        if((localVerticesToSplit.isMarked(dSS) || verticesToSplit.isMarked(dSS))
                        		&& !map.sameVertex(dS,dSS))
367 368 369 370 371 372
                        {
                            map.splitFace(dS,dSS);
                            split=true;
                        }
                        dSS = map.phi1(dSS);
                    } while(!split && dSS!=dS);
373
                    split=true; //go out of the first loop if no split case has been found
374 375 376 377 378 379 380 381 382 383
                }
                dS = map.phi1(dS);
            } while(!split && dS!=d);
        }
    }

    //Step 3 : Find path and split volumes
    TraversorW<typename PFP::MAP> tw(map);
    for(Dart d = tw.begin(); d != tw.end(); d=tw.next()) //Parcours des volumes
    {
384 385
        if(volumesToCut.isMarked(d))
        {
386 387 388 389 390 391 392 393 394
            Traversor3WV<typename PFP::MAP> t3wv(map,d);
            Dart dPath;
            bool found=false;

            //find a vertex of the volume to start the path to split
            for(Dart dd = t3wv.begin(); dd != t3wv.end() && !found; dd=t3wv.next())
            {
                if(localVerticesToSplit.isMarked(dd) || verticesToSplit.isMarked(dd))
                {
395
                    Dart ddd = dd;
396
                    while(!localVerticesToSplit.isMarked(map.phi1(ddd))
397
                    		&& !verticesToSplit.isMarked(map.phi1(ddd)))
398 399 400 401 402 403 404
                        ddd = map.phi1(map.phi2(ddd));
                    found=true;
                    dPath=ddd;
                }
            }
            //define the path to split
            std::vector<Dart> vPath;
405
            vPath.reserve(32);
406
            vPath.push_back(dPath);
Thomas's avatar
Thomas committed
407
            CellMarker<FACE> cmf(map);
408

409

410
            //define the path to split for the whole volume
411 412
            bool pathFound=false;
            for(std::vector<Dart>::iterator it = vPath.begin() ; !pathFound && it != vPath.end(); ++it)
413 414 415
            {
                Dart dd = map.phi1(*it);

416 417 418
                if(map.sameVertex(dd,*vPath.begin()))
                	pathFound=true;
                else
419
                {
420 421 422 423 424 425
                	Dart ddd = map.phi1(map.phi2(dd));

                	while(!localVerticesToSplit.isMarked(map.phi1(ddd)) && !verticesToSplit.isMarked(map.phi1(ddd)))
                		ddd = map.phi1(map.phi2(ddd));

                	vPath.push_back(ddd);
426 427 428 429 430
                }
            }

            assert(vPath.size()>2);
            map.splitVolume(vPath);
431
            vRes.push_back(map.phi2(*vPath.begin()));
432 433 434
        }
    }

435
    return vRes;
436 437
}

untereiner's avatar
untereiner committed
438
template <typename PFP, typename EMBV, typename EMB>
439
void catmullClarkVol(typename PFP::MAP& map, EMBV& attributs)
untereiner's avatar
untereiner committed
440
{
441
	//std::vector<Dart> l_centers;
untereiner's avatar
untereiner committed
442 443 444
	std::vector<Dart> l_vertices;

	//pre-computation : compute the centroid of all volume
445
	VolumeAutoAttribute<EMB> attBary(map);
446
	Volume::Geometry::computeCentroidVolumes<PFP>(map, const_cast<const EMBV&>(attributs), attBary);
untereiner's avatar
untereiner committed
447

448 449 450 451 452
	//subdivision
	//1. cut edges
	DartMarkerNoUnmark mv(map);
	TraversorE<typename PFP::MAP> travE(map);
	for (Dart d = travE.begin(); d != travE.end(); d = travE.next())
untereiner's avatar
untereiner committed
453 454
	{
		//memorize each vertices per volumes
455
		if( !mv.isMarked(d))
untereiner's avatar
untereiner committed
456 457
		{
			l_vertices.push_back(d);
Pierre Kraemer's avatar
Pierre Kraemer committed
458
			mv.markOrbit<PFP::MAP::VERTEX_OF_PARENT>(d);
untereiner's avatar
untereiner committed
459 460
		}

461 462 463
		Dart f = map.phi1(d);
		map.cutEdge(d) ;
		Dart e = map.phi1(d) ;
untereiner's avatar
untereiner committed
464

465 466 467
		attributs[e] =  attributs[d];
		attributs[e] += attributs[f];
		attributs[e] *= 0.5;
untereiner's avatar
untereiner committed
468

untereiner's avatar
untereiner committed
469 470
		travE.skip(d) ;
		travE.skip(e) ;
untereiner's avatar
untereiner committed
471 472
	}

473 474 475
	//2. split faces - quadrangule faces
	TraversorF<typename PFP::MAP> travF(map) ;
	for (Dart d = travF.begin(); d != travF.end(); d = travF.next())
untereiner's avatar
untereiner committed
476
	{
477
		EMB center = Surface::Geometry::faceCentroidGen<PFP,EMBV,EMB>(map,d,attributs);
untereiner's avatar
untereiner committed
478

479 480 481
		Dart dd = map.phi1(d) ;
		Dart next = map.phi1(map.phi1(dd)) ;
		map.splitFace(dd, next) ;
untereiner's avatar
untereiner committed
482

483 484
		Dart ne = map.phi2(map.phi_1(dd)) ;
		map.cutEdge(ne) ;
untereiner's avatar
untereiner committed
485
		travF.skip(dd) ;
untereiner's avatar
untereiner committed
486

487
		attributs[map.phi1(ne)] = center;
untereiner's avatar
untereiner committed
488

489 490 491 492 493
		dd = map.phi1(map.phi1(next)) ;
		while(dd != ne)
		{
			Dart tmp = map.phi1(ne) ;
			map.splitFace(tmp, dd) ;
untereiner's avatar
untereiner committed
494
			travF.skip(tmp) ;
495
			dd = map.phi1(map.phi1(dd)) ;
untereiner's avatar
untereiner committed
496
		}
497

untereiner's avatar
untereiner committed
498
		travF.skip(ne) ;
untereiner's avatar
untereiner committed
499 500
	}

501 502 503 504
	//3. create inside volumes

	std::vector<std::pair<Dart, Dart> > subdividedFaces;
	subdividedFaces.reserve(2048);
untereiner's avatar
untereiner committed
505 506
	for (std::vector<Dart>::iterator it = l_vertices.begin(); it != l_vertices.end(); ++it)
	{
507 508 509
		Dart e = *it;
		std::vector<Dart> v ;

untereiner's avatar
untereiner committed
510 511
		do
		{
512 513 514 515 516
			v.push_back(map.phi1(e));
			v.push_back(map.phi1(map.phi1(e)));

			if(!map.PFP::MAP::ParentMap::isBoundaryEdge(map.phi1(e)))
				subdividedFaces.push_back(std::pair<Dart,Dart>(map.phi1(e),map.phi2(map.phi1(e))));
untereiner's avatar
untereiner committed
517

518 519
			if(!map.PFP::MAP::ParentMap::isBoundaryEdge(map.phi1(map.phi1(e))))
				subdividedFaces.push_back(std::pair<Dart,Dart>(map.phi1(map.phi1(e)),map.phi2(map.phi1(map.phi1(e)))));
untereiner's avatar
untereiner committed
520

521
			e = map.phi2(map.phi_1(e));
untereiner's avatar
untereiner committed
522
		}
523
		while(e != *it);
untereiner's avatar
untereiner committed
524

525 526 527
		//
		// SplitSurfaceInVolume
		//
untereiner's avatar
untereiner committed
528

529 530
		std::vector<Dart> vd2 ;
		vd2.reserve(v.size());
untereiner's avatar
untereiner committed
531

532 533 534 535 536
		// save the edge neighbors darts
		for(std::vector<Dart>::iterator it2 = v.begin() ; it2 != v.end() ; ++it2)
		{
			vd2.push_back(map.phi2(*it2));
		}
untereiner's avatar
untereiner committed
537

538
		assert(vd2.size() == v.size());
untereiner's avatar
untereiner committed
539

540
		map.PFP::MAP::ParentMap::splitSurface(v, true, false);
untereiner's avatar
untereiner committed
541

542 543 544 545 546
		// follow the edge path a second time to embed the vertex, edge and volume orbits
		for(unsigned int i = 0; i < v.size(); ++i)
		{
			Dart dit = v[i];
			Dart dit2 = vd2[i];
untereiner's avatar
untereiner committed
547

548
			// embed the vertex embedded from the origin volume to the new darts
Pierre Kraemer's avatar
Pierre Kraemer committed
549
			if(map.template isOrbitEmbedded<VERTEX>())
550
			{
Pierre Kraemer's avatar
Pierre Kraemer committed
551 552
				map.template copyDartEmbedding<VERTEX>(map.phi2(dit), map.phi1(dit));
				map.template copyDartEmbedding<VERTEX>(map.phi2(dit2), map.phi1(dit2));
553
			}
untereiner's avatar
untereiner committed
554

555
			// embed the edge embedded from the origin volume to the new darts
Pierre Kraemer's avatar
Pierre Kraemer committed
556
			if(map.template isOrbitEmbedded<EDGE>())
557
			{
Pierre Kraemer's avatar
Pierre Kraemer committed
558 559 560
				unsigned int eEmb = map.template getEmbedding<EDGE>(dit) ;
				map.template setDartEmbedding<EDGE>(map.phi2(dit), eEmb);
				map.template setDartEmbedding<EDGE>(map.phi2(dit2), eEmb);
561
			}
untereiner's avatar
untereiner committed
562

563
			// embed the volume embedded from the origin volume to the new darts
Pierre Kraemer's avatar
Pierre Kraemer committed
564
			if(map.template isOrbitEmbedded<VOLUME>())
565
			{
Pierre Kraemer's avatar
Pierre Kraemer committed
566 567
				map.template copyDartEmbedding<VOLUME>(map.phi2(dit), dit);
				map.template copyDartEmbedding<VOLUME>(map.phi2(dit2), dit2);
568 569
			}
		}
untereiner's avatar
untereiner committed
570

571 572 573
		//
		//
		//
untereiner's avatar
untereiner committed
574

575 576 577
		Dart dd = map.phi2(map.phi1(*it));
		Dart next = map.phi1(map.phi1(dd)) ;
		map.PFP::MAP::ParentMap::splitFace(dd, next);
untereiner's avatar
untereiner committed
578

Pierre Kraemer's avatar
Pierre Kraemer committed
579
		if (map.template isOrbitEmbedded<VERTEX>())
580
		{
Pierre Kraemer's avatar
Pierre Kraemer committed
581 582
			map.template copyDartEmbedding<VERTEX>(map.phi_1(next), dd) ;
			map.template copyDartEmbedding<VERTEX>(map.phi_1(dd), next) ;
583 584 585 586
		}

		Dart ne = map.phi2(map.phi_1(dd));
		map.PFP::MAP::ParentMap::cutEdge(ne);
untereiner's avatar
untereiner committed
587

588 589
//		dd = map.phi1(map.phi1(next)) ;
//		while(dd != ne)
untereiner's avatar
untereiner committed
590
//		{
591 592 593
//			Dart tmp = map.phi1(ne) ;
//			map.PFP::MAP::ParentMap::splitFace(tmp, dd);
//
594
//			if (map.isOrbitEmbedded<VERTEX>())
595
//			{
Pierre Kraemer's avatar
Pierre Kraemer committed
596 597
//				map.copyDartEmbedding<VERTEX>(map.phi_1(dd), tmp) ;
//				map.copyDartEmbedding<VERTEX>(map.phi_1(tmp), dd) ;
598 599 600
//			}
//
//			dd = map.phi1(map.phi1(dd)) ;
untereiner's avatar
untereiner committed
601
//		}
602 603
//
	}
untereiner's avatar
untereiner committed
604

605 606 607
//		setCurrentLevel(getMaxLevel()) ;
//		//4 couture des relations precedemment sauvegarde
//		for (std::vector<std::pair<Dart,Dart> >::iterator it = subdividedFaces.begin(); it != subdividedFaces.end(); ++it)
untereiner's avatar
untereiner committed
608
//		{
609 610 611 612 613 614
//			Dart f1 = phi2((*it).first);
//			Dart f2 = phi2((*it).second);
//
//			//if(isBoundaryFace(f1) && isBoundaryFace(f2))
//			if(phi3(f1) == f1 && phi3(f2) == f2)
//				sewVolumes(f1, f2, false);
untereiner's avatar
untereiner committed
615
//		}
616
//		setOrbitEmbedding<VERTEX>(centralDart, getEmbedding<VERTEX>(centralDart));
617 618 619 620 621 622 623 624 625 626
		//attributs[map.phi1(ne)] = attBary[*it];
//
//		setCurrentLevel(getMaxLevel() - 1) ;
//	}
//
//	//A optimiser
//
//	TraversorE<typename PFP::MAP> travE2(map);
//	for (Dart d = travE2.begin(); d != travE2.end(); d = travE2.next())
//	{
627
//		map.setOrbitEmbedding<VERTEX>(map.phi1(d), map.getEmbedding<VERTEX>(map.phi1(d)));
628 629 630 631 632
//	}
//
//	TraversorF<typename PFP::MAP> travF2(map) ;
//	for (Dart d = travF2.begin(); d != travF2.end(); d = travF2.next())
//	{
633
//		map.setOrbitEmbedding<VERTEX>(map.phi2(map.phi1(d)), map.getEmbedding<VERTEX>(map.phi2(map.phi1(d))));
634
//	}
635 636
}

637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657
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 ;
	}
}
658 659

template <typename PFP>
660
void sqrt3Vol(typename PFP::MAP& map, VertexAttribute<typename PFP::VEC3>& position)
661 662
{
	DartMarkerStore m(map);
663

664 665
	DartMarkerStore newBoundaryV(map);

666 667 668
	//
	// 1-4 flip of all tetrahedra
	//
669
	TraversorW<typename PFP::MAP> tW(map);
670 671
	for(Dart dit = tW.begin() ; dit != tW.end() ; dit = tW.next())
	{
672 673 674 675 676 677
		Traversor3WF<typename PFP::MAP> tWF(map, dit);
		for(Dart ditWF = tWF.begin() ; ditWF != tWF.end() ; ditWF = tWF.next())
		{
			if(!map.isBoundaryFace(ditWF))
				m.markOrbit<FACE>(ditWF);
		}
678

679 680 681 682 683 684 685
		typename PFP::VEC3 volCenter(0.0);
		volCenter += position[dit];
		volCenter += position[map.phi1(dit)];
		volCenter += position[map.phi_1(dit)];
		volCenter += position[map.phi_1(map.phi2(dit))];
		volCenter /= 4;

686
		Dart dres = Volume::Modelisation::Tetrahedralization::flip1To4<PFP>(map, dit);
687
		position[dres] = volCenter;
688
	}
689

690 691 692
	//
	// 2-3 swap of all old interior faces
	//
693
	TraversorF<typename PFP::MAP> tF(map);
694 695 696 697 698
	for(Dart dit = tF.begin() ; dit != tF.end() ; dit = tF.next())
	{
		if(m.isMarked(dit))
		{
			m.unmarkOrbit<FACE>(dit);
699
			Volume::Modelisation::Tetrahedralization::swap2To3<PFP>(map, dit);
700 701 702 703 704 705
		}
	}

	//
	// 1-3 flip of all boundary tetrahedra
	//
706
	TraversorW<typename PFP::MAP> tWb(map);
707 708 709 710
	for(Dart dit = tWb.begin() ; dit != tWb.end() ; dit = tWb.next())
	{
		if(map.isBoundaryVolume(dit))
		{
711 712 713 714 715 716 717
			Traversor3WE<typename PFP::MAP> tWE(map, dit);
			for(Dart ditWE = tWE.begin() ; ditWE != tWE.end() ; ditWE = tWE.next())
			{
				if(map.isBoundaryEdge(ditWE))
					m.markOrbit<EDGE>(ditWE);
			}

718 719 720 721 722 723
			typename PFP::VEC3 faceCenter(0.0);
			faceCenter += position[dit];
			faceCenter += position[map.phi1(dit)];
			faceCenter += position[map.phi_1(dit)];
			faceCenter /= 3;

724
			Dart dres = Volume::Modelisation::Tetrahedralization::flip1To3<PFP>(map, dit);
725
			position[dres] = faceCenter;
726 727

			newBoundaryV.markOrbit<VERTEX>(dres);
728 729 730 731 732 733
		}
	}

	//
	// edge-removal on all old boundary edges
	//
734
	TraversorE<typename PFP::MAP> tE(map);
735 736 737 738 739 740
	for(Dart dit = tE.begin() ; dit != tE.end() ; dit = tE.next())
	{
		if(m.isMarked(dit))
		{
			m.unmarkOrbit<EDGE>(dit);
			Dart d = map.phi2(map.phi3(map.findBoundaryFaceOfEdge(dit)));
741
			Volume::Modelisation::Tetrahedralization::swapGen3To2<PFP>(map, d);
742 743 744

		}
	}
745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 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 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 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 875 876 877 878 879 880 881 882 883 884 885

	TraversorV<typename PFP::MAP> tVg(map,selected);
	for(Dart dit = tVg.begin() ; dit != tVg.end() ; dit = tVg.next())
	{
		if(map.isBoundaryVertex(dit) && !newBoundaryV.isMarked(dit))
		{
			typename PFP::VEC3 P = position[dit] ;
			typename PFP::VEC3 newP(0) ;
			unsigned int val = 0 ;
			Dart vit = dit ;
			do
			{
				newP += position[map.phi_1(map.phi2(map.phi1(vit)))] ;
				++val ;
				vit = map.phi2(map.phi_1(vit)) ;
			} while(vit != dit) ;
			typename PFP::REAL K = sqrt3_K(val) ;
			newP *= typename PFP::REAL(3) ;
			newP -= typename PFP::REAL(val) * P ;
			newP *= K / typename PFP::REAL(2 * val) ;
			newP += (typename PFP::REAL(1) - K) * P ;
			position[dit] = newP ;
		}
	}

	//AutoVertexAttribute laplacian qui est une copie de position

//	TraversorV<typename PFP::MAP> tVg2(map,selected);
//	for(Dart dit = tVg2.begin() ; dit != tVg2.end() ; dit = tVg2.next())
//	{
//		if(!map.isBoundaryVertex(dit))
//		{
//			//modifie laplacian
//		}
//	}

	//echange lapaclian et position


	VertexAutoAttribute<typename PFP::VEC3> diffCoord(map, "diffCoord");
	Algo::Volume::Geometry::computeLaplacianTopoVertices<PFP>(map, position, diffCoord) ;

	VertexAutoAttribute<unsigned int> vIndex(map, "vIndex");

	unsigned int nb_vertices = map.template computeIndexCells<VERTEX>(vIndex);


	CellMarker<VERTEX> lockingMarker(map);

	TraversorV<typename PFP::MAP> tv(map);
	for(Dart dit = tv.begin() ; dit != tv.end() ; dit = tv.next())
	{
		if(!lockingMarker.isMarked(dit) && map.isBoundaryVertex(dit))
			lockingMarker.mark(dit);
	}


	NLContext nlContext = nlNewContext();
	nlSolverParameteri(NL_NB_VARIABLES, nb_vertices);
	nlSolverParameteri(NL_LEAST_SQUARES, NL_TRUE);
	nlSolverParameteri(NL_SOLVER, NL_CHOLMOD_EXT);


	nlMakeCurrent(nlContext);
	if(nlGetCurrentState() == NL_STATE_INITIAL)
		nlBegin(NL_SYSTEM) ;

	for(int coord = 0; coord < 3; ++coord)
	{
		LinearSolving::setupVariables<PFP>(map, vIndex, lockingMarker, position, coord) ;
		nlBegin(NL_MATRIX) ;
		LinearSolving::addRowsRHS_Laplacian_Topo<PFP>(map, vIndex, diffCoord, coord) ;
//		LinearSolving::addRowsRHS_Laplacian_Cotan<PFP>(*map, perMap->vIndex, perMap->edgeWeight, perMap->vertexArea, perMap->diffCoord, coord) ;
		nlEnd(NL_MATRIX) ;
		nlEnd(NL_SYSTEM) ;
		nlSolve() ;
		LinearSolving::getResult<PFP>(map, vIndex, position, coord) ;
		nlReset(NL_TRUE) ;
	}

	nlDeleteContext(nlContext);


//
//	float weight = 1.0;

//	LinearSolving::LinearSolver<PFP::REAL> ls(nbV);
//	ls.set_least_squares(true);

//	for(unsigned int coord = 0 ; coord < 3 ; ++coord)
//	{
//		std::cout << "coord " << coord << std::flush;

//		TraversorV<PFP::MAP> tv(map);
//		for(Dart dit = tv.begin() ; dit != tv.end() ; dit = tv.next())
//		{
//			ls.variable(indexV[dit]).set_value(position[dit][coord]);

//			if(map.isBoundaryVertex(dit))
//				ls.variable(indexV[dit]).lock();
//		}
//		std::cout << "... variables set... " << std::flush;

//		ls.begin_system();

//		TraversorV<PFP::MAP> tv2(map);
//		for(Dart dit = tv2.begin() ; dit != tv2.end() ; dit = tv2.next())
//		{
//			if(!map.isBoundaryVertex(dit))
//			{
//				ls.begin_row();

//				float sum = 0;
//				Traversor3VVaE<PFP::MAP> tvvae(map, dit);
//				for(Dart ditvvae = tvvae.begin() ; ditvvae != tvvae.end() ; ditvvae = tvvae.next())
//				{
//					ls.add_coefficient(indexV[ditvvae],weight);
//					sum += weight;
//				}

//				ls.add_coefficient(indexV[dit],-sum);
//				ls.normalize_row();

//				ls.end_row();
//			}
//		}

//		ls.end_system();
//		std::cout << "... system built... " << std::flush;
//		ls.solve();
//		std::cout << "... system solved... " << std::flush;

//		TraversorV<PFP::MAP> tv3(map);
//		for(Dart dit = tv3.begin() ; dit != tv3.end() ; dit = tv3.next())
//		{
//			position[dit][coord] = ls.variable(indexV[dit]).value();
//		}

//		ls.reset(false);
//		std::cout << "... done" << std::endl;
//	}
untereiner's avatar
untereiner committed
886 887
}

888 889 890 891 892 893 894 895 896 897 898 899 900 901 902
template <typename PFP>
void computeDual(typename PFP::MAP& map, VertexAttribute<typename PFP::VEC3>& position)
{
	// VolumeAttribute -> after dual new VertexAttribute
	VolumeAttribute<typename PFP::VEC3> positionV  = map.template getAttribute<typename PFP::VEC3, VOLUME>("position") ;
	if(!positionV.isValid())
		positionV = map.template addAttribute<typename PFP::VEC3, VOLUME>("position") ;

	// Compute Centroid for the volumes
	Algo::Volume::Geometry::computeCentroidVolumes<PFP>(map, position, positionV) ;

	// Compute the Dual mesh
	map.computeDual();
	position = positionV ;
}
untereiner's avatar
untereiner committed
903

904

905
} // namespace Modelisation
untereiner's avatar
untereiner committed
906

907
} // namespace volume
908

909
} // namespace Algo
untereiner's avatar
untereiner committed
910

911
} // namespace CGoGN
untereiner's avatar
untereiner committed
912