subdivision3.hpp 27.8 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
template <typename PFP>
Sylvain Thery's avatar
Sylvain Thery committed
46
bool isHexahedron(typename PFP::MAP& the_map, Dart d)
47 48 49 50
{
    unsigned int nbFaces = 0;

    //Test the number of faces end its valency
Sylvain Thery's avatar
Sylvain Thery committed
51
	Traversor3WF<typename PFP::MAP> travWF(the_map, d, false);
52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
    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>
114
Dart sliceConvexVolume(typename PFP::MAP& map, VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position, Dart d, Geom::Plane3D<typename PFP::REAL> pl)
115
{
116 117 118 119
	typedef typename PFP::MAP MAP;
	typedef typename PFP::VEC3 VEC3;
	typedef typename PFP::REAL REAL;

120 121 122
	Dart dRes=NIL;
	unsigned int nbInter = 0;
	unsigned int nbVertices = 0;
123 124
	CellMarkerStore<MAP, VERTEX> vs(map);			//marker for new vertices from edge cut
	CellMarkerStore<MAP, FACE> cf(map);
125 126
	Dart dPath;

127 128
	MarkerForTraversor<typename MAP::ParentMap, EDGE> mte(map);
	MarkerForTraversor<typename MAP::ParentMap, FACE> mtf(map);
129 130

	//search edges and vertices crossing the plane
131
	Traversor3WE<typename MAP::ParentMap> te(map,d);
132 133 134 135 136 137 138 139 140 141 142 143
	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
			{
144 145 146
				VEC3 interP;
				VEC3 vec(Surface::Geometry::vectorOutOfDart<PFP>(map,dd,position));
				Geom::Intersection inter = Geom::intersectionLinePlane<VEC3, typename Geom::Plane3D<REAL> >(position[dd],vec,pl,interP);
147 148 149 150

				if(inter==Geom::FACE_INTERSECTION)
				{
					Dart dOp = map.phi1(dd);
151 152
					VEC3 v2(interP-position[dd]);
					VEC3 v3(interP-position[dOp]);
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
					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)
	{
179
		Traversor3WF<typename MAP::ParentMap > tf(map,d);
180 181 182 183 184 185 186 187 188 189
		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;

190 191
				do
				{
192 193 194 195 196
					//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
197 198
						do
						{
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 237 238 239 240
							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;
}

241
template <typename PFP>
242
Dart sliceConvexVolume(typename PFP::MAP& map, VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position, Dart d, CellMarker<typename PFP::MAP, EDGE>& edgesToCut, CellMarker<typename PFP::MAP, VERTEX>& verticesToSplit)
243
{
244
	typedef typename PFP::MAP MAP;
245
	typedef typename PFP::VEC3 VEC3;
246
	typedef typename PFP::REAL REAL;
247 248 249 250

	Dart dRes;
	unsigned int nbInter = 0;
	unsigned int nbVertices = 0;
251 252
	CellMarkerStore<MAP, VERTEX> vs(map);			//marker for new vertices from edge cut
	CellMarkerStore<MAP, FACE> cf(map);
253 254
	Dart dPath;

255 256
	MarkerForTraversor<typename MAP::ParentMap, EDGE> mte(map);
	MarkerForTraversor<typename MAP::ParentMap, FACE> mtf(map);
257 258

	//search edges and vertices crossing the plane
259 260
	Traversor3WE<typename MAP::ParentMap > te(map,d);
	for(Dart dd = te.begin(); dd != te.end(); dd = te.next())
261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284
	{
		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);

285
	Traversor3WF<typename MAP::ParentMap > tf(map,d);
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 337 338 339 340 341 342
	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>
343
std::vector<Dart> sliceConvexVolumes(typename PFP::MAP& map, VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position, CellMarker<typename PFP::MAP, VOLUME>& volumesToCut, CellMarker<typename PFP::MAP, EDGE>& edgesToCut, CellMarker<typename PFP::MAP, VERTEX>& verticesToSplit)
344
{
345 346 347 348
	typedef typename PFP::MAP MAP;
	typedef typename PFP::VEC3 VEC3;
	typedef typename PFP::REAL REAL;

Thomas's avatar
Thomas committed
349
    std::vector<Dart> vRes;
350

351
	CellMarker<MAP, VERTEX> localVerticesToSplit(map); //marker for new vertices from edge cut
352 353

    //Step 1: Cut the edges and mark the resulting vertices as vertices to be face-split
354 355
	TraversorE<MAP> te(map);
	CellMarkerStore<MAP, FACE> cf(map);
356 357 358 359 360 361

    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;
362 363

            //turn around the edge and mark for future split face
364
			Traversor3EF<MAP> t3ef(map,d);
365 366
            for(Dart dd = t3ef.begin() ; dd != t3ef.end() ; dd = t3ef.next())
            	cf.mark(dd);			//mark face to split
367 368 369 370 371 372 373 374 375 376

            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
377
	TraversorF<MAP> tf(map);
378 379 380 381 382 383 384 385 386 387
    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
388
                if(localVerticesToSplit.isMarked(dS) || verticesToSplit.isMarked(dS))
389
                {
390 391
                	//start from phi1(phi1()) to avoid the creation of faces of degree 2
                    Dart dSS = map.phi1(map.phi1(dS));
392
                    //search an other new vertex (or an existing vertex to split) in order to split the face
393

394 395
                    do
                    {
396 397
                        if((localVerticesToSplit.isMarked(dSS) || verticesToSplit.isMarked(dSS))
                        		&& !map.sameVertex(dS,dSS))
398 399 400 401 402 403
                        {
                            map.splitFace(dS,dSS);
                            split=true;
                        }
                        dSS = map.phi1(dSS);
                    } while(!split && dSS!=dS);
404
                    split=true; //go out of the first loop if no split case has been found
405 406 407 408 409 410 411
                }
                dS = map.phi1(dS);
            } while(!split && dS!=d);
        }
    }

    //Step 3 : Find path and split volumes
412
	TraversorW<MAP> tw(map);
413 414
    for(Dart d = tw.begin(); d != tw.end(); d=tw.next()) //Parcours des volumes
    {
415 416
        if(volumesToCut.isMarked(d))
        {
417
			Traversor3WV<MAP> t3wv(map,d);
418 419 420 421 422 423 424 425
            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))
                {
426
                    Dart ddd = dd;
427
                    while(!localVerticesToSplit.isMarked(map.phi1(ddd))
428
                    		&& !verticesToSplit.isMarked(map.phi1(ddd)))
429 430 431 432 433 434 435
                        ddd = map.phi1(map.phi2(ddd));
                    found=true;
                    dPath=ddd;
                }
            }
            //define the path to split
            std::vector<Dart> vPath;
436
            vPath.reserve(32);
437
            vPath.push_back(dPath);
438
			CellMarker<MAP, FACE> cmf(map);
439

440
            //define the path to split for the whole volume
441 442
            bool pathFound=false;
            for(std::vector<Dart>::iterator it = vPath.begin() ; !pathFound && it != vPath.end(); ++it)
443 444 445
            {
                Dart dd = map.phi1(*it);

446 447 448
                if(map.sameVertex(dd,*vPath.begin()))
                	pathFound=true;
                else
449
                {
450 451 452 453 454 455
                	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);
456 457 458 459 460
                }
            }

            assert(vPath.size()>2);
            map.splitVolume(vPath);
461
            vRes.push_back(map.phi2(*vPath.begin()));
462 463 464
        }
    }

465
    return vRes;
466 467
}

468
template <typename PFP, typename EMBV>
469
void catmullClarkVol(typename PFP::MAP& map, EMBV& attributs)
untereiner's avatar
untereiner committed
470
{
471
	typedef typename PFP::MAP MAP;
472 473
	typedef typename EMBV::DATA_TYPE EMB;

474
	//std::vector<Dart> l_centers;
untereiner's avatar
untereiner committed
475 476 477
	std::vector<Dart> l_vertices;

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

481 482
	//subdivision
	//1. cut edges
483 484
	DartMarkerNoUnmark<MAP> mv(map);
	TraversorE<MAP> travE(map);
485
	for (Dart d = travE.begin(); d != travE.end(); d = travE.next())
untereiner's avatar
untereiner committed
486 487
	{
		//memorize each vertices per volumes
488
		if( !mv.isMarked(d))
untereiner's avatar
untereiner committed
489 490
		{
			l_vertices.push_back(d);
491
			mv.template markOrbit<PFP::MAP::VERTEX_OF_PARENT>(d);
untereiner's avatar
untereiner committed
492 493
		}

494 495 496
		Dart f = map.phi1(d);
		map.cutEdge(d) ;
		Dart e = map.phi1(d) ;
untereiner's avatar
untereiner committed
497

498 499 500
		attributs[e] =  attributs[d];
		attributs[e] += attributs[f];
		attributs[e] *= 0.5;
untereiner's avatar
untereiner committed
501

untereiner's avatar
untereiner committed
502 503
		travE.skip(d) ;
		travE.skip(e) ;
untereiner's avatar
untereiner committed
504 505
	}

506
	//2. split faces - quadrangule faces
507
	TraversorF<MAP> travF(map) ;
508
	for (Dart d = travF.begin(); d != travF.end(); d = travF.next())
untereiner's avatar
untereiner committed
509
	{
510
		EMB center = Surface::Geometry::faceCentroid<PFP,EMBV>(map,d,attributs);
untereiner's avatar
untereiner committed
511

512 513 514
		Dart dd = map.phi1(d) ;
		Dart next = map.phi1(map.phi1(dd)) ;
		map.splitFace(dd, next) ;
untereiner's avatar
untereiner committed
515

516 517
		Dart ne = map.phi2(map.phi_1(dd)) ;
		map.cutEdge(ne) ;
untereiner's avatar
untereiner committed
518
		travF.skip(dd) ;
untereiner's avatar
untereiner committed
519

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

522 523 524 525 526
		dd = map.phi1(map.phi1(next)) ;
		while(dd != ne)
		{
			Dart tmp = map.phi1(ne) ;
			map.splitFace(tmp, dd) ;
untereiner's avatar
untereiner committed
527
			travF.skip(tmp) ;
528
			dd = map.phi1(map.phi1(dd)) ;
untereiner's avatar
untereiner committed
529
		}
530

untereiner's avatar
untereiner committed
531
		travF.skip(ne) ;
untereiner's avatar
untereiner committed
532 533
	}

534 535 536 537
	//3. create inside volumes

	std::vector<std::pair<Dart, Dart> > subdividedFaces;
	subdividedFaces.reserve(2048);
untereiner's avatar
untereiner committed
538 539
	for (std::vector<Dart>::iterator it = l_vertices.begin(); it != l_vertices.end(); ++it)
	{
540 541 542
		Dart e = *it;
		std::vector<Dart> v ;

untereiner's avatar
untereiner committed
543 544
		do
		{
545 546 547 548 549
			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
550

551 552
			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
553

554
			e = map.phi2(map.phi_1(e));
untereiner's avatar
untereiner committed
555
		}
556
		while(e != *it);
untereiner's avatar
untereiner committed
557

558 559 560
		//
		// SplitSurfaceInVolume
		//
untereiner's avatar
untereiner committed
561

562 563
		std::vector<Dart> vd2 ;
		vd2.reserve(v.size());
untereiner's avatar
untereiner committed
564

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

571
		assert(vd2.size() == v.size());
untereiner's avatar
untereiner committed
572

573
		map.MAP::ParentMap::splitSurface(v, true, false);
untereiner's avatar
untereiner committed
574

575 576 577 578 579
		// 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
580

581
			// embed the vertex embedded from the origin volume to the new darts
582
			if(map.template isOrbitEmbedded<VERTEX>())
583
			{
584 585
				map.template copyDartEmbedding<VERTEX>(map.phi2(dit), map.phi1(dit));
				map.template copyDartEmbedding<VERTEX>(map.phi2(dit2), map.phi1(dit2));
586
			}
untereiner's avatar
untereiner committed
587

588
			// embed the edge embedded from the origin volume to the new darts
589
			if(map.template isOrbitEmbedded<EDGE>())
590
			{
591 592 593
				unsigned int eEmb = map.template getEmbedding<EDGE>(dit) ;
				map.template setDartEmbedding<EDGE>(map.phi2(dit), eEmb);
				map.template setDartEmbedding<EDGE>(map.phi2(dit2), eEmb);
594
			}
untereiner's avatar
untereiner committed
595

596
			// embed the volume embedded from the origin volume to the new darts
597
			if(map.template isOrbitEmbedded<VOLUME>())
598
			{
599 600
				map.template copyDartEmbedding<VOLUME>(map.phi2(dit), dit);
				map.template copyDartEmbedding<VOLUME>(map.phi2(dit2), dit2);
601 602
			}
		}
untereiner's avatar
untereiner committed
603

604 605 606
		//
		//
		//
untereiner's avatar
untereiner committed
607

608 609
		Dart dd = map.phi2(map.phi1(*it));
		Dart next = map.phi1(map.phi1(dd)) ;
610
		map.MAP::ParentMap::splitFace(dd, next);
untereiner's avatar
untereiner committed
611

612
		if (map.template isOrbitEmbedded<VERTEX>())
613
		{
614 615
			map.template copyDartEmbedding<VERTEX>(map.phi_1(next), dd) ;
			map.template copyDartEmbedding<VERTEX>(map.phi_1(dd), next) ;
616 617 618
		}

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

621 622
//		dd = map.phi1(map.phi1(next)) ;
//		while(dd != ne)
untereiner's avatar
untereiner committed
623
//		{
624 625 626
//			Dart tmp = map.phi1(ne) ;
//			map.PFP::MAP::ParentMap::splitFace(tmp, dd);
//
627
//			if (map.isOrbitEmbedded<VERTEX>())
628
//			{
Pierre Kraemer's avatar
Pierre Kraemer committed
629 630
//				map.copyDartEmbedding<VERTEX>(map.phi_1(dd), tmp) ;
//				map.copyDartEmbedding<VERTEX>(map.phi_1(tmp), dd) ;
631 632 633
//			}
//
//			dd = map.phi1(map.phi1(dd)) ;
untereiner's avatar
untereiner committed
634
//		}
635 636
//
	}
untereiner's avatar
untereiner committed
637

638 639 640
//		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
641
//		{
642 643 644 645 646 647
//			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
648
//		}
649
//		setOrbitEmbedding<VERTEX>(centralDart, getEmbedding<VERTEX>(centralDart));
650 651 652 653 654 655 656 657 658 659
		//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())
//	{
660
//		map.setOrbitEmbedding<VERTEX>(map.phi1(d), map.getEmbedding<VERTEX>(map.phi1(d)));
661 662 663 664 665
//	}
//
//	TraversorF<typename PFP::MAP> travF2(map) ;
//	for (Dart d = travF2.begin(); d != travF2.end(); d = travF2.next())
//	{
666
//		map.setOrbitEmbedding<VERTEX>(map.phi2(map.phi1(d)), map.getEmbedding<VERTEX>(map.phi2(map.phi1(d))));
667
//	}
668 669
}

670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690
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 ;
	}
}
691 692

template <typename PFP>
693
void sqrt3Vol(typename PFP::MAP& map, VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position)
694
{
695 696
	typedef typename PFP::MAP MAP;
	typedef typename PFP::VEC3 VEC3;
697

698 699 700
	DartMarkerStore<MAP> m(map);

	DartMarkerStore<MAP> newBoundaryV(map);
701

702 703 704
	//
	// 1-4 flip of all tetrahedra
	//
705
	TraversorW<MAP> tW(map);
706 707
	for(Dart dit = tW.begin() ; dit != tW.end() ; dit = tW.next())
	{
708
		Traversor3WF<MAP> tWF(map, dit);
709 710
		for(Dart ditWF = tWF.begin() ; ditWF != tWF.end() ; ditWF = tWF.next())
		{
711
			if(!map.isBoundaryFace(ditWF) && !m.isMarked(ditWF))
712 713
				m.markOrbit<FACE>(ditWF);
		}
714

715
		VEC3 volCenter(0.0);
716 717 718 719 720 721
		volCenter += position[dit];
		volCenter += position[map.phi1(dit)];
		volCenter += position[map.phi_1(dit)];
		volCenter += position[map.phi_1(map.phi2(dit))];
		volCenter /= 4;

722
		Dart dres = Volume::Modelisation::Tetrahedralization::flip1To4<PFP>(map, dit);
723
		position[dres] = volCenter;
724
	}
725

726 727 728
	//
	// 2-3 swap of all old interior faces
	//
729
	TraversorF<MAP> tF(map);
730 731 732 733 734
	for(Dart dit = tF.begin() ; dit != tF.end() ; dit = tF.next())
	{
		if(m.isMarked(dit))
		{
			m.unmarkOrbit<FACE>(dit);
735
			Volume::Modelisation::Tetrahedralization::swap2To3<PFP>(map, dit);
736 737 738 739 740 741
		}
	}

	//
	// 1-3 flip of all boundary tetrahedra
	//
742
	TraversorW<MAP> tWb(map);
743 744
	for(Dart dit = tWb.begin() ; dit != tWb.end() ; dit = tWb.next())
	{
Sylvain Thery's avatar
Tutos  
Sylvain Thery committed
745
		if(map.isBoundaryAdjacentVolume(dit))
746
		{
747
			Traversor3WE<MAP> tWE(map, dit);
748 749
			for(Dart ditWE = tWE.begin() ; ditWE != tWE.end() ; ditWE = tWE.next())
			{
750
				if(map.isBoundaryEdge(ditWE) && !m.isMarked(ditWE))
751 752 753
					m.markOrbit<EDGE>(ditWE);
			}

754
			VEC3 faceCenter(0.0);
755 756 757 758 759
			faceCenter += position[dit];
			faceCenter += position[map.phi1(dit)];
			faceCenter += position[map.phi_1(dit)];
			faceCenter /= 3;

760
			Dart dres = Volume::Modelisation::Tetrahedralization::flip1To3<PFP>(map, dit);
761
			position[dres] = faceCenter;
762 763

			newBoundaryV.markOrbit<VERTEX>(dres);
764 765 766
		}
	}

767
/*
Sylvain Thery's avatar
Sylvain Thery committed
768
	TraversorV<typename PFP::MAP> tVg(map);
769 770 771 772
	for(Dart dit = tVg.begin() ; dit != tVg.end() ; dit = tVg.next())
	{
		if(map.isBoundaryVertex(dit) && !newBoundaryV.isMarked(dit))
		{
773 774 775
			Dart db = map.findBoundaryFaceOfVertex(dit);

			typename PFP::VEC3 P = position[db] ;
776 777
			typename PFP::VEC3 newP(0) ;
			unsigned int val = 0 ;
778
			Dart vit = db ;
779 780
			do
			{
781 782
				//newP += position[map.phi_1(map.phi2(map.phi1(vit)))] ;
				newP += position[map.phi2(vit)];
783 784
				++val ;
				vit = map.phi2(map.phi_1(vit)) ;
785
			} while(vit != db) ;
786 787 788 789 790
			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 ;
791 792 793
			position[db] = newP ;
		}
	}
794
*/
795

796
/*
797 798 799
	//
	// edge-removal on all old boundary edges
	//
untereiner's avatar
untereiner committed
800
	TraversorE<typename PFP::MAP> tE(map);
801 802 803 804 805 806 807
	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);
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
//	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 ;
//		}
//	}

841 842 843 844 845 846 847 848 849 850 851 852 853 854
	//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


855 856
//	VertexAutoAttribute<typename PFP::VEC3> diffCoord(map, "diffCoord");
//	Algo::Volume::Geometry::computeLaplacianTopoVertices<PFP>(map, position, diffCoord) ;
857
//
858 859 860 861 862 863 864 865 866
//	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())
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
//		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) ;
894
//	}
895 896
//
//	nlDeleteContext(nlContext);
untereiner's avatar
untereiner committed
897 898
}

899 900 901 902

// solving Ax = b

template <typename PFP>
903
void relaxation(typename PFP::MAP& map, VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position)
904
{
905 906 907
	typedef typename PFP::MAP MAP;
	typedef typename PFP::VEC3 VEC3;

908
	VertexAttribute<unsigned int, MAP> indexV = map.template getAttribute<unsigned int, VERTEX>("indexV");
909 910 911 912 913 914 915 916 917 918 919 920 921
	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);

922
//	nlMakeCurrent(nlContext);
923 924 925 926 927 928 929
	if(nlGetCurrentState() == NL_STATE_INITIAL)
		nlBegin(NL_SYSTEM) ;

	for(unsigned int coord = 0; coord < 3; ++coord)
	{
		std::cout << "coord " << coord << std::flush;
		//setup variables
930
		TraversorV<MAP> tv(map);
931 932 933 934 935 936 937 938 939 940 941 942
		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) ;

943 944
		nlEnable(NL_NORMALIZE_ROWS) ;

945
		TraversorV<MAP> tv2(map);
946 947 948 949
		for(Dart dit = tv2.begin() ; dit != tv2.end() ; dit = tv2.next())
		{
			if(!map.isBoundaryVertex(dit))
			{
950
				nlRowParameterd(NL_RIGHT_HAND_SIDE, 0) ; //b[i]
951 952 953 954 955
				//nlRowParameterd(NL_ROW_SCALING, weight) ;

				nlBegin(NL_ROW) ;

				float sum = 0;
956
				Traversor3VVaE<MAP> tvvae(map, dit);
957 958
				for(Dart ditvvae = tvvae.begin() ; ditvvae != tvvae.end() ; ditvvae = tvvae.next())
				{
959
					nlCoefficient(indexV[ditvvae], weight);
960 961 962 963 964 965 966 967
					sum += weight;
				}

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

968 969
		nlDisable(NL_NORMALIZE_ROWS) ;

970 971 972 973 974 975 976 977 978 979 980 981 982 983 984
		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]);
		}

985
		nlReset(NL_TRUE) ;
986 987 988 989 990 991
		std::cout << "... done" << std::endl;
	}

	nlDeleteContext(nlContext);
}

992
template <typename PFP>
993
void computeDual(typename PFP::MAP& map, VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position)
994
{
995 996 997
	typedef typename PFP::MAP MAP;
	typedef typename PFP::VEC3 VEC3;

998
	// VolumeAttribute -> after dual new VertexAttribute
999
	VolumeAttribute<VEC3, MAP> positionV  = map.template getAttribute<VEC3, VOLUME>("position") ;
1000
	if(!positionV.isValid())
1001
		positionV = map.template addAttribute<VEC3, VOLUME>("position") ;
1002 1003 1004 1005 1006 1007 1008 1009

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

1011
} // namespace Modelisation
untereiner's avatar
untereiner committed
1012

1013
} // namespace volume
1014

1015
} // namespace Algo
untereiner's avatar
untereiner committed
1016

1017
} // namespace CGoGN