Changes in src/Tesselation/boundary.cpp [5f7b95:2b7375]
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
src/Tesselation/boundary.cpp
r5f7b95 r2b7375 513 513 * \return volume difference between the non- and the created convex envelope 514 514 */ 515 double ConvexizeNonconvexEnvelope( 516 Tesselation *&TesselStruct, 517 const molecule * const mol, 518 const char * const filename, 519 bool DebugOutputEveryStep) 515 double ConvexizeNonconvexEnvelope(class Tesselation *&TesselStruct, const molecule * const mol, const char * const filename) 520 516 { 521 517 //Info FunctionInfo(__func__); … … 539 535 } 540 536 541 LOG(1, "INFO: Making tesselated surface with " << TesselStruct->TrianglesOnBoundaryCount 542 << " convex ..."); 543 544 // first purge all degenerate triangles 545 TesselStruct->RemoveDegeneratedTriangles(); 546 537 // First step: RemovePointFromTesselatedSurface 547 538 do { 548 539 Concavity = false; 549 550 if (DebugOutputEveryStep) { 551 sprintf(dummy, "-%d", run++); 552 //CalculateConcavityPerBoundaryPoint(TesselStruct); 553 LOG(1, "INFO: Writing " << run << "th tesselation file."); 554 StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, dummy); 555 } 556 557 // first step: remove all full-concave point 540 sprintf(dummy, "-first-%d", run); 541 //CalculateConcavityPerBoundaryPoint(TesselStruct); 542 StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, dummy); 543 558 544 PointRunner = TesselStruct->PointsOnBoundary.begin(); 559 545 PointAdvance = PointRunner; // we need an advanced point, as the PointRunner might get removed … … 561 547 PointAdvance++; 562 548 point = PointRunner->second; 563 LOG(2, "INFO: Current point is " << *point << "."); 564 // check that at least a single line is concave 565 LineMap::iterator LineRunner = point->lines.begin(); 566 for (; LineRunner != point->lines.end(); LineRunner++) { 567 const BoundaryLineSet * line = LineRunner->second; 568 LOG(3, "INFO: Current line of point " << *point << " is " << *line << "."); 569 if (!line->CheckConvexityCriterion()) 549 LOG(1, "INFO: Current point is " << *point << "."); 550 for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) { 551 line = LineRunner->second; 552 LOG(1, "INFO: Current line of point " << *point << " is " << *line << "."); 553 if (!line->CheckConvexityCriterion()) { 554 // remove the point if needed 555 LOG(1, "... point " << *point << " cannot be on convex envelope."); 556 volume += TesselStruct->RemovePointFromTesselatedSurface(point); 557 sprintf(dummy, "-first-%d", ++run); 558 StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, dummy); 559 Concavity = true; 570 560 break; 561 } 571 562 } 572 // remove the point if needed 573 if (LineRunner != point->lines.end()) { 574 const double tmp = TesselStruct->RemoveFullConcavePointFromTesselatedSurface(point); 575 if (tmp > 0.) { 576 volume += tmp; 563 PointRunner = PointAdvance; 564 } 565 566 sprintf(dummy, "-second-%d", run); 567 //CalculateConcavityPerBoundaryPoint(TesselStruct); 568 StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, dummy); 569 570 // second step: PickFarthestofTwoBaselines 571 LineRunner = TesselStruct->LinesOnBoundary.begin(); 572 LineAdvance = LineRunner; // we need an advanced line, as the LineRunner might get removed 573 while (LineRunner != TesselStruct->LinesOnBoundary.end()) { 574 LineAdvance++; 575 line = LineRunner->second; 576 LOG(1, "INFO: Picking farthest baseline for line is " << *line << "."); 577 // take highest of both lines 578 if (TesselStruct->IsConvexRectangle(line) == NULL) { 579 const double tmp = TesselStruct->PickFarthestofTwoBaselines(line); 580 volume += tmp; 581 if (tmp != 0.) { 582 TesselStruct->FlipBaseline(line); 577 583 Concavity = true; 578 584 } 579 585 } 580 PointRunner = PointAdvance;581 }582 583 if (DebugOutputEveryStep) {584 sprintf(dummy, "-%d", run++);585 //CalculateConcavityPerBoundaryPoint(TesselStruct);586 LOG(1, "INFO: Writing " << run << "th tesselation file.");587 StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, dummy);588 }589 590 // second step: flip baselines, i.e. add general tetraeder at concave lines591 // when the tetraeder does not intersect with other already present triangles592 LineRunner = TesselStruct->LinesOnBoundary.begin();593 LineAdvance = LineRunner; // we need an advanced line, as the LineRunner might get removed594 std::map<double, std::pair<BoundaryLineSet *, double> > GainMap;595 while (LineRunner != TesselStruct->LinesOnBoundary.end()) {596 LineAdvance++;597 line = LineRunner->second;598 if (!line->CheckConvexityCriterion()) {599 LOG(2, "INFO: concave line is " << *line << ".");600 // gather the other points601 BoundaryPointSet *BPS[4];602 int m = 0;603 {604 for (TriangleMap::iterator runner = line->triangles.begin(); runner != line->triangles.end(); runner++)605 for (int j = 0; j < 3; j++) // all of their endpoints and baselines606 if (!line->ContainsBoundaryPoint(runner->second->endpoints[j])) // and neither of its endpoints607 BPS[m++] = runner->second->endpoints[j];608 }609 BPS[2] = line->endpoints[0];610 BPS[3] = line->endpoints[1];611 LOG(3, "DEBUG: other line would consist of " << *BPS[0] << " and "612 << *BPS[1] << ".");613 614 // check for already present (third) side of the tetraeder as we then615 // would create a degenerate triangle616 bool TetraederSidePresent = false;617 {618 class TesselPoint *TriangleCandidates[3];619 TriangleCandidates[0] = BPS[0]->node;620 TriangleCandidates[1] = BPS[1]->node;621 TriangleCandidates[2] = BPS[2]->node;622 if ((TesselStruct->GetPresentTriangle(TriangleCandidates) != NULL)) {623 LOG(2, "REJECT: Triangle side " << *TriangleCandidates[0] << ","624 << *TriangleCandidates[1] << "," << *TriangleCandidates[2] << " present.");625 TetraederSidePresent = true;626 }627 TriangleCandidates[2] = BPS[3]->node;628 if ((TesselStruct->GetPresentTriangle(TriangleCandidates) != NULL)) {629 LOG(2, "REJECT: Triangle side " << *TriangleCandidates[0] << ","630 << *TriangleCandidates[1] << "," << *TriangleCandidates[2] << " present.");631 TetraederSidePresent = true;632 }633 }634 635 if ((BPS[0] != BPS[1]) && (m == 2) && (!TetraederSidePresent)) {636 // check whether all adjacent triangles do not intersect with new line637 bool no_line_intersects = true;638 Vector Intersection;639 TriangleSet triangles;640 TriangleSet *firsttriangles = TesselStruct->GetAllTriangles(line->endpoints[0]);641 TriangleSet *secondtriangles = TesselStruct->GetAllTriangles(line->endpoints[1]);642 triangles.insert( firsttriangles->begin(), firsttriangles->end() );643 triangles.insert( secondtriangles->begin(), secondtriangles->end() );644 delete firsttriangles;645 delete secondtriangles;646 for (TriangleSet::const_iterator triangleiter = triangles.begin();647 triangleiter != triangles.end(); ++triangleiter) {648 const BoundaryTriangleSet * triangle = *triangleiter;649 bool line_intersects = triangle->GetIntersectionInsideTriangle(650 BPS[0]->node->getPosition(),651 BPS[1]->node->getPosition(),652 Intersection);653 // switch result when coinciding with endpoint654 bool concave_adjacent_line = false;655 bool intersection_is_endnode = false;656 for (int j=0;j<2;++j) {657 if (Intersection.DistanceSquared(BPS[j]->node->getPosition()) < MYEPSILON) {658 intersection_is_endnode = true;659 // check whether its an adjacent triangle and if it's concavely connected660 // only then are we in danger of cutting through it and need to check661 // sign of normal vector and intersecting line662 for (int i=0;i<2;++i)663 for (int lineindex=0;lineindex < 3;++lineindex)664 if ((triangle->lines[lineindex]->ContainsBoundaryPoint(line->endpoints[i]))665 && (triangle->lines[lineindex]->ContainsBoundaryPoint(BPS[j]))) {666 concave_adjacent_line = !triangle->lines[lineindex]->CheckConvexityCriterion();667 }668 if (concave_adjacent_line) {669 const Vector intersector =670 BPS[(j+1)%2]->node->getPosition() - Intersection;671 if (triangle->NormalVector.ScalarProduct(intersector) >= -MYEPSILON) {672 LOG(4, "ACCEPT: Intersection coincides with first endpoint "673 << *BPS[j] << ".");674 line_intersects = false;675 } else {676 LOG(4, "REJECT: Intersection ends on wrong side of triangle.");677 }678 } else {679 LOG(4, "ACCEPT: Intersection coincides with first endpoint "680 << *BPS[j] << ".");681 line_intersects = false;682 }683 }684 }685 // if we have an intersection, check that it is within either686 // endpoint, i.e. check that scalar product between vectors going687 // from intersction to either endpoint has negative sign (both688 // vectors point in opposite directions)689 if (!intersection_is_endnode && line_intersects) {690 const Vector firstvector = BPS[0]->node->getPosition() - Intersection;691 const Vector secondvector = BPS[1]->node->getPosition() - Intersection;692 if (firstvector.ScalarProduct(secondvector) >= 0)693 line_intersects = false;694 }695 no_line_intersects &= !line_intersects;696 }697 698 if (no_line_intersects) {699 // calculate the volume700 const double tmp = line->CalculateConvexity();701 const double gain =702 CalculateVolumeofGeneralTetraeder(703 BPS[0]->node->getPosition(),704 BPS[1]->node->getPosition(),705 BPS[2]->node->getPosition(),706 BPS[3]->node->getPosition());707 708 GainMap.insert(std::make_pair(tmp, std::make_pair(line,gain) ));709 LOG(2, "DEBUG: Adding concave line " << *line << " with gain of "710 << gain << ".");711 } else {712 // if 2 or 3 don't713 LOG(2, "DEBUG: We don't added concave line " << *line714 << " as other line intersects with adjacent triangles.");715 }716 }717 }718 586 LineRunner = LineAdvance; 719 587 } 720 // flip line with most gain 721 if (!GainMap.empty()) { 722 line = GainMap.begin()->second.first; 723 const double tmp = GainMap.begin()->second.second; 724 volume += tmp; 725 726 // GainMap.clear(); 727 728 // and flip the line 729 LOG(1, "INFO: Flipping current most concave line " << *line << " with gain of " 730 << tmp << "."); 731 TesselStruct->FlipBaseline(line); 732 Concavity = true; 733 } 734 } while ((Concavity)); // && (run < 100) 588 run++; 589 } while (Concavity); 590 //CalculateConcavityPerBoundaryPoint(TesselStruct); 591 //StoreTrianglesinFile(mol, filename, "-third"); 592 593 // third step: IsConvexRectangle 594 // LineRunner = TesselStruct->LinesOnBoundary.begin(); 595 // LineAdvance = LineRunner; // we need an advanced line, as the LineRunner might get removed 596 // while (LineRunner != TesselStruct->LinesOnBoundary.end()) { 597 // LineAdvance++; 598 // line = LineRunner->second; 599 // LOG(1, "INFO: Current line is " << *line << "."); 600 // //if (LineAdvance != TesselStruct->LinesOnBoundary.end()) 601 // //LOG(1, "INFO: Next line will be " << *(LineAdvance->second) << "."); 602 // if (!line->CheckConvexityCriterion(out)) { 603 // LOG(1, "INFO: ... line " << *line << " is concave, flipping it."); 604 // 605 // // take highest of both lines 606 // point = TesselStruct->IsConvexRectangle(line); 607 // if (point != NULL) 608 // volume += TesselStruct->RemovePointFromTesselatedSurface(point); 609 // } 610 // LineRunner = LineAdvance; 611 // } 735 612 736 613 CalculateConcavityPerBoundaryPoint(TesselStruct); … … 738 615 739 616 // end 740 LOG(0, " RESULT: Added volume in convexizationis " << volume << ".");617 LOG(0, "Volume is " << volume << "."); 741 618 return volume; 742 619 }; … … 889 766 // status = CheckListOfBaselines(TesselStruct); 890 767 891 //cout << "before correction" << endl;768 cout << "before correction" << endl; 892 769 893 770 // store before correction … … 902 779 // write final envelope 903 780 CalculateConcavityPerBoundaryPoint(TesselStruct); 904 //cout << "after correction" << endl;781 cout << "after correction" << endl; 905 782 StoreTrianglesinFile(mol, TesselStruct, filename, ""); 906 783
Note:
See TracChangeset
for help on using the changeset viewer.