00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #include "g-map-vertex.hh"
00026 #include "plane.hh"
00027 #include <list>
00028 using namespace std;
00029 using namespace GMap3d;
00030
00031 CDart * CGMapVertex::triangulateEdge(CDart * ADart)
00032 {
00033 assert(ADart!=NULL);
00034
00035
00036
00037
00038
00039 return CGMapGeneric::triangulateEdge(ADart);
00040 }
00041
00042 CDart * CGMapVertex::triangulateFace(CDart * ADart)
00043 {
00044 assert(ADart!=NULL);
00045
00046 CVertex bary = barycenter(ADart, ORBIT_01);
00047 CDart * center= CGMapGeneric::triangulateFace(ADart);
00048 setVertex(center, bary);
00049
00050 return center;
00051 }
00052
00053 CDart * CGMapVertex::triangulateVolume(CDart * ADart)
00054 {
00055 assert(ADart!=NULL);
00056
00057 CVertex bary = barycenter(ADart, ORBIT_012);
00058 CDart * center= CGMapGeneric::triangulateVolume(ADart);
00059 setVertex(center, bary);
00060
00061 return center;
00062 }
00063
00064 #define VTX(d) (AVertexDI < 0 ? findVertex(d) : \
00065 (CAttributeVertex*)getDirectInfo(d, AVertexDI))
00066
00067 void CGMapVertex::triangulateGeoFace(CDart * AFace, int ANewEdgesMark,
00068 int AVertexDI)
00069 {
00070 if ( alpha101(AFace) == alpha010(AFace) )
00071 return;
00072
00073 CVertex v1, v2, v3, pt, normal;
00074 CDart *face, *d1, *d2, *d3, *iter, *best;
00075 TCoordinate coef, best_coef;
00076 int mark = getNewMark();
00077 list<CDart*> new_edges;
00078
00079 normal = faceNormalVector(AFace);
00080
00081 face = AFace;
00082 while (alpha101(face) != alpha010(face)) {
00083
00084 d1 = face;
00085 best = NULL;
00086 best_coef = 0;
00087 do {
00088 d2 = alpha01(d1);
00089 d3 = alpha10(d1);
00090
00091 v1 = *VTX(d1);
00092 v2 = *VTX(d2);
00093 v3 = *VTX(d3);
00094
00095 if (((v2 - v1) * (v3 - v1)).dot(normal) > 0.0) {
00096 markOrbit(d1, ORBIT_VERTEX, mark);
00097 markOrbit(d2, ORBIT_VERTEX, mark);
00098 markOrbit(d3, ORBIT_VERTEX, mark);
00099 iter = alpha01(d2);
00100 while (iter != d3 &&
00101 (isMarked(iter, mark) ||
00102 !isPointInTriangle(*VTX(iter), v1, v2, v3, normal)))
00103 iter = alpha01(iter);
00104 unmarkOrbit(d1, ORBIT_VERTEX, mark);
00105 unmarkOrbit(d2, ORBIT_VERTEX, mark);
00106 unmarkOrbit(d3, ORBIT_VERTEX, mark);
00107 if (iter == d3) {
00108 coef = getTriangleCoef(v1, v2, v3);
00109 if (coef > best_coef) {
00110 best_coef = coef;
00111 best = d1;
00112 }
00113 }
00114 }
00115
00116 d1 = d2;
00117 }
00118 while(d1 != face);
00119
00120 if (!best) {
00121 unmarkAll(1);
00122 setMark(face, 1);
00123 save("triangulateGeoFace.moka", BinaryFormat);
00124 }
00125
00126 assert(best != NULL);
00127
00128
00129 d1 = best;
00130 d2 = alpha01(d1);
00131 d3 = alpha10(d1);
00132 insertEdge(alpha1(d2), d3);
00133
00134 new_edges.push_back(alpha1(d3));
00135 if (ANewEdgesMark >= 0)
00136 markOrbit(alpha1(d3), ORBIT_EDGE, ANewEdgesMark);
00137 if (AVertexDI >= 0) {
00138 pointDirectInfoToAttributeVertex(AVertexDI, d2);
00139 pointDirectInfoToAttributeVertex(AVertexDI, d3);
00140 }
00141 face = d2;
00142 normal = faceNormalVector(face);
00143 }
00144
00145 while (!new_edges.empty()) {
00146 if (shouldSwapEdge(new_edges.front(), AVertexDI))
00147 swapEdge(new_edges.front(), AVertexDI);
00148 new_edges.pop_front();
00149 }
00150
00151 freeMark(mark);
00152 }
00153
00154 void CGMapVertex::triangulateMarkedFaces(int AMark,
00155 int ANewEdgesMark,
00156 int AVertexDI)
00157 {
00158 bool closed;
00159 int mark = getNewMark();
00160
00161 for (CDynamicCoverageAll dca(this); dca.cont(); ++dca)
00162 {
00163 if (isMarked(*dca, AMark) )
00164 {
00165 if (!isMarked(*dca, mark))
00166 {
00167 closed = true;
00168
00169 for (CDynamicCoverageFace dcf(this, *dca); dcf.cont(); ++dcf)
00170 {
00171 setMark(*dcf, mark);
00172 if (isFree0(*dcf) || isFree1(*dcf))
00173 closed = false;
00174 }
00175
00176 if (closed)
00177 {
00178 triangulateGeoFace(*dca, ANewEdgesMark, AVertexDI);
00179 }
00180 }
00181 }
00182 else
00183 setMark(*dca, mark);
00184 }
00185
00186 unmarkAll(mark);
00187 freeMark(mark);
00188 }
00189
00190 bool CGMapVertex::shouldSwapEdge(CDart * AEdge, int AVertexDI)
00191 {
00192 CVertex p[4];
00193 bool result = false;
00194
00195 p[0] = *VTX(AEdge);
00196 p[1] = *VTX(alpha210(AEdge));
00197 p[2] = *VTX(alpha0(AEdge));
00198 p[3] = *VTX(alpha10(AEdge));
00199
00200 CVertex v01 = p[1] - p[0];
00201 CVertex v03 = p[3] - p[0];
00202 CVertex v23 = p[3] - p[2];
00203 CVertex v21 = p[1] - p[2];
00204
00205 CVertex v013 = v01 * v03;
00206 CVertex v231 = v23 * v21;
00207
00208 CVertex n013 = CGeometry::getNormalVector(p[0],p[1], p[3]);
00209 CVertex n123 = CGeometry::getNormalVector(p[1], p[2], p[3]);
00210
00211 if (v013.dot(v231) > 0 && !n013.isNull() && !n123.isNull()) {
00212 CVertex v120 = v21 * v01;
00213 CVertex v302 = v03 * v23;
00214
00215 TCoordinate diff1 = fabs(v013.norm() - v231.norm());
00216 TCoordinate diff2 = fabs(v120.norm() - v302.norm());
00217
00218 result = diff2 > diff1;
00219 }
00220
00221 return result;
00222 }
00223
00224 bool CGMapVertex::swapEdge(CDart * AEdge, int AVertexDI)
00225 {
00226 if ( isFree2(AEdge) || isFree1(AEdge) || isFree1(alpha0(AEdge)) ||
00227 isFree1(alpha2(AEdge)) || isFree1(alpha20(AEdge)) )
00228 return false;
00229
00230 CDart *d[2] = {AEdge, alpha20(AEdge)};
00231 CDart *old[2] = {alpha1(d[0]), alpha1(d[1])};
00232 CDart *d1, *d2;
00233
00234 for (int i = 0; i < 2; i++) {
00235 d1 = alpha1(d[i]);
00236 d2 = alpha21(d[i]);
00237
00238 unsew1(d1);
00239 unsew1(d2);
00240 sew1(d1, d2);
00241
00242 d1 = alpha0(d2);
00243 d2 = alpha1(d1);
00244
00245 unsew1(d1);
00246 sew1(d[i], d1);
00247 sew1(alpha2(d[i]), d2);
00248 }
00249
00250 if (AVertexDI >= 0)
00251 for (int i = 0; i < 2; i++)
00252 {
00253 pointDirectInfoToAttributeVertex(AVertexDI, d[i]);
00254 pointDirectInfoToAttributeVertex(AVertexDI, old[i]);
00255 }
00256
00257 return true;
00258 }
00259
00260 bool CGMapVertex::isPointInTriangle(const CVertex & APoint,
00261 const CVertex & AVertex1,
00262 const CVertex & AVertex2,
00263 const CVertex & AVertex3,
00264 const CVertex & ANormal)
00265 {
00266 CPlane plane1(ANormal * (AVertex2 - AVertex1), AVertex1);
00267 CPlane plane2(ANormal * (AVertex3 - AVertex2), AVertex2);
00268 CPlane plane3(ANormal * (AVertex1 - AVertex3), AVertex3);
00269 return (isPositive(plane1.pointDistance(APoint)) &&
00270 isPositive(plane2.pointDistance(APoint)) &&
00271 isPositive(plane3.pointDistance(APoint)));
00272 }
00273
00274 TCoordinate CGMapVertex::getTriangleCoef(const CVertex & AVertex1,
00275 const CVertex & AVertex2,
00276 const CVertex & AVertex3)
00277 {
00278 if ( CGeometry::areColinear(AVertex1-AVertex2, AVertex2-AVertex3) )
00279 return 0.0;
00280
00281 CVertex AB = AVertex2 - AVertex1;
00282 CVertex AC = AVertex3 - AVertex1;
00283 CVertex BC = AVertex3 - AVertex2;
00284
00285 TCoordinate a = BC.norm();
00286 TCoordinate b = AC.norm();
00287 TCoordinate c = AB.norm();
00288
00289 CPlane plane(AB, (AVertex1 + AVertex2) / 2.0);
00290 CVertex line = (AB * AC) * BC;
00291 CVertex pt;
00292
00293 if (!plane.getLineIntersection((AVertex2 + AVertex3) / 2.0, line, &pt)) {
00294 cerr << "CTriangulationTools::getTriangleCoef: "
00295 << "Erreur dans le calcul du coefficient !" << endl
00296 << AVertex1 << endl
00297 << AVertex2 << endl
00298 << AVertex3 << endl;
00299 assert(false);
00300 }
00301
00302 TCoordinate R = (AVertex1 - pt).norm();
00303 TCoordinate coef = (a*b*c) / (4 * M_PI * R*R*R);
00304 assert(coef >= 0.0 && coef < 1.0);
00305 return coef;
00306 }
00307
00308 #undef VTX
00309
00310