Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed
  • src/Tesselation/boundary.cpp

    r5f7b95 r2b7375  
    513513 * \return volume difference between the non- and the created convex envelope
    514514 */
    515 double ConvexizeNonconvexEnvelope(
    516     Tesselation *&TesselStruct,
    517     const molecule * const mol,
    518     const char * const filename,
    519     bool DebugOutputEveryStep)
     515double ConvexizeNonconvexEnvelope(class Tesselation *&TesselStruct, const molecule * const mol, const char * const filename)
    520516{
    521517        //Info FunctionInfo(__func__);
     
    539535  }
    540536
    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
    547538  do {
    548539    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
    558544    PointRunner = TesselStruct->PointsOnBoundary.begin();
    559545    PointAdvance = PointRunner; // we need an advanced point, as the PointRunner might get removed
     
    561547      PointAdvance++;
    562548      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;
    570560          break;
     561        }
    571562      }
    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);
    577583          Concavity = true;
    578584        }
    579585      }
    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 lines
    591     // when the tetraeder does not intersect with other already present triangles
    592     LineRunner = TesselStruct->LinesOnBoundary.begin();
    593     LineAdvance = LineRunner;  // we need an advanced line, as the LineRunner might get removed
    594     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 points
    601         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 baselines
    606               if (!line->ContainsBoundaryPoint(runner->second->endpoints[j])) // and neither of its endpoints
    607                 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 then
    615         // would create a degenerate triangle
    616         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 line
    637           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 endpoint
    654             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 connected
    660                 // only then are we in danger of cutting through it and need to check
    661                 // sign of normal vector and intersecting line
    662                 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 either
    686             // endpoint, i.e. check that scalar product between vectors going
    687             // from intersction to either endpoint has negative sign (both
    688             // 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 volume
    700             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't
    713             LOG(2, "DEBUG: We don't added concave line " << *line
    714                 << " as other line intersects with adjacent triangles.");
    715           }
    716         }
    717       }
    718586      LineRunner = LineAdvance;
    719587    }
    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//  }
    735612
    736613  CalculateConcavityPerBoundaryPoint(TesselStruct);
     
    738615
    739616  // end
    740   LOG(0, "RESULT: Added volume in convexization is " << volume << ".");
     617  LOG(0, "Volume is " << volume << ".");
    741618  return volume;
    742619};
     
    889766//  status = CheckListOfBaselines(TesselStruct);
    890767
    891 //  cout << "before correction" << endl;
     768  cout << "before correction" << endl;
    892769
    893770  // store before correction
     
    902779  // write final envelope
    903780  CalculateConcavityPerBoundaryPoint(TesselStruct);
    904 //  cout << "after correction" << endl;
     781  cout << "after correction" << endl;
    905782  StoreTrianglesinFile(mol, TesselStruct, filename, "");
    906783
Note: See TracChangeset for help on using the changeset viewer.