subdivision3.hpp 27.5 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
#include "NL/nl.h"
30
//#include "Algo/LinearSolving/basic.h"
31
#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
namespace Modelisation
{

45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
template <typename PFP>
bool isHexahedron(typename PFP::MAP& the_map, Dart d, unsigned int thread)
{
    unsigned int nbFaces = 0;

    //Test the number of faces end its valency
    Traversor3WF<typename PFP::MAP> travWF(the_map, d, false, thread);
    for(Dart dit = travWF.begin() ; dit != travWF.end(); dit = travWF.next())
    {
        //increase the number of faces
        nbFaces++;
        if(nbFaces > 6)	//too much faces
            return false;

        //test the valency of this face
        if(the_map.faceDegree(dit) != 4)
            return false;
    }

    return true;
}

untereiner's avatar
untereiner committed
67 68 69
template <typename PFP>
Dart cut3Ear(typename PFP::MAP& map, Dart d)
{
Pierre Kraemer's avatar
Pierre Kraemer committed
70 71 72 73 74 75 76 77
	Dart e = d;
	int nb = 0;
	Dart dNew;

	Dart dRing;
	Dart dRing2;

	//count the valence of the vertex
untereiner's avatar
untereiner committed
78 79
	do
	{
Pierre Kraemer's avatar
Pierre Kraemer committed
80 81 82
		nb++;
		e = map.phi1(map.phi2(e));
	} while (e != d);
untereiner's avatar
untereiner committed
83

Pierre Kraemer's avatar
Pierre Kraemer committed
84 85 86 87 88 89 90
	if(nb < 3)
	{
		CGoGNout << "Warning : cannot cut 2 volumes without creating a degenerated face " << CGoGNendl;
		return d;
	}
	else
	{
91 92
		std::vector<Dart> vPath;

Pierre Kraemer's avatar
Pierre Kraemer committed
93 94 95 96 97 98
		//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
99

100
			dRing = map.phi2(map.phi1(e));
untereiner's avatar
untereiner committed
101

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

Pierre Kraemer's avatar
Pierre Kraemer committed
104 105
			e = dN;
		} while (e != d);
untereiner's avatar
untereiner committed
106

107
		map.splitVolume(vPath);
Pierre Kraemer's avatar
Pierre Kraemer committed
108
	}
untereiner's avatar
untereiner committed
109

Pierre Kraemer's avatar
Pierre Kraemer committed
110 111
	return map.phi2(dRing);
}
untereiner's avatar
untereiner committed
112

113
template <typename PFP>
Thomas's avatar
Thomas committed
114
Dart sliceConvexVolume(typename PFP::MAP& map, VertexAttribute<typename PFP::VEC3>& position, Dart d, Geom::Plane3D<typename PFP::REAL > pl)
115 116 117 118
{
	Dart dRes=NIL;
	unsigned int nbInter = 0;
	unsigned int nbVertices = 0;
Thomas's avatar
Thomas committed
119 120
	CellMarkerStore<VERTEX> vs(map);			//marker for new vertices from edge cut
	CellMarkerStore<FACE> cf(map);
121 122
	Dart dPath;

Thomas's avatar
Thomas committed
123 124
	MarkerForTraversor<typename PFP::MAP::ParentMap, EDGE > mte(map);
	MarkerForTraversor<typename PFP::MAP::ParentMap, FACE > mtf(map);
125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140

	//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;
141
				typename PFP::VEC3 vec(Surface::Geometry::vectorOutOfDart<PFP>(map,dd,position));
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 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185
				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;

186 187
				do
				{
188 189 190 191 192
					//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
193 194
						do
						{
195 196 197 198 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 225 226 227 228 229 230 231 232 233 234 235 236
							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;
}

237
template <typename PFP>
Thomas's avatar
Thomas committed
238
Dart sliceConvexVolume(typename PFP::MAP& map, VertexAttribute<typename PFP::VEC3>& position, Dart d, CellMarker<EDGE>& edgesToCut, CellMarker<VERTEX>& verticesToSplit)
239 240 241 242 243 244
{
	typedef typename PFP::VEC3 VEC3;

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

Thomas's avatar
Thomas committed
249 250
	MarkerForTraversor<typename PFP::MAP::ParentMap, EDGE > mte(map);
	MarkerForTraversor<typename PFP::MAP::ParentMap, FACE > mtf(map);
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 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336

	//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
337
std::vector<Dart> sliceConvexVolumes(typename PFP::MAP& map, VertexAttribute<typename PFP::VEC3>& position,CellMarker<VOLUME>& volumesToCut, CellMarker<EDGE>& edgesToCut, CellMarker<VERTEX>& verticesToSplit)
338
{
Thomas's avatar
Thomas committed
339
    std::vector<Dart> vRes;
340 341

    typedef typename PFP::VEC3 VEC3;
Thomas's avatar
Thomas committed
342
    CellMarker<VERTEX> localVerticesToSplit(map); //marker for new vertices from edge cut
343 344 345

    //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
346
    CellMarkerStore<FACE> cf(map);
347 348 349 350 351 352

    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;
353 354 355 356 357

            //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
358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378

            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
379
                if(localVerticesToSplit.isMarked(dS) || verticesToSplit.isMarked(dS))
380
                {
381 382
                	//start from phi1(phi1()) to avoid the creation of faces of degree 2
                    Dart dSS = map.phi1(map.phi1(dS));
383
                    //search an other new vertex (or an existing vertex to split) in order to split the face
384

385 386
                    do
                    {
387 388
                        if((localVerticesToSplit.isMarked(dSS) || verticesToSplit.isMarked(dSS))
                        		&& !map.sameVertex(dS,dSS))
389 390 391 392 393 394
                        {
                            map.splitFace(dS,dSS);
                            split=true;
                        }
                        dSS = map.phi1(dSS);
                    } while(!split && dSS!=dS);
395
                    split=true; //go out of the first loop if no split case has been found
396 397 398 399 400 401 402 403 404 405
                }
                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
    {
406 407
        if(volumesToCut.isMarked(d))
        {
408 409 410 411 412 413 414 415 416
            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))
                {
417
                    Dart ddd = dd;
418
                    while(!localVerticesToSplit.isMarked(map.phi1(ddd))
419
                    		&& !verticesToSplit.isMarked(map.phi1(ddd)))
420 421 422 423 424 425 426
                        ddd = map.phi1(map.phi2(ddd));
                    found=true;
                    dPath=ddd;
                }
            }
            //define the path to split
            std::vector<Dart> vPath;
427
            vPath.reserve(32);
428
            vPath.push_back(dPath);
Thomas's avatar
Thomas committed
429
            CellMarker<FACE> cmf(map);
430

431

432
            //define the path to split for the whole volume
433 434
            bool pathFound=false;
            for(std::vector<Dart>::iterator it = vPath.begin() ; !pathFound && it != vPath.end(); ++it)
435 436 437
            {
                Dart dd = map.phi1(*it);

438 439 440
                if(map.sameVertex(dd,*vPath.begin()))
                	pathFound=true;
                else
441
                {
442 443 444 445 446 447
                	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);
448 449 450 451 452
                }
            }

            assert(vPath.size()>2);
            map.splitVolume(vPath);
453
            vRes.push_back(map.phi2(*vPath.begin()));
454 455 456
        }
    }

457
    return vRes;
458 459
}

460
template <typename PFP, typename EMBV>
461
void catmullClarkVol(typename PFP::MAP& map, EMBV& attributs)
untereiner's avatar
untereiner committed
462
{
463 464
	typedef typename EMBV::DATA_TYPE EMB;

465
	//std::vector<Dart> l_centers;
untereiner's avatar
untereiner committed
466 467 468
	std::vector<Dart> l_vertices;

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

472 473 474 475 476
	//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
477 478
	{
		//memorize each vertices per volumes
479
		if( !mv.isMarked(d))
untereiner's avatar
untereiner committed
480 481
		{
			l_vertices.push_back(d);
Pierre Kraemer's avatar
Pierre Kraemer committed
482
			mv.markOrbit<PFP::MAP::VERTEX_OF_PARENT>(d);
untereiner's avatar
untereiner committed
483 484
		}

485 486 487
		Dart f = map.phi1(d);
		map.cutEdge(d) ;
		Dart e = map.phi1(d) ;
untereiner's avatar
untereiner committed
488

489 490 491
		attributs[e] =  attributs[d];
		attributs[e] += attributs[f];
		attributs[e] *= 0.5;
untereiner's avatar
untereiner committed
492

untereiner's avatar
untereiner committed
493 494
		travE.skip(d) ;
		travE.skip(e) ;
untereiner's avatar
untereiner committed
495 496
	}

497 498 499
	//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
500
	{
501
		EMB center = Surface::Geometry::faceCentroid<PFP,EMBV>(map,d,attributs);
untereiner's avatar
untereiner committed
502

503 504 505
		Dart dd = map.phi1(d) ;
		Dart next = map.phi1(map.phi1(dd)) ;
		map.splitFace(dd, next) ;
untereiner's avatar
untereiner committed
506

507 508
		Dart ne = map.phi2(map.phi_1(dd)) ;
		map.cutEdge(ne) ;
untereiner's avatar
untereiner committed
509
		travF.skip(dd) ;
untereiner's avatar
untereiner committed
510

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

513 514 515 516 517
		dd = map.phi1(map.phi1(next)) ;
		while(dd != ne)
		{
			Dart tmp = map.phi1(ne) ;
			map.splitFace(tmp, dd) ;
untereiner's avatar
untereiner committed
518
			travF.skip(tmp) ;
519
			dd = map.phi1(map.phi1(dd)) ;
untereiner's avatar
untereiner committed
520
		}
521

untereiner's avatar
untereiner committed
522
		travF.skip(ne) ;
untereiner's avatar
untereiner committed
523 524
	}

525 526 527 528
	//3. create inside volumes

	std::vector<std::pair<Dart, Dart> > subdividedFaces;
	subdividedFaces.reserve(2048);
untereiner's avatar
untereiner committed
529 530
	for (std::vector<Dart>::iterator it = l_vertices.begin(); it != l_vertices.end(); ++it)
	{
531 532 533
		Dart e = *it;
		std::vector<Dart> v ;

untereiner's avatar
untereiner committed
534 535
		do
		{
536 537 538 539 540
			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
541

542 543
			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
544

545
			e = map.phi2(map.phi_1(e));
untereiner's avatar
untereiner committed
546
		}
547
		while(e != *it);
untereiner's avatar
untereiner committed
548

549 550 551
		//
		// SplitSurfaceInVolume
		//
untereiner's avatar
untereiner committed
552

553 554
		std::vector<Dart> vd2 ;
		vd2.reserve(v.size());
untereiner's avatar
untereiner committed
555

556 557 558 559 560
		// 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
561

562
		assert(vd2.size() == v.size());
untereiner's avatar
untereiner committed
563

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

566 567 568 569 570
		// 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
571

572
			// embed the vertex embedded from the origin volume to the new darts
Pierre Kraemer's avatar
Pierre Kraemer committed
573
			if(map.template isOrbitEmbedded<VERTEX>())
574
			{
Pierre Kraemer's avatar
Pierre Kraemer committed
575 576
				map.template copyDartEmbedding<VERTEX>(map.phi2(dit), map.phi1(dit));
				map.template copyDartEmbedding<VERTEX>(map.phi2(dit2), map.phi1(dit2));
577
			}
untereiner's avatar
untereiner committed
578

579
			// embed the edge embedded from the origin volume to the new darts
Pierre Kraemer's avatar
Pierre Kraemer committed
580
			if(map.template isOrbitEmbedded<EDGE>())
581
			{
Pierre Kraemer's avatar
Pierre Kraemer committed
582 583 584
				unsigned int eEmb = map.template getEmbedding<EDGE>(dit) ;
				map.template setDartEmbedding<EDGE>(map.phi2(dit), eEmb);
				map.template setDartEmbedding<EDGE>(map.phi2(dit2), eEmb);
585
			}
untereiner's avatar
untereiner committed
586

587
			// embed the volume embedded from the origin volume to the new darts
Pierre Kraemer's avatar
Pierre Kraemer committed
588
			if(map.template isOrbitEmbedded<VOLUME>())
589
			{
Pierre Kraemer's avatar
Pierre Kraemer committed
590 591
				map.template copyDartEmbedding<VOLUME>(map.phi2(dit), dit);
				map.template copyDartEmbedding<VOLUME>(map.phi2(dit2), dit2);
592 593
			}
		}
untereiner's avatar
untereiner committed
594

595 596 597
		//
		//
		//
untereiner's avatar
untereiner committed
598

599 600 601
		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
602

Pierre Kraemer's avatar
Pierre Kraemer committed
603
		if (map.template isOrbitEmbedded<VERTEX>())
604
		{
Pierre Kraemer's avatar
Pierre Kraemer committed
605 606
			map.template copyDartEmbedding<VERTEX>(map.phi_1(next), dd) ;
			map.template copyDartEmbedding<VERTEX>(map.phi_1(dd), next) ;
607 608 609 610
		}

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

612 613
//		dd = map.phi1(map.phi1(next)) ;
//		while(dd != ne)
untereiner's avatar
untereiner committed
614
//		{
615 616 617
//			Dart tmp = map.phi1(ne) ;
//			map.PFP::MAP::ParentMap::splitFace(tmp, dd);
//
618
//			if (map.isOrbitEmbedded<VERTEX>())
619
//			{
Pierre Kraemer's avatar
Pierre Kraemer committed
620 621
//				map.copyDartEmbedding<VERTEX>(map.phi_1(dd), tmp) ;
//				map.copyDartEmbedding<VERTEX>(map.phi_1(tmp), dd) ;
622 623 624
//			}
//
//			dd = map.phi1(map.phi1(dd)) ;
untereiner's avatar
untereiner committed
625
//		}
626 627
//
	}
untereiner's avatar
untereiner committed
628

629 630 631
//		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
632
//		{
633 634 635 636 637 638
//			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
639
//		}
640
//		setOrbitEmbedding<VERTEX>(centralDart, getEmbedding<VERTEX>(centralDart));
641 642 643 644 645 646 647 648 649 650
		//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())
//	{
651
//		map.setOrbitEmbedding<VERTEX>(map.phi1(d), map.getEmbedding<VERTEX>(map.phi1(d)));
652 653 654 655 656
//	}
//
//	TraversorF<typename PFP::MAP> travF2(map) ;
//	for (Dart d = travF2.begin(); d != travF2.end(); d = travF2.next())
//	{
657
//		map.setOrbitEmbedding<VERTEX>(map.phi2(map.phi1(d)), map.getEmbedding<VERTEX>(map.phi2(map.phi1(d))));
658
//	}
659 660
}

661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681
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 ;
	}
}
682 683

template <typename PFP>
684
void sqrt3Vol(typename PFP::MAP& map, VertexAttribute<typename PFP::VEC3>& position)
685 686
{
	DartMarkerStore m(map);
687

688 689
	DartMarkerStore newBoundaryV(map);

690 691 692
	//
	// 1-4 flip of all tetrahedra
	//
693
	TraversorW<typename PFP::MAP> tW(map);
694 695
	for(Dart dit = tW.begin() ; dit != tW.end() ; dit = tW.next())
	{
696 697 698
		Traversor3WF<typename PFP::MAP> tWF(map, dit);
		for(Dart ditWF = tWF.begin() ; ditWF != tWF.end() ; ditWF = tWF.next())
		{
699
			if(!map.isBoundaryFace(ditWF) && !m.isMarked(ditWF))
700 701
				m.markOrbit<FACE>(ditWF);
		}
702

703 704 705 706 707 708 709
		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;

710
		Dart dres = Volume::Modelisation::Tetrahedralization::flip1To4<PFP>(map, dit);
711
		position[dres] = volCenter;
712
	}
713

714 715 716
	//
	// 2-3 swap of all old interior faces
	//
717
	TraversorF<typename PFP::MAP> tF(map);
718 719 720 721 722
	for(Dart dit = tF.begin() ; dit != tF.end() ; dit = tF.next())
	{
		if(m.isMarked(dit))
		{
			m.unmarkOrbit<FACE>(dit);
723
			Volume::Modelisation::Tetrahedralization::swap2To3<PFP>(map, dit);
724 725 726 727 728 729
		}
	}

	//
	// 1-3 flip of all boundary tetrahedra
	//
730
	TraversorW<typename PFP::MAP> tWb(map);
731 732 733 734
	for(Dart dit = tWb.begin() ; dit != tWb.end() ; dit = tWb.next())
	{
		if(map.isBoundaryVolume(dit))
		{
735 736 737
			Traversor3WE<typename PFP::MAP> tWE(map, dit);
			for(Dart ditWE = tWE.begin() ; ditWE != tWE.end() ; ditWE = tWE.next())
			{
738
				if(map.isBoundaryEdge(ditWE) && !m.isMarked(ditWE))
739 740 741
					m.markOrbit<EDGE>(ditWE);
			}

742 743 744 745 746 747
			typename PFP::VEC3 faceCenter(0.0);
			faceCenter += position[dit];
			faceCenter += position[map.phi1(dit)];
			faceCenter += position[map.phi_1(dit)];
			faceCenter /= 3;

748
			Dart dres = Volume::Modelisation::Tetrahedralization::flip1To3<PFP>(map, dit);
749
			position[dres] = faceCenter;
750 751

			newBoundaryV.markOrbit<VERTEX>(dres);
752 753 754
		}
	}

755
/*
Sylvain Thery's avatar
Sylvain Thery committed
756
	TraversorV<typename PFP::MAP> tVg(map);
757 758 759 760
	for(Dart dit = tVg.begin() ; dit != tVg.end() ; dit = tVg.next())
	{
		if(map.isBoundaryVertex(dit) && !newBoundaryV.isMarked(dit))
		{
761 762 763
			Dart db = map.findBoundaryFaceOfVertex(dit);

			typename PFP::VEC3 P = position[db] ;
764 765
			typename PFP::VEC3 newP(0) ;
			unsigned int val = 0 ;
766
			Dart vit = db ;
767 768
			do
			{
769 770
				//newP += position[map.phi_1(map.phi2(map.phi1(vit)))] ;
				newP += position[map.phi2(vit)];
771 772
				++val ;
				vit = map.phi2(map.phi_1(vit)) ;
773
			} while(vit != db) ;
774 775 776 777 778
			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 ;
779 780 781
			position[db] = newP ;
		}
	}
782
*/
783

784
/*
785 786 787
	//
	// edge-removal on all old boundary edges
	//
untereiner's avatar
untereiner committed
788
	TraversorE<typename PFP::MAP> tE(map);
789 790 791 792 793 794 795
	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)));
			Volume::Modelisation::Tetrahedralization::swapGen3To2<PFP>(map, d);
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
//	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))
//		{
//			Dart db = map.findBoundaryFaceOfVertex(dit);
//
//			typename PFP::VEC3 P = position[db] ;
//			typename PFP::VEC3 newP(0) ;
//			unsigned int val = 0 ;
//			Dart vit = db ;
//			do
//			{
//				newP += position[map.phi_1(map.phi2(map.phi1(vit)))] ;
//				++val ;
//				vit = map.phi2(map.phi_1(vit)) ;
//			} while(vit != db) ;
//			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[db] = newP ;
//		}
//	}

829 830 831 832 833 834 835 836 837 838 839 840 841 842
	//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


843 844
//	VertexAutoAttribute<typename PFP::VEC3> diffCoord(map, "diffCoord");
//	Algo::Volume::Geometry::computeLaplacianTopoVertices<PFP>(map, position, diffCoord) ;
845
//
846 847 848 849 850 851 852 853 854
//	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())
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
//		if(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) ;
882
//	}
883 884
//
//	nlDeleteContext(nlContext);
untereiner's avatar
untereiner committed
885 886
}

887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906

// solving Ax = b

template <typename PFP>
void relaxation(typename PFP::MAP& map, VertexAttribute<typename PFP::VEC3>& position)
{
	VertexAttribute<unsigned int> indexV = map.template getAttribute<unsigned int, VERTEX>("indexV");
	if(!indexV.isValid())
		indexV = map.template addAttribute<unsigned int, VERTEX>("indexV");

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

	//uniform weight
	float weight = 1.0;

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

907
//	nlMakeCurrent(nlContext);
908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927
	if(nlGetCurrentState() == NL_STATE_INITIAL)
		nlBegin(NL_SYSTEM) ;

	for(unsigned int coord = 0; coord < 3; ++coord)
	{
		std::cout << "coord " << coord << std::flush;
		//setup variables
		TraversorV<typename PFP::MAP> tv(map);
		for(Dart dit = tv.begin() ; dit != tv.end() ; dit = tv.next())
		{
			nlSetVariable(indexV[dit], (position[dit])[coord]);

			if(map.isBoundaryVertex(dit))
				nlLockVariable(indexV[dit]);
		}

		std::cout << "... variables set... " << std::flush;

		nlBegin(NL_MATRIX) ;

928 929
		nlEnable(NL_NORMALIZE_ROWS) ;

930 931 932 933 934
		TraversorV<typename PFP::MAP> tv2(map);
		for(Dart dit = tv2.begin() ; dit != tv2.end() ; dit = tv2.next())
		{
			if(!map.isBoundaryVertex(dit))
			{
935
				nlRowParameterd(NL_RIGHT_HAND_SIDE, 0) ; //b[i]
936 937 938 939 940 941 942 943
				//nlRowParameterd(NL_ROW_SCALING, weight) ;

				nlBegin(NL_ROW) ;

				float sum = 0;
				Traversor3VVaE<typename PFP::MAP> tvvae(map, dit);
				for(Dart ditvvae = tvvae.begin() ; ditvvae != tvvae.end() ; ditvvae = tvvae.next())
				{
944
					nlCoefficient(indexV[ditvvae], weight);
945 946 947 948 949 950 951 952
					sum += weight;
				}

				nlCoefficient(indexV[dit], -sum) ;
				nlEnd(NL_ROW) ;
			}
		}

953 954
		nlDisable(NL_NORMALIZE_ROWS) ;

955 956 957 958 959 960 961 962 963 964 965 966 967 968 969
		nlEnd(NL_MATRIX) ;

		nlEnd(NL_SYSTEM) ;
		std::cout << "... system built... " << std::flush;

		nlSolve();
		std::cout << "... system solved... " << std::flush;

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

970
		nlReset(NL_TRUE) ;
971 972 973 974 975 976
		std::cout << "... done" << std::endl;
	}

	nlDeleteContext(nlContext);
}

977 978 979 980 981 982 983 984 985 986 987 988 989 990 991
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
992

993

994
} // namespace Modelisation
untereiner's avatar
untereiner committed
995

996
} // namespace volume
997

998
} // namespace Algo
untereiner's avatar
untereiner committed
999

1000
} // namespace CGoGN
untereiner's avatar
untereiner committed
1001