voronoiDiagrams.hpp 13 KB
Newer Older
Basile Sauvage's avatar
Basile Sauvage committed
1 2 3 4 5 6
namespace CGoGN
{

namespace Algo
{

Sylvain Thery's avatar
Sylvain Thery committed
7 8 9
namespace Surface
{

Basile Sauvage's avatar
Basile Sauvage committed
10 11 12
namespace Geometry
{

13 14 15 16
/***********************************************************
 * class VoronoiDiagram
 ***********************************************************/

Basile Sauvage's avatar
Basile Sauvage committed
17
template <typename PFP>
18
VoronoiDiagram<PFP>::VoronoiDiagram (MAP& m, const EdgeAttribute<REAL, MAP>& p, VertexAttribute<unsigned int, MAP>& r) :
Pierre Kraemer's avatar
Pierre Kraemer committed
19 20 21 22
	map(m),
	edgeCost (p),
	regions (r),
	vmReached(m)
Basile Sauvage's avatar
Basile Sauvage committed
23
{
Sauvage's avatar
Sauvage committed
24
	vertexInfo = map.template addAttribute<VertexInfo, VERTEX, typename PFP::MAP>("vertexInfo");
Basile Sauvage's avatar
Basile Sauvage committed
25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42
}

template <typename PFP>
VoronoiDiagram<PFP>::~VoronoiDiagram ()
{
	map.removeAttribute(vertexInfo);
}

template <typename PFP>
void VoronoiDiagram<PFP>::clear ()
{
	regions.setAllValues(0);
	border.clear();
	front.clear();
	vmReached.unmarkAll();
}

template <typename PFP>
Basile Sauvage's avatar
Basile Sauvage committed
43
void VoronoiDiagram<PFP>::setSeeds_fromVector (const std::vector<Dart>& s)
Basile Sauvage's avatar
Basile Sauvage committed
44
{
Basile Sauvage's avatar
Basile Sauvage committed
45
	seeds.clear();
Basile Sauvage's avatar
Basile Sauvage committed
46 47 48 49
	seeds = s;
}

template <typename PFP>
Basile Sauvage's avatar
Basile Sauvage committed
50
void VoronoiDiagram<PFP>::setSeeds_random (unsigned int nseeds)
Basile Sauvage's avatar
Basile Sauvage committed
51
{
Basile Sauvage's avatar
Basile Sauvage committed
52
	seeds.clear();
Basile Sauvage's avatar
Basile Sauvage committed
53
	srand ( time(NULL) );
54 55
	const unsigned int nbv = map.getNbCells(VERTEX);

Pierre Kraemer's avatar
Pierre Kraemer committed
56
	std::set<unsigned int> myVertices ;
57 58 59 60 61 62 63
	while (myVertices.size() < nseeds)
	{
		myVertices.insert(rand() % nbv);
	}

	std::set<unsigned int>::iterator it = myVertices.begin();
	unsigned int n = 0;
Pierre Kraemer's avatar
Pierre Kraemer committed
64
	TraversorV<MAP> tv (map);
65 66 67 68 69
	Dart dit = tv.begin();

	while (it != myVertices.end())
	{
		while(n<*it)
Basile Sauvage's avatar
Basile Sauvage committed
70
		{
71 72
			dit = tv.next();
			++n;
Basile Sauvage's avatar
Basile Sauvage committed
73
		}
74 75 76 77 78
		seeds.push_back(dit);
		it++;
	}

	// random permutation = un-sort the seeds
Pierre Kraemer's avatar
Pierre Kraemer committed
79
	for (unsigned int i = 0; i < nseeds; i++)
80 81 82 83 84
	{
		unsigned int j = i + rand() % (nseeds - i);
		Dart d = seeds[i];
		seeds[i] = seeds[j];
		seeds[j] = d;
Basile Sauvage's avatar
Basile Sauvage committed
85 86 87 88 89 90
	}
}

template <typename PFP>
void VoronoiDiagram<PFP>::initFrontWithSeeds ()
{
Basile Sauvage's avatar
Basile Sauvage committed
91 92
//	vmReached.unmarkAll();
	clear();
Basile Sauvage's avatar
Basile Sauvage committed
93 94 95 96 97 98
	for (unsigned int i = 0; i < seeds.size(); i++)
	{
		Dart d = seeds[i];
		vmReached.mark(d);
		vertexInfo[d].it = front.insert(std::pair<float,Dart>(0.0, d));
		vertexInfo[d].valid = true;
99 100
		regions[d] = i;
		vertexInfo[d].pathOrigin = d;
Basile Sauvage's avatar
Basile Sauvage committed
101 102 103 104
	}
}

template <typename PFP>
105
void VoronoiDiagram<PFP>::setCost (const EdgeAttribute<typename PFP::REAL,typename PFP::MAP>& c)
Pierre Kraemer's avatar
Pierre Kraemer committed
106
{
Basile Sauvage's avatar
Basile Sauvage committed
107 108 109
	edgeCost = c;
}

110
template <typename PFP>
Pierre Kraemer's avatar
Pierre Kraemer committed
111 112
void VoronoiDiagram<PFP>::collectVertexFromFront(Dart e)
{
113
	regions[e] = regions[vertexInfo[e].pathOrigin];
114 115
	front.erase(vertexInfo[e].it);
	vertexInfo[e].valid=false;
116 117 118
}

template <typename PFP>
Pierre Kraemer's avatar
Pierre Kraemer committed
119 120
void VoronoiDiagram<PFP>::addVertexToFront(Dart f, float d)
{
121 122
	VertexInfo& vi (vertexInfo[f]);
	vi.it = front.insert(std::pair<float,Dart>(d + edgeCost[f], f));
Pierre Kraemer's avatar
Pierre Kraemer committed
123
	vi.valid = true;
124 125 126 127 128
	vi.pathOrigin = map.phi2(f);
	vmReached.mark(f);
}

template <typename PFP>
Pierre Kraemer's avatar
Pierre Kraemer committed
129 130
void VoronoiDiagram<PFP>::updateVertexInFront(Dart f, float d)
{
131 132 133 134 135 136 137 138 139
	VertexInfo& vi (vertexInfo[f]);
	float dist = d + edgeCost[f];
	if (dist < vi.it->first)
	{
		front.erase(vi.it);
		vi.it = front.insert(std::pair<float,Dart>(dist, f));
		vi.pathOrigin = map.phi2(f);
	}
}
Basile Sauvage's avatar
Basile Sauvage committed
140 141

template <typename PFP>
Basile Sauvage's avatar
Basile Sauvage committed
142
Dart VoronoiDiagram<PFP>::computeDiagram ()
Basile Sauvage's avatar
Basile Sauvage committed
143 144 145
{
	initFrontWithSeeds();

Basile Sauvage's avatar
Basile Sauvage committed
146
	Dart e;
Basile Sauvage's avatar
Basile Sauvage committed
147 148
	while ( !front.empty() )
	{
Basile Sauvage's avatar
Basile Sauvage committed
149
		e = front.begin()->second;
Basile Sauvage's avatar
Basile Sauvage committed
150
		float d = front.begin()->first;
151 152

		collectVertexFromFront(e);
Basile Sauvage's avatar
Basile Sauvage committed
153

Pierre Kraemer's avatar
Pierre Kraemer committed
154
		Traversor2VVaE<MAP> tv (map, e);
Basile Sauvage's avatar
Basile Sauvage committed
155 156 157
		for (Dart f = tv.begin(); f != tv.end(); f=tv.next())
		{
			if (vmReached.isMarked(f))
158 159 160 161
			{ // f has been reached
				if (vertexInfo[f].valid) // f is in the front : update
					updateVertexInFront(f,d);
				else // f is not in the front any more (already collected) : detect a border edge
Pierre Kraemer's avatar
Pierre Kraemer committed
162 163 164 165
				{
					if ( regions[f] != regions[e] )
						border.push_back(f);
				}
Basile Sauvage's avatar
Basile Sauvage committed
166 167
			}
			else
168 169
			{ // f has not been reached : add it to the front
				addVertexToFront(f,d);
Basile Sauvage's avatar
Basile Sauvage committed
170 171 172
			}
		}
	}
Basile Sauvage's avatar
Basile Sauvage committed
173 174 175 176 177 178 179 180 181 182 183 184
	return e;
}

template <typename PFP>
void VoronoiDiagram<PFP>::computeDiagram_incremental (unsigned int nseeds)
{
	seeds.clear();

	// first seed
	srand ( time(NULL) );
	unsigned int s = rand() % map.getNbCells(VERTEX);
	unsigned int n = 0;
Pierre Kraemer's avatar
Pierre Kraemer committed
185
	TraversorV<MAP> tv (map);
Basile Sauvage's avatar
Basile Sauvage committed
186
	Dart dit = tv.begin();
Pierre Kraemer's avatar
Pierre Kraemer committed
187
	while(n < s)
Basile Sauvage's avatar
Basile Sauvage committed
188 189 190 191 192 193 194 195 196
	{
		dit = tv.next();
		++n;
	}
	seeds.push_back(dit);

	// add other seeds one by one
	Dart e = computeDiagram();

Pierre Kraemer's avatar
Pierre Kraemer committed
197
	for(unsigned int i = 1; i < nseeds ; i++)
Basile Sauvage's avatar
Basile Sauvage committed
198 199 200 201
	{
		seeds.push_back(e);
		e = computeDiagram();
	}
Basile Sauvage's avatar
Basile Sauvage committed
202 203
}

204 205
template <typename PFP>
void VoronoiDiagram<PFP>::computeDistancesWithinRegion (Dart seed)
Basile Sauvage's avatar
Basile Sauvage committed
206
{
207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223
	// init
	front.clear();
	vmReached.unmarkAll();

	vmReached.mark(seed);
	vertexInfo[seed].it = front.insert(std::pair<float,Dart>(0.0, seed));
	vertexInfo[seed].valid = true;
	vertexInfo[seed].pathOrigin = seed;

	//compute
	while ( !front.empty() )
	{
		Dart e = front.begin()->second;
		float d = front.begin()->first;

		collectVertexFromFront(e);

Pierre Kraemer's avatar
Pierre Kraemer committed
224
		Traversor2VVaE<MAP> tv (map, e);
225 226 227 228 229 230 231 232 233 234 235 236 237 238
		for (Dart f = tv.begin(); f != tv.end(); f=tv.next())
		{
			if (vmReached.isMarked(f))
			{ // f has been reached
				if (vertexInfo[f].valid) updateVertexInFront(f,d); // f is in the front : update
			}
			else
			{ // f has not been reached : add it to the front
				if ( regions[f] == regions[e] ) addVertexToFront(f,d);
			}
		}
	}
}

239 240 241 242 243
/***********************************************************
 * class CentroidalVoronoiDiagram
 ***********************************************************/

template <typename PFP>
Pierre Kraemer's avatar
Pierre Kraemer committed
244 245
CentroidalVoronoiDiagram<PFP>::CentroidalVoronoiDiagram (
	MAP& m,
246 247 248 249 250
	const EdgeAttribute<REAL, MAP>& c,
	VertexAttribute<unsigned int, MAP>& r,
	VertexAttribute<REAL, MAP>& d,
	VertexAttribute<Dart, MAP>& o,
	VertexAttribute<REAL, MAP>& a) :
Pierre Kraemer's avatar
Pierre Kraemer committed
251
		VoronoiDiagram<PFP>(m,c,r), distances(d), pathOrigins(o), areaElts(a)
252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267
{
}

template <typename PFP>
CentroidalVoronoiDiagram<PFP>::~CentroidalVoronoiDiagram ()
{
}

template <typename PFP>
void CentroidalVoronoiDiagram<PFP>::clear ()
{
	VoronoiDiagram<PFP>::clear();
	distances.setAllValues(0.0);
}

template <typename PFP>
Pierre Kraemer's avatar
Pierre Kraemer committed
268 269
void CentroidalVoronoiDiagram<PFP>::collectVertexFromFront(Dart e)
{
270
	distances[e] = this->vertexInfo[e].it->first;
271 272 273 274 275
	pathOrigins[e] = this->vertexInfo[e].pathOrigin;

	VoronoiDiagram<PFP>::collectVertexFromFront(e);
}

276

Basile Sauvage's avatar
Basile Sauvage committed
277
template <typename PFP>
Basile Sauvage's avatar
Basile Sauvage committed
278 279 280 281 282 283 284 285
void CentroidalVoronoiDiagram<PFP>::setSeeds_fromVector (const std::vector<Dart>& s)
{
	VoronoiDiagram<PFP>::setSeeds_fromVector (s);
	energyGrad.resize(this->seeds.size());
}

template <typename PFP>
void CentroidalVoronoiDiagram<PFP>::setSeeds_random (unsigned int nseeds)
Basile Sauvage's avatar
Basile Sauvage committed
286
{
Basile Sauvage's avatar
Basile Sauvage committed
287
	VoronoiDiagram<PFP>::setSeeds_random (nseeds);
Basile Sauvage's avatar
Basile Sauvage committed
288 289
	energyGrad.resize(this->seeds.size());
}
290

291
template <typename PFP>
Basile Sauvage's avatar
Basile Sauvage committed
292
void CentroidalVoronoiDiagram<PFP>::computeDiagram_incremental (unsigned int nseeds)
Basile Sauvage's avatar
Basile Sauvage committed
293
{
Basile Sauvage's avatar
Basile Sauvage committed
294
	VoronoiDiagram<PFP>::computeDiagram_incremental (nseeds);
Basile Sauvage's avatar
Basile Sauvage committed
295 296 297 298
	energyGrad.resize(this->seeds.size());
}

template <typename PFP>
Pierre Kraemer's avatar
Pierre Kraemer committed
299 300
void CentroidalVoronoiDiagram<PFP>::cumulateEnergy()
{
Basile Sauvage's avatar
Basile Sauvage committed
301
	globalEnergy = 0.0;
302 303
	for (unsigned int i = 0; i < this->seeds.size(); i++)
	{
Basile Sauvage's avatar
Basile Sauvage committed
304 305
		cumulateEnergyFromRoot(this->seeds[i]);
		globalEnergy += distances[this->seeds[i]];
306 307 308 309
	}
}

template <typename PFP>
Pierre Kraemer's avatar
Pierre Kraemer committed
310 311
void CentroidalVoronoiDiagram<PFP>::cumulateEnergyAndGradients()
{
Basile Sauvage's avatar
Basile Sauvage committed
312
	globalEnergy = 0.0;
Basile Sauvage's avatar
Basile Sauvage committed
313
	for (unsigned int i = 0; i < this->seeds.size(); i++)
Basile Sauvage's avatar
Basile Sauvage committed
314 315 316 317 318 319 320
	{
		cumulateEnergyAndGradientFromSeed(i);
		globalEnergy += distances[this->seeds[i]];
	}
}

template <typename PFP>
Pierre Kraemer's avatar
Pierre Kraemer committed
321 322
unsigned int CentroidalVoronoiDiagram<PFP>::moveSeedsOneEdgeNoCheck()
{
Basile Sauvage's avatar
Basile Sauvage committed
323 324
	unsigned int m = 0;
	for (unsigned int i = 0; i < this->seeds.size(); i++)
Basile Sauvage's avatar
Basile Sauvage committed
325
	{
326
		Dart oldSeed = this->seeds[i];
Basile Sauvage's avatar
Basile Sauvage committed
327 328 329 330 331 332 333 334
		Dart newSeed = selectBestNeighborFromSeed(i);

		// move the seed
		if (newSeed != oldSeed)
		{
			this->seeds[i] = newSeed;
			m++;
		}
Basile Sauvage's avatar
Basile Sauvage committed
335 336 337 338
	}
	return m;
}

Basile Sauvage's avatar
Basile Sauvage committed
339
template <typename PFP>
Pierre Kraemer's avatar
Pierre Kraemer committed
340 341
unsigned int CentroidalVoronoiDiagram<PFP>::moveSeedsOneEdgeCheck()
{
Basile Sauvage's avatar
Basile Sauvage committed
342 343 344 345 346 347 348 349 350 351 352 353 354 355 356
	unsigned int m = 0;
	for (unsigned int i = 0; i < this->seeds.size(); i++)
	{
		Dart oldSeed = this->seeds[i];
		Dart newSeed = selectBestNeighborFromSeed(i);

		// move the seed
		if (newSeed != oldSeed)
		{
			REAL regionEnergy = distances[oldSeed];
			this->seeds[i] = newSeed;
			this->computeDistancesWithinRegion(newSeed);
			cumulateEnergyAndGradientFromSeed(i);
			if (distances[newSeed] < regionEnergy)
				m++;
357
			else
Basile Sauvage's avatar
Basile Sauvage committed
358
				this->seeds[i] = oldSeed;
359 360 361 362
		}
	}
	return m;
}
Basile Sauvage's avatar
Basile Sauvage committed
363

Basile Sauvage's avatar
Basile Sauvage committed
364
template <typename PFP>
Pierre Kraemer's avatar
Pierre Kraemer committed
365 366
unsigned int CentroidalVoronoiDiagram<PFP>::moveSeedsToMedioid()
{
Basile Sauvage's avatar
Basile Sauvage committed
367
	unsigned int m = 0;
Basile Sauvage's avatar
Basile Sauvage committed
368
	for (unsigned int i = 0; i < this->seeds.size(); i++)
Basile Sauvage's avatar
Basile Sauvage committed
369
	{
Basile Sauvage's avatar
Basile Sauvage committed
370
		Dart oldSeed, newSeed;
Basile Sauvage's avatar
Basile Sauvage committed
371
		unsigned int seedMoved = 0;
Basile Sauvage's avatar
Basile Sauvage committed
372
		REAL regionEnergy;
Basile Sauvage's avatar
Basile Sauvage committed
373

Basile Sauvage's avatar
Basile Sauvage committed
374
		do
Basile Sauvage's avatar
Basile Sauvage committed
375
		{
Basile Sauvage's avatar
Basile Sauvage committed
376 377 378 379 380 381 382
			oldSeed = this->seeds[i];
			regionEnergy = distances[oldSeed];

			newSeed = selectBestNeighborFromSeed(i);
			this->seeds[i] = newSeed;
			this->computeDistancesWithinRegion(newSeed);
			cumulateEnergyAndGradientFromSeed(i);
Basile Sauvage's avatar
Basile Sauvage committed
383
			if (distances[newSeed] < regionEnergy)
Basile Sauvage's avatar
Basile Sauvage committed
384 385
				seedMoved = 1;
			else
Basile Sauvage's avatar
Basile Sauvage committed
386
			{
Basile Sauvage's avatar
Basile Sauvage committed
387 388
				this->seeds[i] = oldSeed;
				newSeed = oldSeed;
Basile Sauvage's avatar
Basile Sauvage committed
389
			}
Basile Sauvage's avatar
Basile Sauvage committed
390 391 392 393 394

		} while (newSeed != oldSeed);

		m += seedMoved;
	}
Basile Sauvage's avatar
Basile Sauvage committed
395 396
	return m;
}
397

Basile Sauvage's avatar
Basile Sauvage committed
398
template <typename PFP>
Pierre Kraemer's avatar
Pierre Kraemer committed
399 400
typename PFP::REAL CentroidalVoronoiDiagram<PFP>::cumulateEnergyFromRoot(Dart e)
{
Basile Sauvage's avatar
Basile Sauvage committed
401 402
	REAL distArea = areaElts[e] * distances[e];
	distances[e] = areaElts[e] * distances[e] * distances[e];
Basile Sauvage's avatar
Basile Sauvage committed
403

Pierre Kraemer's avatar
Pierre Kraemer committed
404
	Traversor2VVaE<MAP> tv (this->map, e);
405 406 407 408
	for (Dart f = tv.begin(); f != tv.end(); f=tv.next())
	{
		if ( pathOrigins[f] == this->map.phi2(f))
		{
Basile Sauvage's avatar
Basile Sauvage committed
409
			distArea += cumulateEnergyFromRoot(f);
410 411 412
			distances[e] += distances[f];
		}
	}
Basile Sauvage's avatar
Basile Sauvage committed
413
	return distArea;
Basile Sauvage's avatar
Basile Sauvage committed
414 415
}

Basile Sauvage's avatar
Basile Sauvage committed
416
template <typename PFP>
Pierre Kraemer's avatar
Pierre Kraemer committed
417 418
void CentroidalVoronoiDiagram<PFP>::cumulateEnergyAndGradientFromSeed(unsigned int numSeed)
{
Basile Sauvage's avatar
Basile Sauvage committed
419
	// precondition : energyGrad.size() > numSeed
Basile Sauvage's avatar
Basile Sauvage committed
420 421 422 423 424 425 426 427 428 429
	Dart e = this->seeds[numSeed];

	std::vector<Dart> v;
	v.reserve(8);

	std::vector<float> da;
	da.reserve(8);

	distances[e] = 0.0;

Pierre Kraemer's avatar
Pierre Kraemer committed
430
	Traversor2VVaE<MAP> tv (this->map, e);
Basile Sauvage's avatar
Basile Sauvage committed
431 432 433 434 435 436 437 438 439 440 441 442 443 444
	for (Dart f = tv.begin(); f != tv.end(); f=tv.next())
	{
		if ( pathOrigins[f] == this->map.phi2(f))
		{
			float distArea = cumulateEnergyFromRoot(f);
			da.push_back(distArea);
			distances[e] += distances[f];
			v.push_back(f);
		}
	}

	// compute the gradient
	// TODO : check if the computation of grad and proj is still valid for other edgeCost than geodesic distances
	VEC3 grad (0.0);
Sauvage's avatar
Sauvage committed
445
	const VertexAttribute<VEC3, MAP>& pos = this->map.template getAttribute<VEC3, VERTEX, typename PFP::MAP>("position");
Basile Sauvage's avatar
Basile Sauvage committed
446

Pierre Kraemer's avatar
Pierre Kraemer committed
447
	for (unsigned int j = 0; j < v.size(); ++j)
Basile Sauvage's avatar
Basile Sauvage committed
448 449 450 451 452 453 454 455 456 457 458 459 460 461 462
	{
		Dart f = v[j];
		VEC3 edgeV = pos[f] - pos[this->map.phi2(f)];
		edgeV.normalize();
		grad += da[j] * edgeV;
	}
	grad /= 2.0;
	energyGrad[numSeed] = grad;
}

template <typename PFP>
Dart CentroidalVoronoiDiagram<PFP>::selectBestNeighborFromSeed(unsigned int numSeed)
{
	Dart e = this->seeds[numSeed];
	Dart newSeed = e;
Sauvage's avatar
Sauvage committed
463
	const VertexAttribute<VEC3, MAP>& pos = this->map.template getAttribute<VEC3,VERTEX,typename PFP::MAP>("position");
Basile Sauvage's avatar
Basile Sauvage committed
464 465 466

	// TODO : check if the computation of grad and proj is still valid for other edgeCost than geodesic distances
	float maxProj = 0.0;
Pierre Kraemer's avatar
Pierre Kraemer committed
467
	Traversor2VVaE<MAP> tv (this->map, e);
Basile Sauvage's avatar
Basile Sauvage committed
468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508
	for (Dart f = tv.begin(); f != tv.end(); f=tv.next())
	{
		if ( pathOrigins[f] == this->map.phi2(f))
		{
			VEC3 edgeV = pos[f] - pos[this->map.phi2(f)];
	//		edgeV.normalize();
			float proj = edgeV * energyGrad[numSeed];
			if (proj > maxProj)
			{
				maxProj = proj;
				newSeed = f;
			}
		}
	}
	return newSeed;
}

/*
template <typename PFP>
unsigned int CentroidalVoronoiDiagram<PFP>::moveSeed(unsigned int numSeed){

	// collect energy and compute the gradient
//	cumulateEnergyAndGradientFromSeed(numSeed);

	// select the new seed
	Dart newSeed = selectBestNeighborFromSeed(numSeed);

	// move the seed
	Dart oldSeed = this->seeds[numSeed];
	unsigned int seedMoved = 0;
	if (newSeed != oldSeed)
	{
		this->seeds[numSeed] = newSeed;
		seedMoved = 1;
	}

	return seedMoved;
}
*/

/*
Basile Sauvage's avatar
Basile Sauvage committed
509
template <typename PFP>
Basile Sauvage's avatar
Basile Sauvage committed
510 511 512
unsigned int CentroidalVoronoiDiagram<PFP>::moveSeed(unsigned int numSeed){
	Dart e = this->seeds[numSeed];
	unsigned int seedMoved = 0;
Basile Sauvage's avatar
Basile Sauvage committed
513 514 515 516

	std::vector<Dart> v;
	v.reserve(8);

Basile Sauvage's avatar
Basile Sauvage committed
517 518
	std::vector<float> da;
	da.reserve(8);
Basile Sauvage's avatar
Basile Sauvage committed
519 520 521 522 523 524 525 526

	distances[e] = 0.0;

	Traversor2VVaE<typename PFP::MAP> tv (this->map, e);
	for (Dart f = tv.begin(); f != tv.end(); f=tv.next())
	{
		if ( pathOrigins[f] == this->map.phi2(f))
		{
Basile Sauvage's avatar
Basile Sauvage committed
527 528
			float distArea = cumulateEnergyFromRoot(f);
			da.push_back(distArea);
Basile Sauvage's avatar
Basile Sauvage committed
529 530 531 532 533
			distances[e] += distances[f];
			v.push_back(f);
		}
	}

Basile Sauvage's avatar
Basile Sauvage committed
534
	// TODO : check if the computation of grad and proj is still valid for other edgeCost than geodesic distances
Basile Sauvage's avatar
Basile Sauvage committed
535 536 537
	VEC3 grad (0.0);
	const VertexAttribute<VEC3>& pos = this->map.template getAttribute<VEC3,VERTEX>("position");

Basile Sauvage's avatar
Basile Sauvage committed
538
	// compute the gradient
Basile Sauvage's avatar
Basile Sauvage committed
539 540 541 542 543
	for (unsigned int j = 0; j<v.size(); ++j)
	{
		Dart f = v[j];
		VEC3 edgeV = pos[f] - pos[this->map.phi2(f)];
		edgeV.normalize();
Basile Sauvage's avatar
Basile Sauvage committed
544
		grad += da[j] * edgeV;
Basile Sauvage's avatar
Basile Sauvage committed
545
	}
Basile Sauvage's avatar
Basile Sauvage committed
546
	grad /= 2.0;
Basile Sauvage's avatar
Basile Sauvage committed
547 548

	float maxProj = 0.0;
Basile Sauvage's avatar
Basile Sauvage committed
549
//	float memoForTest = 0.0;
Basile Sauvage's avatar
Basile Sauvage committed
550 551 552 553 554 555
	for (unsigned int j = 0; j<v.size(); ++j)
	{
		Dart f = v[j];
		VEC3 edgeV = pos[f] - pos[this->map.phi2(f)];
//		edgeV.normalize();
		float proj = edgeV * grad;
Basile Sauvage's avatar
Basile Sauvage committed
556
//		proj -= areaElts[e] * this->edgeCost[f] * this->edgeCost[f];
Basile Sauvage's avatar
Basile Sauvage committed
557 558
		if (proj > maxProj)
		{
Basile Sauvage's avatar
Basile Sauvage committed
559 560
//			if (numSeed==1) memoForTest = (edgeV * grad) / (areaElts[e] * this->edgeCost[f] * this->edgeCost[f]);
//				CGoGNout << (edgeV * grad) / (areaElts[e] * this->edgeCost[f] * this->edgeCost[f]) * this->seeds.size() << CGoGNendl;
Basile Sauvage's avatar
Basile Sauvage committed
561
//				CGoGNout << edgeV * grad << "\t - \t" << areaElts[e] * this->edgeCost[f] * this->edgeCost[f] << CGoGNendl;
Basile Sauvage's avatar
Basile Sauvage committed
562 563 564 565 566
			maxProj = proj;
			seedMoved = 1;
			this->seeds[numSeed] = v[j];
		}
	}
Basile Sauvage's avatar
Basile Sauvage committed
567

Basile Sauvage's avatar
Basile Sauvage committed
568
	return seedMoved;
569
}
Basile Sauvage's avatar
Basile Sauvage committed
570
*/
Pierre Kraemer's avatar
Pierre Kraemer committed
571 572 573 574 575 576 577 578

} // namespace Geometry

} // namespace Surface

} // namespace Algo

} // namespace CGoGN