00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 #include "corefine-3d-tools.hh"
00025
00026
00027 #include "message-macros.hh"
00028
00029 using namespace std;
00030 using namespace GMap3d;
00031
00032
00033
00034 #define a0 FMap->alpha0
00035 #define a1 FMap->alpha1
00036 #define a2 FMap->alpha2
00037 #define a3 FMap->alpha3
00038
00039 #define VTX(d) (AVertexDI < 0 ? FMap->findVertex(d) : \
00040 (CAttributeVertex*)FMap->getDirectInfo(d, AVertexDI))
00041
00042
00043
00044 CCorefine3dTools::CCorefine3dTools(CGMapVertex * AMap,
00045 TCoordinate AEpsilon)
00046 : CCorefine2dTools(AMap, AEpsilon)
00047 {
00048 }
00049
00050
00051
00052 bool CCorefine3dTools::isVectorInSector(const CVertex & AVector,
00053 const CVertex & ASectorVector1,
00054 const CVertex & ASectorVector2,
00055 const CPlane & APlane,
00056 bool AThickVector1,
00057 bool AThickVector2)
00058 {
00059 TCoordinate pv = (ASectorVector1 * ASectorVector2).dot(APlane.getNormal());
00060 TCoordinate pv1 = (ASectorVector1 * AVector).dot(APlane.getNormal());
00061 TCoordinate pv2 = (ASectorVector2 * AVector).dot(APlane.getNormal());
00062 bool result = pv < 0.0 ? pv1 >= 0.0 || pv2 <= 0.0 : pv1 >= 0.0 && pv2 <= 0.0;
00063
00064 if (!result && AThickVector1 && pv1 >= -FEps &&
00065 AVector.dot(ASectorVector1) > 0.0) result = true;
00066 if (!result && AThickVector2 && pv2 <= FEps &&
00067 AVector.dot(ASectorVector2) > 0.0) result = true;
00068
00069 return result;
00070 }
00071
00072
00073
00074 bool CCorefine3dTools::isVectorInSector(const CVertex & AVector,
00075 CDart * ASector,
00076 const CPlane & APlane,
00077 bool ATestBorders,
00078 int AVertexDI)
00079 {
00080 CVertex pt = *VTX(ASector);
00081
00082 if (a2(ASector) != a1(ASector))
00083 if (!ATestBorders)
00084 return isVectorInSector(AVector,
00085 *VTX(a0(ASector)) - pt,
00086 *VTX(a0(a1(ASector))) - pt,
00087 APlane);
00088 else
00089 return (isVectorInSector(AVector,
00090 *VTX(a0(ASector)) - pt,
00091 *VTX(a0(a1(ASector))) - pt,
00092 APlane) ||
00093 isPointOnLine(*VTX(a0(ASector)), pt,
00094 pt + AVector) ||
00095 isPointOnLine(*VTX(a0(a1(ASector))), pt,
00096 pt + AVector));
00097 else
00098 return true;
00099 }
00100
00101
00102
00103 bool CCorefine3dTools::areVectorsColinear(const CVertex & AVector1,
00104 const CVertex & AVector2)
00105 {
00106 return CGeneralTools::areVectorsColinear(AVector1, AVector2);
00107 }
00108
00109
00110
00111 TCoordinate CCorefine3dTools::pointParameterOnLine(const CVertex & APoint,
00112 const CVertex & ALineVertex1,
00113 const CVertex & ALineVertex2)
00114 {
00115 return CGeneralTools::pointParameterOnLine(APoint, ALineVertex1,
00116 ALineVertex2);
00117 }
00118
00119
00120
00121 CDart * CCorefine3dTools::findSectorOfVector(const CVertex & AVector,
00122 CDart * AVertex,
00123 const CPlane & APlane,
00124 int AVertexDI)
00125 {
00126 CDart *d = AVertex;
00127 CDart *best = NULL;
00128
00129 do {
00130 if (isVectorInSector(AVector, d, APlane, false, AVertexDI))
00131 best = d;
00132
00133 d = a2(a1(d));
00134 }
00135 while (d != AVertex && best == NULL);
00136
00137 if (best == NULL || !FMap->isSameOrbit(best, AVertex, ORBIT_01)) {
00138 best = AVertex;
00139 WARN_BB("<findSectorOfVector> Aucun secteur valide sur la face !");
00140 }
00141
00142 return best;
00143 }
00144
00145
00146
00147 bool CCorefine3dTools::areFacesCoplanar(const CVertex & AVertex1,
00148 const CPlane & APlane1,
00149 const CVertex & AVertex2,
00150 const CPlane & APlane2)
00151 {
00152 CVertex normal1 = APlane1.getNormal();
00153 CVertex normal2 = APlane2.getNormal();
00154
00155 normal1 = normalizeVector(normal1);
00156 normal2 = normalizeVector(normal2);
00157
00158 return (isVectorNull(normal1 * normal2) &&
00159 isPointOnPlane(AVertex2, APlane1) &&
00160 isPointOnPlane(AVertex1, APlane2));
00161 }
00162
00163
00164
00165 CDart * CCorefine3dTools::findWellOrientedDart(CDart * AVertex, int AVertexDI)
00166 {
00167 TCoordinate n1[4], n2[4];
00168 volumeNormalVector(AVertex, n1, AVertexDI);
00169 volumeNormalVector(a3(AVertex), n2, AVertexDI);
00170
00171 if ((n1[3] < 0.0 && n2[3] < 0.0) ||
00172 (n1[3] > 0.0 && n2[3] > 0.0)) {
00173
00174
00175
00176 CBoundingBox bb1, bb2;
00177
00178 bb1 = orbitBoundingBox(AVertex, ORBIT_VOLUME, AVertexDI);
00179 bb2 = orbitBoundingBox(a3(AVertex), ORBIT_VOLUME, AVertexDI);
00180
00181 if (bb1.getSurface() <= bb2.getSurface()) {
00182
00183 if (n1[3] < 0.0) {
00184 AVertex = a1(AVertex);
00185 }
00186 }
00187 else {
00188
00189 if (n1[3] > 0.0) {
00190 AVertex = a1(AVertex);
00191 }
00192 }
00193 }
00194 else if (n1[3] < 0.0) {
00195 AVertex = a1(AVertex);
00196 }
00197
00198 return AVertex;
00199 }
00200
00201
00202
00203 TPositionOnEdge
00204 CCorefine3dTools::localizeEdgesIntersection(const CVertex & ALinePoint1,
00205 const CVertex & ALinePoint2,
00206 const CVertex & AEdgePoint1,
00207 const CVertex & AEdgePoint2,
00208 TCoordinate * ALineParam,
00209 TCoordinate * AEdgeParam)
00210 {
00211 assert(AEdgeParam != NULL && ALineParam != NULL);
00212
00213 if (localizePointOnEdge(AEdgePoint2, ALinePoint1,
00214 ALinePoint2, ALineParam) != EP_OutOfEdge) {
00215
00216 *AEdgeParam = 1.0;
00217 return EP_OnSecondVertex;
00218 }
00219
00220 if (localizePointOnEdge(AEdgePoint1, ALinePoint1,
00221 ALinePoint2, ALineParam) != EP_OutOfEdge) {
00222 *AEdgeParam = 0.0;
00223 return EP_OnFirstVertex;
00224 }
00225
00226 CPlane plane;
00227
00228 if (sqrt(squareDistanceBetweenLines(ALinePoint1, ALinePoint2,
00229 AEdgePoint1, AEdgePoint2, &plane)) > FEps)
00230 return EP_OutOfEdge;
00231
00232 if (isVectorNull(plane.getNormal()))
00233 return EP_OutOfEdge;
00234
00235 if (!getLinesIntersection(ALinePoint1, ALinePoint2, AEdgePoint1, AEdgePoint2,
00236 ALineParam, AEdgeParam, plane))
00237 return EP_OutOfEdge;
00238
00239 if (*AEdgeParam >= 0.0 && *AEdgeParam <= 1.0)
00240 return EP_OnEdge;
00241
00242 return EP_OutOfEdge;
00243 }
00244
00245
00246
00247 TPositionOnEdge
00248 CCorefine3dTools::localizeEdgeAndPlaneIntersection(CVertex AEdgePoint1,
00249 CVertex AEdgePoint2,
00250 const CPlane & APlane,
00251 TCoordinate * AEdgeParam)
00252 {
00253 assert(AEdgeParam != NULL);
00254
00255 CVertex v = AEdgePoint2 - AEdgePoint1;
00256
00257 if (isPointOnPlane(AEdgePoint2, APlane)) {
00258 *AEdgeParam = 1.0;
00259 return EP_OnSecondVertex;
00260 }
00261
00262 if (isPointOnPlane(AEdgePoint1, APlane)) {
00263 *AEdgeParam = 0.0;
00264 return EP_OnFirstVertex;
00265 }
00266
00267 if (!APlane.getLineIntersection(AEdgePoint1, v, AEdgeParam))
00268 return EP_OutOfEdge;
00269
00270
00271
00272
00273
00274 if (*AEdgeParam >= 0.0 && *AEdgeParam <= 1.0)
00275 return EP_OnEdge;
00276
00277 return EP_OutOfEdge;
00278 }
00279
00280
00281
00282 CEdgeIntersection
00283 CCorefine3dTools::findNearestIntersectionInFace(const CVertex & AOrigin,
00284 const CVertex & ADirection,
00285 CDart * AFace,
00286 const CPlane & APlane,
00287 bool AExcludeOrigin,
00288 int AVertexDI)
00289 {
00290 ENTER;
00291
00292 TProjection proj = APlane.getBestProjection();
00293
00294 CVertex projOrigin = APlane.projectPoint(AOrigin, proj);
00295 CVertex projDirection = APlane.projectPoint(ADirection, proj);
00296
00297 CDart *d = AFace;
00298 CEdgeIntersection inter;
00299 TPositionOnEdge pos;
00300 int edge_mark = FMap->getNewMark();
00301 TCoordinate t1 = 0.0, t2 = 0.0;
00302 CVertex pt1, pt2, pt3;
00303 CVertex extremity = projOrigin + projDirection;
00304
00305 if (projOrigin != extremity) {
00306 MSG("Parcours des sommets de la face");
00307 do {
00308 pt1 = APlane.projectPoint(*VTX(d), proj);
00309
00310 if (CCorefine2dTools::isPointOnLine(pt1, projOrigin, extremity)) {
00311 if (!AExcludeOrigin || !arePointsEqual(projOrigin, pt1)) {
00312 t1 = CCorefine2dTools::pointParameterOnLine(pt1, projOrigin,
00313 extremity);
00314
00315 if (t1 >= 0.0 && (t1 <= inter.getParameter() ||
00316 inter.getCell() == NULL)) {
00317 inter.setCell(d);
00318 inter.setCellDimension(0);
00319 inter.setPosition(EP_OnEdge);
00320 inter.setParameter(t1);
00321 inter.setPoint(*VTX(d));
00322
00323 MSG("une intersection trouvée :" << endl << inter);
00324 }
00325 }
00326
00327
00328 FMap->setMark(d, edge_mark);
00329 FMap->setMark(a0(d), edge_mark);
00330 FMap->setMark(a1(d), edge_mark);
00331 FMap->setMark(a0(a1(d)), edge_mark);
00332 }
00333
00334 d = a1(a0(d));
00335 }
00336 while (d != AFace);
00337
00338 MSG("Parcours des arêtes de la face");
00339 do {
00340 pt1 = APlane.projectPoint(*VTX(d), proj);
00341 pt2 = APlane.projectPoint(*VTX(a0(d)), proj);
00342
00343 if (FMap->isMarked(d, edge_mark)) {
00344 FMap->unsetMark(d, edge_mark);
00345 FMap->unsetMark(a0(d), edge_mark);
00346 }
00347 else if (fabs((pt2 - pt1).norm()) > FEps / 2.0 &&
00348 (!AExcludeOrigin ||
00349 !CCorefine2dTools::isPointOnLine(projOrigin, pt1, pt2))) {
00350 pos = CCorefine2dTools::localizeEdgesIntersection(projOrigin, extremity,
00351 pt1, pt2, &t1, &t2);
00352
00353
00354
00355 if (pos != EP_OutOfEdge &&
00356
00357 t1 >= 0.0 && (t1 <= inter.getParameter() ||
00358 inter.getCell() == NULL)) {
00359 inter.setCell(d);
00360 inter.setCellDimension(1);
00361 inter.setPosition(EP_OnEdge);
00362 inter.setParameter(t1);
00363 inter.setPoint(*VTX(d) + t2 * (*VTX(a0(d)) - *VTX(d)));
00364
00365 MSG("une intersection trouvée :" << endl << inter);
00366 }
00367 }
00368
00369 d = a1(a0(d));
00370 }
00371 while (d != AFace);
00372
00373
00374
00375 }
00376
00377 FMap->freeMark(edge_mark);
00378
00379 EXIT;
00380
00381 return inter;
00382 }
00383
00384
00385
00386 CEdgeIntersection
00387 CCorefine3dTools::findNearestIntersectionInFace(const CVertex & AOrigin,
00388 const CVertex & ADirection,
00389 CDart * AFace,
00390 const CPlane & AFacePlane,
00391 const CPlane & ACutPlane,
00392 bool AExcludeOrigin,
00393 int AVertexDI)
00394 {
00395 DEBUG_FUNCTION;
00396
00397 CDart *d = AFace;
00398 CEdgeIntersection inter;
00399 TCoordinate dist, t;
00400 CVertex pt1, pt2, pt3, dir, extremity = AOrigin + ADirection;
00401
00402 int pos_mark = FMap->getNewMark();
00403 int neg_mark = FMap->getNewMark();
00404 int face_mark = FMap->getNewMark();
00405
00406 FMap->markOrbit(AFace, ORBIT_01, face_mark);
00407
00408 MSG("Parcours des sommets de la face");
00409 do {
00410 pt1 = *VTX(d);
00411 MSG("sommet testé : " << pt1);
00412 dist = ACutPlane.pointDistance(pt1);
00413 if (dist < -FEps) {
00414 MSG("le sommet se trouve du côté négatif");
00415 for (CDynamicCoverageVertex dcv(FMap, d); dcv.cont(); ++dcv)
00416 if (FMap->isMarked(*dcv, face_mark)) FMap->setMark(*dcv, neg_mark);
00417 }
00418 else if (dist > FEps) {
00419 MSG("le sommet se trouve du côté positif");
00420 for (CDynamicCoverageVertex dcv(FMap, d); dcv.cont(); ++dcv)
00421 if (FMap->isMarked(*dcv, face_mark)) FMap->setMark(*dcv, pos_mark);
00422 }
00423 else if (!AExcludeOrigin || !arePointsEqual(pt1, AOrigin)) {
00424 t = pointParameterOnLine(pt1, AOrigin, extremity);
00425 if (t >= 0.0 && (t <= inter.getParameter() || !inter.getCell())) {
00426 MSG("le sommet se trouve sur la ligne de coupe");
00427 inter.setCell(d);
00428 inter.setCellDimension(0);
00429 inter.setPosition(EP_OnEdge);
00430 inter.setParameter(t);
00431 inter.setPoint(pt1);
00432 }
00433 }
00434 d = FMap->alpha01(d);
00435 }
00436 while (d != AFace);
00437
00438 MSG("Parcours des arêtes de la face");
00439 do {
00440 pt1 = *VTX(d);
00441 pt2 = *VTX(a0(d));
00442 MSG("arête testée : [" << pt1 << ";" << pt2 << "]");
00443 dir = pt2 - pt1;
00444 if ((FMap->isMarked(d, pos_mark) && FMap->isMarked(a0(d), neg_mark)) ||
00445 (FMap->isMarked(d, neg_mark) && FMap->isMarked(a0(d), pos_mark))) {
00446 if (!AExcludeOrigin ||
00447 !(abs(CPlane(AFacePlane.getNormal() * (pt2-pt1), pt1).
00448 pointDistance(AOrigin)) <= FEps)) {
00449 if (ACutPlane.getLineIntersection(pt1, dir, &t)) {
00450 assert(t > 0.0 && t < 1.0);
00451 pt3 = pt1 + dir * t;
00452 t = pointParameterOnLine(pt3, AOrigin, extremity);
00453 if (t >= 0.0 && (t <= inter.getParameter() || !inter.getCell())) {
00454 MSG("la ligne de coupe intersecte l'arête en " << pt3);
00455 inter.setCell(d);
00456 inter.setCellDimension(1);
00457 inter.setPosition(EP_OnEdge);
00458 inter.setParameter(t);
00459 inter.setPoint(pt3);
00460 }
00461 }
00462 else
00463 assert(false);
00464 }
00465 }
00466 FMap->unsetMark(d, face_mark);
00467 FMap->unsetMark(a0(d), face_mark);
00468 for (CDynamicCoverageEdge dce(FMap, d); dce.cont(); ++dce) {
00469 FMap->unsetMark(*dce, pos_mark);
00470 FMap->unsetMark(*dce, neg_mark);
00471 }
00472 d = FMap->alpha01(d);
00473 }
00474 while (d != AFace);
00475
00476
00477
00478
00479 FMap->freeMark(pos_mark);
00480 FMap->freeMark(neg_mark);
00481 FMap->freeMark(face_mark);
00482
00483 return inter;
00484 }
00485
00486
00487
00488 CEdgeIntersection
00489 CCorefine3dTools::findNearestIntersectionInOrbit(const CVertex & AOrigin,
00490 const CVertex & ADirection,
00491 CDart * ADart, TOrbit AOrbit,
00492 int AVertexDI)
00493 {
00494 CEdgeIntersection inter;
00495 TCoordinate t1 = 0.0, t2 = 0.0;
00496 CVertex *pt1, *pt2, pt, extremity = AOrigin + ADirection;
00497 CPlane plane;
00498 CDart *face;
00499 int mark = FMap->getNewMark();
00500 int edge_mark = FMap->getNewMark();
00501 int face_mark = FMap->getNewMark();
00502 int orient_mark = FMap->getNewMark();
00503
00504 TOrbit orbit_vertex = AND_ORBIT(ORBIT_VERTEX, AOrbit);
00505 TOrbit orbit_edge = AND_ORBIT(ORBIT_EDGE, AOrbit);
00506 TOrbit orbit_face = AND_ORBIT(ORBIT_FACE, AOrbit);
00507
00508 CCoverage *cov = FMap->getDynamicCoverage(ADart, AOrbit);
00509 CCoverage *dcf, *tmp_cov;
00510
00511 FMap->halfMarkOrbit(ADart, AOrbit, orient_mark);
00512
00513 inter.setCell(NULL);
00514
00515
00516 for (cov->reinit(); cov->cont(); ++(*cov)) {
00517 dcf = FMap->getDynamicCoverage(**cov, orbit_face);
00518
00519 for (; dcf->cont(); ++(*dcf)) {
00520 if (!FMap->isMarked(**dcf, mark) && FMap->isMarked(**dcf, orient_mark)) {
00521
00522 FMap->markOrbit(**dcf, orbit_vertex, mark);
00523
00524 pt = *VTX(**dcf);
00525
00526 if (isPointOnLine(pt, AOrigin, extremity)) {
00527 t2 = pointParameterOnLine(pt, AOrigin, extremity);
00528
00529 if (t2 > 0.0) {
00530 if (t2 < inter.getParameter() || inter.getCell() == NULL) {
00531 inter.setCell(**dcf);
00532 inter.setCellDimension(0);
00533 inter.setPoint(pt);
00534 inter.setParameter(t2);
00535 }
00536
00537 tmp_cov = FMap->getDynamicCoverage(**dcf, orbit_vertex);
00538
00539 for (; tmp_cov->cont(); ++(*tmp_cov)) {
00540
00541 if (!FMap->isMarked(**tmp_cov, edge_mark))
00542 FMap->markOrbit(**tmp_cov, orbit_edge, edge_mark);
00543
00544
00545 if (!FMap->isMarked(**tmp_cov, face_mark))
00546 FMap->markOrbit(**tmp_cov, orbit_face, face_mark);
00547 }
00548
00549 delete tmp_cov;
00550 }
00551 }
00552 }
00553 }
00554
00555 delete dcf;
00556 }
00557
00558
00559 for (cov->reinit(); cov->cont(); ++(*cov)) {
00560 dcf = FMap->getDynamicCoverage(**cov, orbit_face);
00561
00562 for (; dcf->cont(); ++(*dcf)) {
00563 if (FMap->isMarked(**dcf, edge_mark)) {
00564 FMap->unmarkOrbit(**dcf, orbit_edge, edge_mark);
00565 FMap->unmarkOrbit(**dcf, orbit_edge, mark);
00566 }
00567 else if (FMap->isMarked(**dcf, mark) &&
00568 FMap->isMarked(**dcf, orient_mark)) {
00569 FMap->unmarkOrbit(**dcf, orbit_edge, mark);
00570
00571 pt1 = VTX(**dcf);
00572 pt2 = VTX(a0(**dcf));
00573
00574 if (sqrt(squareDistanceBetweenLines(*pt1, *pt2, AOrigin, extremity,
00575 &plane)) <= FEps) {
00576
00577 if (!isVectorNull(plane.getNormal()) &&
00578 getLinesIntersection(*pt1, *pt2, AOrigin, extremity,
00579 &t1, &t2, plane)) {
00580
00581 pt = AOrigin + t2 * ADirection;
00582
00583 if (t1 >= 0.0 && t1 <= 1.0 && t2 > 0) {
00584 if (t2 < inter.getParameter() || inter.getCell() == NULL) {
00585 inter.setCell(**dcf);
00586 inter.setCellDimension(1);
00587 inter.setPoint(pt);
00588 inter.setParameter(t2);
00589 }
00590
00591 tmp_cov = FMap->getDynamicCoverage(**dcf, orbit_edge);
00592
00593 for (; tmp_cov->cont(); ++(*tmp_cov)) {
00594
00595 if (!FMap->isMarked(**tmp_cov, face_mark))
00596 FMap->markOrbit(**tmp_cov, orbit_face, face_mark);
00597 }
00598
00599 delete tmp_cov;
00600 }
00601 }
00602 }
00603 }
00604 }
00605
00606 delete dcf;
00607 }
00608
00609
00610 for (cov->reinit(); cov->cont(); ++(*cov)) {
00611 if (FMap->isMarked(**cov, face_mark)) {
00612 FMap->unmarkOrbit(**cov, orbit_face, face_mark);
00613 FMap->markOrbit(**cov, orbit_face, mark);
00614 }
00615 else if (!FMap->isMarked(**cov, mark)) {
00616 face = (FMap->isMarked(**cov, orient_mark)) ? **cov : a3(**cov);
00617 FMap->markOrbit(face, orbit_face, mark);
00618
00619 plane = facePlane(face, AVertexDI);
00620
00621 if (plane.getLineIntersection(AOrigin, ADirection, &t2)) {
00622 if (t2 > 0.0 &&
00623 (t2 < inter.getParameter() || inter.getCell() == NULL)) {
00624 pt = AOrigin + t2 * ADirection;
00625
00626 if (isPointInFace(pt, face, &plane, AVertexDI)) {
00627 inter.setCell(face);
00628 inter.setCellDimension(2);
00629 inter.setPoint(pt);
00630 inter.setParameter(t2);
00631 }
00632 }
00633 }
00634 }
00635 }
00636
00637 for (cov->reinit(); cov->cont(); ++(*cov)) {
00638 FMap->unsetMark(**cov, mark);
00639 FMap->unsetMark(**cov, orient_mark);
00640 }
00641
00642 delete cov;
00643
00644 FMap->freeMark(mark);
00645 FMap->freeMark(edge_mark);
00646 FMap->freeMark(face_mark);
00647 FMap->freeMark(orient_mark);
00648
00649 return inter;
00650 }
00651
00652
00653
00654 bool CCorefine3dTools::isPointInFace(const CVertex & APoint, CDart * AFace,
00655 const CPlane * APlane, int AVertexDI)
00656 {
00657 bool result;
00658 CPlane plane;
00659
00660 if (APlane == NULL)
00661 plane = facePlane(AFace, AVertexDI);
00662 else
00663 plane = *APlane;
00664
00665 if (APoint == *VTX(AFace))
00666 result = true;
00667 else {
00668 CVertex direction = *VTX(AFace) - APoint;
00669 CEdgeIntersection inter;
00670
00671 inter = findNearestIntersectionInFace(APoint, direction, AFace, plane,
00672 false, AVertexDI);
00673
00674 if (inter.getCell() == NULL) {
00675 result = false;
00676 }
00677 else {
00678 if (inter.getCellDimension() == 0) {
00679 result = isVectorInSector(-direction, inter.getCell(), plane, false,
00680 AVertexDI);
00681 }
00682 else {
00683 CVertex v = edgeVector(inter.getCell(), AVertexDI);
00684 result = (plane.getNormal() * v).dot(direction) < 0.0;
00685 }
00686 }
00687 }
00688
00689 return result;
00690 }
00691
00692
00693 #define LINK(d1, d2) (FMap->setDirectInfo(d1, ALinkDirectInfo, d2), \
00694 FMap->setDirectInfo(d2, ALinkDirectInfo, d1), \
00695 FMap->setDirectInfo(a0(d1), ALinkDirectInfo, a0(d2)), \
00696 FMap->setDirectInfo(a0(d2), ALinkDirectInfo, a0(d1)))
00697
00698 #define UNLINK(d) (FMap->setDirectInfo(d, ALinkDirectInfo, NULL), \
00699 FMap->setDirectInfo(a0(d), ALinkDirectInfo, NULL))
00700
00701 #define ALPHA2(d) ((FMap->getDirectInfo(d, ALinkDirectInfo) != NULL) ? \
00702 (CDart*)FMap->getDirectInfo(d, ALinkDirectInfo) : \
00703 alpha2(d))
00704
00705 void CCorefine3dTools::sortFacesAroundEdges_Naive(CDart * AEdge1,
00706 CDart * AEdge2,
00707 int ALinkDirectInfo,
00708 int AVertexDI)
00709 {
00710 CPlane plane(edgeVector(AEdge1, AVertexDI), *VTX(AEdge1));
00711
00712 CDart *d1, *d2;
00713 CVertex v, v1, v2;
00714 CDart *best, *tmp1, *tmp2;
00715
00716 d1 = AEdge1;
00717 do {
00718 v = faceNormalVector(d1, AVertexDI);
00719
00720 best = NULL;
00721 d2 = AEdge2;
00722 v2 = faceNormalVector(d2, AVertexDI);
00723 do {
00724 tmp2 = a3(alpha2(d2));
00725
00726 v1 = v2;
00727 v2 = faceNormalVector(tmp2, AVertexDI);
00728
00729 if (isVectorInSector(v, v1, v2, plane))
00730 best = d2;
00731
00732 d2 = tmp2;
00733 }
00734 while (d2 != AEdge2 && best == NULL);
00735
00736 if (best == NULL) {
00737 WARN_BB("<sortFacesAroundEdges_Naive> Aucun secteur trouvé !");
00738 d2 = AEdge2;
00739 }
00740 else
00741 d2 = best;
00742
00743
00744 tmp1 = alpha2(d1);
00745 if (FMap->getDirectInfo(tmp1, ALinkDirectInfo) == d2) {
00746 UNLINK(tmp1);
00747 UNLINK(d2);
00748 }
00749 if (FMap->getDirectInfo(d2, ALinkDirectInfo) == NULL)
00750 LINK(d2, a3(d1));
00751
00752
00753 tmp1 = alpha2(a3(d1));
00754 tmp2 = alpha2(d2);
00755 if (FMap->getDirectInfo(tmp1, ALinkDirectInfo) == tmp2) {
00756 UNLINK(tmp1);
00757 UNLINK(tmp2);
00758 }
00759 if (FMap->getDirectInfo(tmp2, ALinkDirectInfo) == NULL)
00760 LINK(d1, tmp2);
00761
00762 d1 = a3(alpha2(d1));
00763 }
00764 while (d1 != AEdge1);
00765 }
00766
00767
00768
00769 void CCorefine3dTools::sortFacesAroundEdges_NaiveBis(CDart * AEdge1,
00770 CDart * AEdge2,
00771 int ALinkDirectInfo,
00772 int AVertexDI)
00773 {
00774 CPlane plane(edgeVector(AEdge1, AVertexDI), *VTX(AEdge1));
00775
00776 CDart *d1, *d2;
00777 CVertex v, v1, v2;
00778 CDart *best, *tmp;
00779 int mark = FMap->getNewMark();
00780
00781 for (CDynamicCoverageEdge dce1(FMap, AEdge1); dce1.cont(); dce1++)
00782 FMap->setDirectInfo(*dce1, ALinkDirectInfo, NULL);
00783 for (CDynamicCoverageEdge dce2(FMap, AEdge2); dce2.cont(); dce2++)
00784 FMap->setDirectInfo(*dce2, ALinkDirectInfo, NULL);
00785
00786 d1 = AEdge1;
00787 d2 = AEdge2;
00788 do {
00789 v = faceNormalVector(d1, AVertexDI);
00790
00791 best = NULL;
00792 v2 = faceNormalVector(d2, AVertexDI);
00793 do {
00794 FMap->setMark(d2, mark);
00795
00796
00797 tmp = a3(ALPHA2(d2));
00798
00799 v1 = v2;
00800 v2 = faceNormalVector(tmp, AVertexDI);
00801
00802 if (isVectorInSector(v, v1, v2, plane))
00803 best = d2;
00804
00805 d2 = tmp;
00806 }
00807 while (best == NULL && !FMap->isMarked(d2, mark));
00808
00809 FMap->unmarkOrbit(AEdge1, ORBIT_23, mark);
00810 FMap->unmarkOrbit(AEdge2, ORBIT_23, mark);
00811
00812 if (best == NULL) {
00813 WARN_BB("<sortFacesAroundEdges_NaiveBis> Aucun secteur trouvé !");
00814 d2 = AEdge2;
00815 }
00816 else
00817 d2 = best;
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832 if (FMap->getDirectInfo(d2, ALinkDirectInfo) != NULL) {
00833 tmp = (CDart*)FMap->getDirectInfo(d2, ALinkDirectInfo);
00834 UNLINK(tmp);
00835 UNLINK(d2);
00836 }
00837 else
00838 tmp = alpha2(d2);
00839
00840 LINK(d2, a3(d1));
00841 LINK(d1, tmp);
00842
00843 d1 = a3(alpha2(d1));
00844 }
00845 while (d1 != AEdge1);
00846
00847 FMap->freeMark(mark);
00848 }
00849
00850
00851
00852 bool CCorefine3dTools::sortFacesAroundEdges_Optimized(CDart * AEdge1,
00853 CDart * AEdge2,
00854 int ALinkDirectInfo,
00855 int AVertexDI)
00856 {
00857 int mark = FMap->getNewMark();
00858
00859 CPlane plane(edgeVector(AEdge1, AVertexDI), *VTX(AEdge1));
00860
00861 CDart *d = AEdge1;
00862 CVertex v = faceNormalVector(d, AVertexDI);
00863
00864 CDart *tmp = a3(alpha2(AEdge2));
00865 CVertex v1 = faceNormalVector(AEdge2, AVertexDI);
00866 CVertex v2 = faceNormalVector(tmp, AVertexDI);
00867
00868
00869
00870
00871 while (!isVectorInSector(v, v1, v2, plane)) {
00872
00873 v1 = v2;
00874 AEdge2 = tmp;
00875
00876 tmp = a3(alpha2(AEdge2));
00877 v2 = faceNormalVector(tmp, AVertexDI);
00878 }
00879
00880 CDart *d1 = a3(alpha2(AEdge1));
00881 CDart *d2 = a3(alpha2(AEdge2));
00882
00883 v1 = faceNormalVector(d1, AVertexDI);
00884 v2 = faceNormalVector(d2, AVertexDI);
00885
00886 do {
00887 FMap->setMark(d, mark);
00888
00889 if (areVectorsColinear(v1, v2)) {
00890 if (!FMap->isMarked(d2, mark)) {
00891 LINK(d, a3(d2));
00892 FMap->setMark(d2, mark);
00893 d = d2;
00894 }
00895
00896 LINK(d, a3(d1));
00897
00898 d = d1;
00899 v = v1;
00900
00901 d1 = a3(alpha2(d1));
00902 d2 = a3(alpha2(d2));
00903
00904 v1 = faceNormalVector(d1, AVertexDI);
00905 v2 = faceNormalVector(d2, AVertexDI);
00906 }
00907 else if (isVectorInSector(v2, v, v1, plane)) {
00908 LINK(d, a3(d2));
00909
00910 d = d2;
00911 v = v2;
00912
00913 d2 = a3(alpha2(d2));
00914 v2 = faceNormalVector(d2, AVertexDI);
00915 }
00916 else {
00917 LINK(d, a3(d1));
00918
00919 d = d1;
00920 v = v1;
00921
00922 d1 = a3(alpha2(d1));
00923 v1 = faceNormalVector(d1, AVertexDI);
00924 }
00925 }
00926 while (!FMap->isMarked(d, mark));
00927
00928 bool problem = false;
00929
00930
00931 d2 = AEdge2;
00932 do {
00933 if (FMap->isMarked(d2, mark))
00934 FMap->unsetMark(d2, mark);
00935 else {
00936 WARN_BR("<sortFacesAroundEdges> Brin de la deuxième arête non marqué !");
00937 problem = true;
00938 }
00939
00940 d2 = a3(alpha2(d2));
00941 }
00942 while (d2 != AEdge2);
00943
00944
00945 d1 = AEdge1;
00946 do {
00947 if (FMap->isMarked(d1, mark))
00948 FMap->unsetMark(d1, mark);
00949 else {
00950 WARN_BR("<sortFacesAroundEdges> Brin de la première arête non marqué !");
00951 problem = true;
00952 }
00953
00954 d1 = a3(alpha2(d1));
00955 }
00956 while (d1 != AEdge1);
00957
00958 FMap->freeMark(mark);
00959
00960 if (problem) {
00961 d1 = AEdge1;
00962 do {
00963 FMap->setDirectInfo(d1, ALinkDirectInfo, NULL);
00964 d1 = a3(alpha2(d1));
00965 }
00966 while (d1 != AEdge1);
00967
00968 d2 = AEdge2;
00969 do {
00970 FMap->setDirectInfo(d2, ALinkDirectInfo, NULL);
00971 d2 = a3(alpha2(d2));
00972 }
00973 while (d2 != AEdge2);
00974 }
00975
00976 return !problem;
00977 }
00978
00979
00980
00981 class CFaceComparator
00982 {
00983 public:
00984 CFaceComparator(CCorefine3dTools & ATools, const CVertex & AAxis,
00985 int AVertexDI = -1)
00986 : FTools(ATools), FAxis(AAxis), FVertexDI(AVertexDI) {}
00987
00988 bool operator () (CDart * AFace1, CDart * AFace2)
00989 {
00990 CVertex v1 = FTools.faceNormalVector(AFace1, FVertexDI);
00991 CVertex v2 = FTools.faceNormalVector(AFace2, FVertexDI);
00992 return (v1 * v2).dot(FAxis) > 0.0;
00993 }
00994
00995 private:
00996 CCorefine3dTools & FTools;
00997 CVertex FAxis;
00998 int FVertexDI;
00999 };
01000
01001
01002
01003 void CCorefine3dTools::sortFacesAroundEdges_SuperNaive(CDart * AEdge1,
01004 CDart * AEdge2,
01005 int ALinkDirectInfo,
01006 int AVertexDI)
01007 {
01008 CDart *d1, *d2;
01009 CVertex axis = edgeVector(AEdge1, AVertexDI);
01010 CFaceComparator comp(*this, axis, AVertexDI);
01011 multiset<CDart*, CFaceComparator> faces(comp);
01012 multiset<CDart*, CFaceComparator>::iterator it, tmp_it;
01013
01014 for (CDynamicCoverageEdge dce1(FMap, AEdge1); dce1.cont(); ++dce1)
01015 FMap->setDirectInfo(*dce1, ALinkDirectInfo, NULL);
01016 for (CDynamicCoverageEdge dce2(FMap, AEdge2); dce2.cont(); ++dce2)
01017 FMap->setDirectInfo(*dce2, ALinkDirectInfo, NULL);
01018
01019 d1 = AEdge1;
01020 do {
01021 faces.insert(d1);
01022 d1 = a3(alpha2(d1));
01023 }
01024 while (d1 != AEdge1);
01025
01026 d2 = AEdge2;
01027 do {
01028 faces.insert(d2);
01029 d2 = a3(alpha2(d2));
01030 }
01031 while (d2 != AEdge2);
01032
01033 for (it = faces.begin(); it != faces.end(); ) {
01034 tmp_it = it++;
01035
01036 d1 = *tmp_it;
01037 d2 = (it != faces.end()) ? *it : *faces.begin();
01038
01039 LINK(d1, a3(d2));
01040 }
01041 }
01042
01043 #undef LINK
01044 #undef UNLINK
01045 #undef ALPHA2
01046
01047
01048 TCoordinate CCorefine3dTools::matrixDet(TCoordinate AMat[3][3])
01049 {
01050 return (AMat[0][0] * (AMat[1][1] * AMat[2][2] - AMat[2][1] * AMat[1][2]) -
01051 AMat[0][1] * (AMat[1][0] * AMat[2][2] - AMat[2][0] * AMat[1][2]) -
01052 AMat[0][2] * (AMat[1][0] * AMat[2][1] - AMat[2][0] * AMat[1][1]));
01053 }
01054
01055
01056
01057 void CCorefine3dTools::vectProduct4d(const CVertex & AV1,
01058 const CVertex & AV2,
01059 const CVertex & AV3,
01060 TCoordinate ANormal[4])
01061 {
01062 TCoordinate mat[3][3];
01063
01064
01065
01066
01067
01068
01069
01070
01071
01072
01073
01074
01075
01076
01077
01078
01079
01080
01081
01082
01083
01084
01085
01086
01087
01088
01089
01090
01091
01092
01093
01094
01095
01096
01097
01098
01099
01100
01101
01102
01103
01104
01105
01106 mat[0][0] = AV1.getX();
01107 mat[0][1] = AV1.getY();
01108 mat[0][2] = AV1.getZ();
01109
01110 mat[1][0] = AV2.getX();
01111 mat[1][1] = AV2.getY();
01112 mat[1][2] = AV2.getZ();
01113
01114 mat[2][0] = AV3.getX();
01115 mat[2][1] = AV3.getY();
01116 mat[2][2] = AV3.getZ();
01117
01118 ANormal[3] = matrixDet(mat);
01119 }
01120
01121
01122
01123 void CCorefine3dTools::volumeNormalVector(CDart * AVolume,
01124 TCoordinate ANormal[4],
01125 int AVertexDI)
01126 {
01127 int orient_mark = FMap->getNewMark();
01128 CVertex v1, v2, v3;
01129 TCoordinate tmp[4];
01130
01131 ANormal[0] = ANormal[1] = ANormal[2] = ANormal[3] = 0.0;
01132
01133 FMap->halfMarkOrbit(AVolume, ORBIT_VOLUME, orient_mark);
01134
01135 for (CDynamicCoverageVolume dcv(FMap, AVolume); dcv.cont(); dcv++) {
01136 if (FMap->isMarked(*dcv, orient_mark)) {
01137 FMap->unsetMark(*dcv, orient_mark);
01138
01139 if (!FMap->isFree2(*dcv)) {
01140 v1 = edgeVector(a1(a2(*dcv)), AVertexDI);
01141 v2 = edgeVector(*dcv, AVertexDI);
01142 v3 = edgeVector(a1(*dcv), AVertexDI);
01143
01144 vectProduct4d(v1, v2, v3, tmp);
01145
01146
01147
01148 ANormal[3] += tmp[3];
01149 }
01150 }
01151 }
01152
01153 FMap->freeMark(orient_mark);
01154 }
01155
01156
01157