26 |
//#include "utils.hpp" |
//#include "utils.hpp" |
27 |
#endif /* STANDALONE */ |
#endif /* STANDALONE */ |
28 |
|
|
29 |
|
#include <cstdlib> |
30 |
|
|
31 |
|
int myrand(const int width) { |
32 |
|
return (int)(static_cast<double>(width)*rand()/(RAND_MAX+1.0)); |
33 |
|
} |
34 |
|
|
35 |
|
#include <map> |
36 |
|
#include <set> |
37 |
/************************************************************/ |
/************************************************************/ |
38 |
/******************** Fonctions locales *********************/ |
/******************** Fonctions locales *********************/ |
39 |
/************************************************************/ |
/************************************************************/ |
48 |
// return Vertex((tr[1]-tr[0]) ^ (tr[2]-tr[0]))[2]; |
// return Vertex((tr[1]-tr[0]) ^ (tr[2]-tr[0]))[2]; |
49 |
// } |
// } |
50 |
|
|
51 |
inline bool checkInter(const Vertex &P1, const Vertex &P2, const Vertex &Q1, const Vertex &Q2) { |
inline bool checkInter(const Vertex &P1, const Vertex &P2, const Vertex &Q1, const Vertex &Q2) { |
52 |
// retourne true si P1P2 inter Q1Q2 != vide |
// retourne true si P1P2 inter Q1Q2 != vide |
53 |
assert(!(P1 == P2 || P1 == Q1 || P1 == Q2 || P2 == Q1 || P2 == Q2 || Q1 == Q2)); |
assert(!(P1 == P2 || P1 == Q1 || P1 == Q2 || P2 == Q1 || P2 == Q2 || Q1 == Q2)); |
54 |
|
return (SDet(P1,P2,Q1)*SDet(P1,P2,Q2) < 0) && (SDet(Q1,Q2,P1)*SDet(Q1,Q2,P2) < 0); // version optimisée |
|
// const double |
|
|
// s1 = SDet(P1,P2,Q1)*SDet(P1,P2,Q2), |
|
|
// s2 = SDet(Q1,Q2,P1)*SDet(Q1,Q2,P2); |
|
|
// if (s1 == 0 || s2 == 0) ffout(4) << "warning s=0\n"; |
|
|
// return (s1 <= 0 && s2 <= 0); |
|
|
|
|
|
const double eps = 1e-15; |
|
|
if (std::abs(SDet(P1,P2,Q1)) < eps || |
|
|
std::abs(SDet(P1,P2,Q2)) < eps || |
|
|
std::abs(SDet(Q1,Q2,P1)) < eps || |
|
|
std::abs(SDet(Q1,Q2,P2)) < eps) { |
|
|
//ffout(0) << " !!!!! NO !!!!!!!" << std::endl; |
|
|
return false; |
|
|
} |
|
|
|
|
|
return (SDet(P1,P2,Q1)*SDet(P1,P2,Q2) < 0) && (SDet(Q1,Q2,P1)*SDet(Q1,Q2,P2) < 0); // version optimisée |
|
55 |
} |
} |
56 |
|
|
57 |
inline double Angle(const Vertex &V1, const Vertex &V2, const Vertex &V3) { // cos de l'angle * ||V2-V1|| (V1 point pivot) |
inline double Angle(const Vertex &V1, const Vertex &V2, const Vertex &V3) { // cos de l'angle * ||V2-V1|| (V1 point pivot) |
135 |
|
|
136 |
#endif /* ALTERNATE_INCIRCLE */ |
#endif /* ALTERNATE_INCIRCLE */ |
137 |
|
|
138 |
void Triangulation::permutation(TriangleIndex num_tri1, const unsigned k) { |
void Triangulation::permutation(TriangleIndex num_tri1, const unsigned k1) { |
139 |
// realise la permutation d'une arête commune a deux triangles. |
// realise la permutation d'une arête commune a deux triangles. |
140 |
// TriangleIndex = pointeur du triangle. |
// TriangleIndex = pointeur du triangle. |
141 |
// k = numero du triangle voisin. |
// k1 = numero du triangle voisin. |
142 |
// S = numero du sommet |
// S = numero du sommet |
143 |
// V = numero des triangles voisins. |
// V = numero des triangles voisins. |
144 |
|
|
145 |
// S1 S1 |
// S1 S1 |
146 |
// * * |
// * * |
147 |
// V2 * * V1 V2 *** V1 |
// V2 * * V1 V2 *** V1 |
148 |
// * * * * * |
// * * * * * |
149 |
// * T1 * * * * |
// * T1 * * * * |
150 |
// S2*********S4 S2 * T1*T2 *S4 |
// S2*********S4 S2 * T1*T2 *S4 |
151 |
// * T2 * * * * |
// * T2 * * * * |
152 |
// * * * * * |
// * * * * * |
153 |
// V3 * * V4 V3 *** V4 |
// V3 * * V4 V3 *** V4 |
154 |
// * * |
// * * |
155 |
// S3 S3 |
// S3 S3 |
156 |
|
|
157 |
// => modification de V1 qui deivient T2 et V3 qui devient T1 |
// => modification de V1 qui deivient T2 et V3 qui devient T1 |
158 |
|
|
159 |
|
const TriangleIndex num_tri2 = num_tri1->neigh(k1); |
160 |
|
assert(num_tri2 != NULL); |
161 |
|
const unsigned k2 = num_tri2->localizeNeigh(num_tri1); |
162 |
|
|
163 |
|
const PointIndex |
164 |
|
S1=num_tri1->base()[k1], |
165 |
|
S2=num_tri1->base()[(k1+1)%3], |
166 |
|
S4=num_tri1->base()[(k1+2)%3]; |
167 |
|
|
168 |
|
const unsigned tag01 = num_tri1->tag; |
169 |
|
const unsigned tag02 = num_tri2->tag; |
170 |
|
// say() << *num_tri1 << " " << k1 << " " << tag01 << " " << (tag01 & HSide[k1]) << "\n"; |
171 |
|
// say() << *num_tri2 << " " << k2 << " " << tag02 << " " << (tag02 & HSide[k2]) << "\n"; |
172 |
|
|
173 |
|
assert((tag01 & HSide[k1]) == 0); // Ne flip pas des arêtes marquées |
174 |
|
assert((tag02 & HSide[k2]) == 0); |
175 |
|
num_tri1->tag = ((tag01&HSide[(k1+2)%3])?HSide[2]:0)|((tag02&HSide[(k2+1)%3])?HSide[0]:0); |
176 |
|
num_tri2->tag = ((tag01&HSide[(k1+1)%3])?HSide[1]:0)|((tag02&HSide[(k2+2)%3])?HSide[0]:0); |
177 |
|
|
178 |
|
const TriangleIndex |
179 |
|
V1=num_tri1->neigh((k1+1)%3), |
180 |
|
V2=num_tri1->neigh((k1+2)%3); |
181 |
|
|
182 |
|
PointIndex S3; |
183 |
|
TriangleIndex V3,V4; |
184 |
|
|
185 |
|
// recherche du sommet dans le triangle num_tri2 opposé au coté commun (ie opposé à k1 dans num_tri1) |
186 |
|
for(unsigned k=0;k<3;++k) |
187 |
|
if(num_tri2->neigh(k) == num_tri1) { |
188 |
|
S3=num_tri2->base()[k]; |
189 |
|
V3=num_tri2->neigh((k+1)%3); // Ici, un bug est mort dans d'attroces souffrances |
190 |
|
V4=num_tri2->neigh((k+2)%3); |
191 |
|
break; |
192 |
|
} |
193 |
|
|
194 |
const TriangleIndex num_tri2 = num_tri1->neigh(k); |
for(unsigned k=0;k<3;++k) { |
195 |
assert(num_tri2 != NULL); |
if(V1 != NULL) |
196 |
|
if(V1->neigh(k) == num_tri1) |
197 |
PointIndex |
V1->neigh(k) = num_tri2; |
198 |
S1=num_tri1->base()[k], |
if(V3 != NULL) |
199 |
S2=num_tri1->base()[(k+1)%3], |
if(V3->neigh(k) == num_tri2) |
200 |
S4=num_tri1->base()[(k+2)%3]; |
V3->neigh(k) = num_tri1; |
201 |
|
} |
|
TriangleIndex |
|
|
V1=num_tri1->neigh((k+1)%3), |
|
|
V2=num_tri1->neigh((k+2)%3); |
|
|
|
|
|
PointIndex S3; |
|
|
TriangleIndex V3,V4; |
|
202 |
|
|
203 |
// recherche du sommet dans le triangle num_tri2 opposé au coté commun (ie opposé à k dans num_tri1) |
num_tri1->setVertices(S1,S2,S3); |
204 |
for(unsigned k=0;k<3;++k) |
num_tri2->setVertices(S1,S3,S4); |
205 |
if(num_tri2->neigh(k) == num_tri1) { |
num_tri1->setNeighs(V3,num_tri2,V2); |
206 |
S3=num_tri2->base()[k]; |
num_tri2->setNeighs(V4,V1,num_tri1); |
|
V3=num_tri2->neigh((k+1)%3); // Ici, un bug est mort dans d'attroces souffrances |
|
|
V4=num_tri2->neigh((k+2)%3); |
|
|
break; |
|
|
} |
|
|
|
|
|
for(unsigned k=0;k<3;++k) { |
|
|
if(V1 != NULL) |
|
|
if(V1->neigh(k) == num_tri1) |
|
|
V1->neigh(k) = num_tri2; |
|
|
if(V3 != NULL) |
|
|
if(V3->neigh(k) == num_tri2) |
|
|
V3->neigh(k) = num_tri1; |
|
|
} |
|
|
|
|
|
num_tri1->setVertices(S1,S2,S3); |
|
|
num_tri2->setVertices(S1,S3,S4); |
|
|
num_tri1->setNeighs(V3,num_tri2,V2); |
|
|
num_tri2->setNeighs(V4,V1,num_tri1); |
|
207 |
} |
} |
208 |
|
|
209 |
|
|
210 |
void Triangulation::elimination_tri(TriangleIndex num_tri) { |
void Triangulation::elimination_tri(TriangleIndex num_tri) { |
211 |
// elimination fu triangle pointe par neigh_index |
// elimination fu triangle pointe par neigh_index |
356 |
i->tag = 0; |
i->tag = 0; |
357 |
} |
} |
358 |
|
|
359 |
|
|
360 |
bool Triangulation::checkLine(const CurveVertex &L, const bool closed) { |
bool Triangulation::checkLine(const CurveVertex &L, const bool closed) { |
361 |
assert(L.size() > 1); |
assert(L.size() > 1); |
362 |
if (closed) |
assert(closed); // !closed not implemented |
363 |
assert(L.front() == L.back()); |
assert(!closed || L.front() == L.back()); |
364 |
|
|
365 |
|
typedef TinyVector<2,PointIndex> Edge; |
366 |
|
typedef std::map<Edge,TinyVector<2,TriangleIndex> > EdgeMapping; |
367 |
|
|
368 |
|
typedef std::set<Edge> BorderEdges; |
369 |
|
BorderEdges borderEdges; |
370 |
|
PointIndex lastPoint = L.front(); |
371 |
|
|
372 |
|
for(CurveVertex::const_iterator i=L.begin(); ++i != L.end();) { |
373 |
|
Edge edge; |
374 |
|
edge[0]=lastPoint; |
375 |
|
edge[1]=*i; |
376 |
|
if(edge[0] > edge[1]) std::swap(edge[0],edge[1]); // arithmétique des pointeurs ou entiers |
377 |
|
borderEdges.insert(edge); |
378 |
|
lastPoint = *i; |
379 |
|
} |
380 |
|
|
381 |
bool change = false; |
EdgeMapping edgeMapping; |
382 |
|
do { |
383 |
|
edgeMapping.clear(); |
384 |
|
for(Triangles::iterator tr = triangles.begin();tr !=triangles.end(); ++tr) { |
385 |
|
for(unsigned j=0; j<3; ++j) { |
386 |
|
PointIndex p1 = tr->base()[(j+1)%3]; |
387 |
|
PointIndex p2 = tr->base()[(j+2)%3]; |
388 |
|
if(p1 > p2) std::swap(p1,p2); // arithmétique des pointeurs ou des entiers |
389 |
|
Edge edge; |
390 |
|
edge[0]=p1; |
391 |
|
edge[1]=p2; |
392 |
|
const TriangleIndex tr2 = tr->neigh(j); |
393 |
|
|
394 |
|
if ((tr2 != NULL) && !(tr->tag & HSide[j]) && edgeMapping.find(edge) == edgeMapping.end()) { |
395 |
|
Edge edgetemp; |
396 |
|
edgetemp[0]=p1; |
397 |
|
edgetemp[1]=p2; |
398 |
|
TinyVector<2,TriangleIndex> trtemp; |
399 |
|
trtemp[0]=&*tr; |
400 |
|
trtemp[1]=tr2; |
401 |
|
edgeMapping.insert(EdgeMapping::value_type(edgetemp,trtemp)); |
402 |
|
} |
403 |
|
} |
404 |
|
} |
405 |
|
|
406 |
CurveVertex::const_iterator position = L.begin(); |
for(BorderEdges::const_iterator i = borderEdges.begin(); i != borderEdges.end();) { |
407 |
const CurveVertex::const_iterator end = L.end(); |
EdgeMapping::iterator foundi = edgeMapping.find(*i); |
408 |
|
if (foundi != edgeMapping.end()) { |
409 |
PointIndex |
// say() << "Edge " << *(*i)[0] << " " << *(*i)[1] << " found\n"; |
410 |
curVertex = *position++, |
TriangleIndex tr0 = foundi->second[0]; |
411 |
curTarget = *position; |
TriangleIndex tr1 = foundi->second[1]; |
412 |
|
assert(tr0 != NULL && tr1 != NULL); |
413 |
// Etape initiale : recherche d'un triangle possédant en sommet le premier point de L |
tr0->tag |= HSide[tr0->localizeNeigh(tr1)]; |
414 |
// On peut y mettre une recherche quadtree ou par adjacente optimale (cf find_P_in_elt) |
tr1->tag |= HSide[tr1->localizeNeigh(tr0)]; |
415 |
TriangleIndex curTriangle = findVertex(curVertex); |
edgeMapping.erase(foundi); |
416 |
|
BorderEdges::const_iterator i_to_rm = i++; |
417 |
bool found = false; |
borderEdges.erase(i_to_rm); |
418 |
while (!found) { |
} else ++i; |
419 |
assert(curTriangle->base()[0] == curVertex || |
} |
420 |
curTriangle->base()[1] == curVertex || |
|
421 |
curTriangle->base()[2] == curVertex); |
if (!borderEdges.empty()) { |
422 |
|
bool done = false; |
423 |
unsigned localVertex=curTriangle->localizeVertex(curVertex); // position local de curVertex |
while (!done) { |
424 |
|
EdgeMapping::iterator e = edgeMapping.begin(); |
425 |
if (curTriangle->base()[(localVertex+1)%3] == curTarget) { |
for(unsigned n = myrand(edgeMapping.size());n>0;--n) ++e; |
426 |
// On a trouvé curTarget dans le triangle courant, on passe au suivant. |
const PointIndex P1 = e->first[0]; |
427 |
// On marque le coté comme frontière si closed == true ie c'est pas une ligne interne |
const PointIndex P2 = e->first[1]; |
428 |
if (closed) { |
TriangleIndex tr0 = e->second[0]; |
429 |
// coté entre localVertex et LocalVertex+1 => neigh[localVertex+2] |
TriangleIndex tr1 = e->second[1]; |
430 |
curTriangle->tag |= HSide[(localVertex+2)%3]; |
assert(tr0 != NULL && tr1 != NULL); |
431 |
TriangleIndex NTri = curTriangle->neigh((localVertex+2)%3); // forcement non null |
const unsigned l0 = tr0->localizeNeigh(tr1); |
432 |
assert(NTri != NULL); // Car tout est inclus dans la boite |
const unsigned l1 = tr1->localizeNeigh(tr0); |
433 |
NTri->tag |= HSide[NTri->localizeNeigh(curTriangle)]; |
const PointIndex Q1 = tr0->base()[l0]; |
434 |
} |
const PointIndex Q2 = tr1->base()[l1]; |
435 |
curVertex = curTarget; |
if (checkInter(*P1,*P2,*Q1,*Q2)) { |
436 |
found = (++position == end); |
// say() << "Flip edge " << *P1 << " " << *P2 << "\n"; |
437 |
curTarget = *position; |
done = true; |
438 |
|
permutation(tr0,tr0->localizeNeigh(tr1)); |
439 |
} else if (curTriangle->base()[(localVertex+2)%3] == curTarget) { |
} |
|
// On a trouvé curTarget dans le triangle courant, on passe au suivant. |
|
|
// On marque le coté comme frontière si closed == true ie c'est pas une ligne interne |
|
|
if (closed) { |
|
|
// coté entre localVertex et LocalVertex+1 => neigh[localVertex+2] |
|
|
curTriangle->tag |= HSide[(localVertex+1)%3]; |
|
|
TriangleIndex NTri = curTriangle->neigh((localVertex+1)%3); // forcement non null |
|
|
assert(NTri != NULL); // Car tout est inclus dans la boite |
|
|
NTri->tag |= HSide[NTri->localizeNeigh(curTriangle)]; |
|
|
} |
|
|
curVertex = curTarget; |
|
|
found = (++position == end); |
|
|
curTarget = *position; |
|
|
|
|
|
} else if(checkInter(*curVertex,*curTarget, |
|
|
(*curTriangle)((localVertex+1)%3), |
|
|
(*curTriangle)((localVertex+2)%3))) { |
|
|
// Ok : found a problem on the border |
|
|
// gestion spécifique du cnflit arète-frontière |
|
|
bool atteint = false; |
|
|
while(!atteint) { |
|
|
assert(curTriangle != NULL); |
|
|
|
|
|
// description globale du coté opposée à localvertex dans curTriangle ie l'entrée |
|
|
const PointIndex |
|
|
g1 = curTriangle->base()[(localVertex+1)%3], |
|
|
g2 = curTriangle->base()[(localVertex+2)%3]; |
|
|
|
|
|
TriangleIndex new_t = curTriangle->neigh(localVertex); |
|
|
assert(new_t != NULL); |
|
|
|
|
|
unsigned l=new_t->localizeNeigh(curTriangle); |
|
|
|
|
|
const PointIndex from = curTriangle->base()[localVertex]; |
|
|
const PointIndex to = new_t->base()[l]; |
|
|
|
|
|
// Test de convexité |
|
|
if (checkInter(*from,*to,*g1,*g2)) { // c'est convexe ou fusinnable |
|
|
#ifndef STEP_DETAILS |
|
|
// Test de fusion |
|
|
// g1 sur [from,to] |
|
|
if (SDet(*from,*to,*g1) == 0) |
|
|
ffout(4) << "G1 is on the border\n"; |
|
|
// g2 sur [from,to] |
|
|
if (SDet(*from,*to,*g2) == 0) |
|
|
ffout(4) << "G2 is on the border\n"; |
|
|
// from sur [g1,g2] |
|
|
if (SDet(*g1,*g2,*from) == 0) |
|
|
ffout(4) << "from is on the border\n"; |
|
|
// tp sur [g1,g2] |
|
|
if (SDet(*g1,*g2,*to) == 0) |
|
|
ffout(4) << "to is on the border\n"; |
|
|
#endif /* STEP_DETAILS */ |
|
|
|
|
|
permutation(curTriangle,localVertex); |
|
|
change = true; // changement |
|
|
|
|
|
// le cas où g3 serait curTarget |
|
|
localVertex=curTriangle->localizeVertex(from); |
|
|
l=new_t->localizeVertex(from); |
|
|
|
|
|
if (to == curTarget) { // passer au point suivant. |
|
|
#warning ! |
|
|
// if (closed && from == curVertex) { |
|
|
// // on a l'arête complète |
|
|
// curTriangle->tag |= HSide[curTriangle->localizeNeigh(new_t)]; |
|
|
// new_t->tag |= HSide[new_t->localizeNeigh(curTriangle)]; |
|
|
// } |
|
|
#warning end of ! |
|
|
curVertex = curTarget; |
|
|
found = (++position == end); |
|
|
curTarget = *position; |
|
|
atteint = true; |
|
|
} else { |
|
|
// on pourrait y mettre les noms globaux que l'on y connait |
|
|
if (checkInter(*curVertex,*curTarget, |
|
|
(*curTriangle)((localVertex+1)%3), |
|
|
(*curTriangle)((localVertex+2)%3))) { |
|
|
// curTriangle : OK |
|
|
// localVertex : OK |
|
|
} else { |
|
|
curTriangle = new_t; |
|
|
localVertex = l; |
|
|
assert(checkInter(*curVertex,*curTarget, |
|
|
(*curTriangle)((l+1)%3), |
|
|
(*curTriangle)((l+2)%3))); |
|
440 |
} |
} |
|
} |
|
|
} else { |
|
|
// c'est pas convexe et ça coupe encore l'arète |
|
|
curTriangle = new_t; |
|
|
// test de continuation |
|
|
if (checkInter(*curVertex,*curTarget, |
|
|
(*curTriangle)(l), |
|
|
(*curTriangle)((l+1)%3))) { |
|
|
// sinon c'est l'autre |
|
|
// c'est sur la droite [l, l+1] opp à l+2 |
|
|
localVertex = (l+2)%3; |
|
|
} else { |
|
|
// c'est sur la droite [l,l+2] opp à l+1 |
|
|
assert(checkInter(*curVertex,*curTarget, |
|
|
(*curTriangle)(l), |
|
|
(*curTriangle)((l+2)%3))); |
|
|
localVertex = (l+1)%3; |
|
|
} |
|
441 |
} |
} |
442 |
} |
} while (!borderEdges.empty()); |
|
} else { |
|
|
// A fortment optimiser, peut coupable avec le calcul d'intersection |
|
|
const double theta1 = Angle(*curVertex,*curTarget, |
|
|
(*curTriangle)((localVertex+1)%3)); |
|
|
const double theta2 = Angle(*curVertex,*curTarget, |
|
|
(*curTriangle)((localVertex+2)%3)); |
|
|
const int sens = (theta1 > theta2)?2:1; |
|
|
// ffout(4) << "angles : " << theta1 << " " << theta2 << " sens : " << sens << "\n"; |
|
|
curTriangle = curTriangle->neigh((localVertex+sens)%3); |
|
|
} |
|
|
} // end of while(!found) |
|
443 |
|
|
444 |
return change; |
return false; |
445 |
} |
} |
446 |
|
|
447 |
bool Triangulation::colorize(const Points & points) { |
bool Triangulation::colorize(const Points & points) { |