Changes in src/Tesselation/boundary.cpp [5f7b95:2b7375]

File:

• src/Tesselation/boundary.cpp (modified) (6 diffs)

src/Tesselation/boundary.cpp

@@ -513 +513 @@
  * \return volume difference between the non- and the created convex envelope
  */
-double ConvexizeNonconvexEnvelope(
-    Tesselation *&TesselStruct,
-    const molecule * const mol,
-    const char * const filename,
-    bool DebugOutputEveryStep)
+double ConvexizeNonconvexEnvelope(class Tesselation *&TesselStruct, const molecule * const mol, const char * const filename)
 {
         //Info FunctionInfo(__func__);
@@ -539 +535 @@
   }
 
-  LOG(1, "INFO: Making tesselated surface with " << TesselStruct->TrianglesOnBoundaryCount
-      << " convex ...");
-
-  // first purge all degenerate triangles
-  TesselStruct->RemoveDegeneratedTriangles();
-
+  // First step: RemovePointFromTesselatedSurface
   do {
     Concavity = false;
-
-    if (DebugOutputEveryStep) {
-      sprintf(dummy, "-%d", run++);
-      //CalculateConcavityPerBoundaryPoint(TesselStruct);
-      LOG(1, "INFO: Writing " << run << "th tesselation file.");
-      StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, dummy);
-    }
-
-    // first step: remove all full-concave point
+    sprintf(dummy, "-first-%d", run);
+    //CalculateConcavityPerBoundaryPoint(TesselStruct);
+    StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, dummy);
+
     PointRunner = TesselStruct->PointsOnBoundary.begin();
     PointAdvance = PointRunner; // we need an advanced point, as the PointRunner might get removed
@@ -561 +547 @@
       PointAdvance++;
       point = PointRunner->second;
-      LOG(2, "INFO: Current point is " << *point << ".");
-      // check that at least a single line is concave
-      LineMap::iterator LineRunner = point->lines.begin();
-      for (; LineRunner != point->lines.end(); LineRunner++) {
-        const BoundaryLineSet * line = LineRunner->second;
-        LOG(3, "INFO: Current line of point " << *point << " is " << *line << ".");
-        if (!line->CheckConvexityCriterion())
+      LOG(1, "INFO: Current point is " << *point << ".");
+      for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
+        line = LineRunner->second;
+        LOG(1, "INFO: Current line of point " << *point << " is " << *line << ".");
+        if (!line->CheckConvexityCriterion()) {
+          // remove the point if needed
+          LOG(1, "... point " << *point << " cannot be on convex envelope.");
+          volume += TesselStruct->RemovePointFromTesselatedSurface(point);
+          sprintf(dummy, "-first-%d", ++run);
+          StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, dummy);
+          Concavity = true;
           break;
+        }
       }
-      // remove the point if needed
-      if (LineRunner != point->lines.end()) {
-        const double tmp = TesselStruct->RemoveFullConcavePointFromTesselatedSurface(point);
-        if (tmp > 0.) {
-          volume += tmp;
+      PointRunner = PointAdvance;
+    }
+
+    sprintf(dummy, "-second-%d", run);
+    //CalculateConcavityPerBoundaryPoint(TesselStruct);
+    StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, dummy);
+
+    // second step: PickFarthestofTwoBaselines
+    LineRunner = TesselStruct->LinesOnBoundary.begin();
+    LineAdvance = LineRunner;  // we need an advanced line, as the LineRunner might get removed
+    while (LineRunner != TesselStruct->LinesOnBoundary.end()) {
+      LineAdvance++;
+      line = LineRunner->second;
+      LOG(1, "INFO: Picking farthest baseline for line is " << *line << ".");
+      // take highest of both lines
+      if (TesselStruct->IsConvexRectangle(line) == NULL) {
+        const double tmp = TesselStruct->PickFarthestofTwoBaselines(line);
+        volume += tmp;
+        if (tmp != 0.) {
+          TesselStruct->FlipBaseline(line);
           Concavity = true;
         }
       }
-      PointRunner = PointAdvance;
-    }
-
-    if (DebugOutputEveryStep) {
-      sprintf(dummy, "-%d", run++);
-      //CalculateConcavityPerBoundaryPoint(TesselStruct);
-      LOG(1, "INFO: Writing " << run << "th tesselation file.");
-      StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, dummy);
-    }
-
-    // second step: flip baselines, i.e. add general tetraeder at concave lines
-    // when the tetraeder does not intersect with other already present triangles
-    LineRunner = TesselStruct->LinesOnBoundary.begin();
-    LineAdvance = LineRunner;  // we need an advanced line, as the LineRunner might get removed
-    std::map<double, std::pair<BoundaryLineSet *, double> > GainMap;
-    while (LineRunner != TesselStruct->LinesOnBoundary.end()) {
-      LineAdvance++;
-      line = LineRunner->second;
-      if (!line->CheckConvexityCriterion()) {
-        LOG(2, "INFO: concave line is " << *line << ".");
-        // gather the other points
-        BoundaryPointSet *BPS[4];
-        int m = 0;
-        {
-          for (TriangleMap::iterator runner = line->triangles.begin(); runner != line->triangles.end(); runner++)
-            for (int j = 0; j < 3; j++) // all of their endpoints and baselines
-              if (!line->ContainsBoundaryPoint(runner->second->endpoints[j])) // and neither of its endpoints
-                BPS[m++] = runner->second->endpoints[j];
-        }
-        BPS[2] = line->endpoints[0];
-        BPS[3] = line->endpoints[1];
-        LOG(3, "DEBUG: other line would consist of " << *BPS[0] << " and "
-            << *BPS[1] << ".");
-
-        // check for already present (third) side of the tetraeder as we then
-        // would create a degenerate triangle
-        bool TetraederSidePresent = false;
-        {
-          class TesselPoint *TriangleCandidates[3];
-          TriangleCandidates[0] = BPS[0]->node;
-          TriangleCandidates[1] = BPS[1]->node;
-          TriangleCandidates[2] = BPS[2]->node;
-          if ((TesselStruct->GetPresentTriangle(TriangleCandidates) != NULL)) {
-            LOG(2, "REJECT: Triangle side " << *TriangleCandidates[0] << ","
-                << *TriangleCandidates[1] << "," << *TriangleCandidates[2] << " present.");
-            TetraederSidePresent = true;
-          }
-          TriangleCandidates[2] = BPS[3]->node;
-          if ((TesselStruct->GetPresentTriangle(TriangleCandidates) != NULL)) {
-            LOG(2, "REJECT: Triangle side " << *TriangleCandidates[0] << ","
-                << *TriangleCandidates[1] << "," << *TriangleCandidates[2] << " present.");
-            TetraederSidePresent = true;
-          }
-        }
-
-        if ((BPS[0] != BPS[1]) && (m == 2) && (!TetraederSidePresent)) {
-          // check whether all adjacent triangles do not intersect with new line
-          bool no_line_intersects = true;
-          Vector Intersection;
-          TriangleSet triangles;
-          TriangleSet *firsttriangles = TesselStruct->GetAllTriangles(line->endpoints[0]);
-          TriangleSet *secondtriangles = TesselStruct->GetAllTriangles(line->endpoints[1]);
-          triangles.insert( firsttriangles->begin(), firsttriangles->end() );
-          triangles.insert( secondtriangles->begin(), secondtriangles->end() );
-          delete firsttriangles;
-          delete secondtriangles;
-          for (TriangleSet::const_iterator triangleiter = triangles.begin();
-              triangleiter != triangles.end(); ++triangleiter) {
-            const BoundaryTriangleSet * triangle = *triangleiter;
-            bool line_intersects = triangle->GetIntersectionInsideTriangle(
-                BPS[0]->node->getPosition(),
-                BPS[1]->node->getPosition(),
-                Intersection);
-            // switch result when coinciding with endpoint
-            bool concave_adjacent_line = false;
-            bool intersection_is_endnode = false;
-            for (int j=0;j<2;++j) {
-              if (Intersection.DistanceSquared(BPS[j]->node->getPosition()) < MYEPSILON) {
-                intersection_is_endnode = true;
-                // check whether its an adjacent triangle and if it's concavely connected
-                // only then are we in danger of cutting through it and need to check
-                // sign of normal vector and intersecting line
-                for (int i=0;i<2;++i)
-                  for (int lineindex=0;lineindex < 3;++lineindex)
-                    if ((triangle->lines[lineindex]->ContainsBoundaryPoint(line->endpoints[i]))
-                        && (triangle->lines[lineindex]->ContainsBoundaryPoint(BPS[j]))) {
-                      concave_adjacent_line = !triangle->lines[lineindex]->CheckConvexityCriterion();
-                  }
-                if (concave_adjacent_line) {
-                  const Vector intersector =
-                      BPS[(j+1)%2]->node->getPosition() - Intersection;
-                  if (triangle->NormalVector.ScalarProduct(intersector) >= -MYEPSILON) {
-                    LOG(4, "ACCEPT: Intersection coincides with first endpoint "
-                        << *BPS[j] << ".");
-                    line_intersects = false;
-                  } else {
-                    LOG(4, "REJECT: Intersection ends on wrong side of triangle.");
-                  }
-                } else {
-                  LOG(4, "ACCEPT: Intersection coincides with first endpoint "
-                      << *BPS[j] << ".");
-                  line_intersects = false;
-                }
-              }
-            }
-            // if we have an intersection, check that it is within either
-            // endpoint, i.e. check that scalar product between vectors going
-            // from intersction to either endpoint has negative sign (both
-            // vectors point in opposite directions)
-            if (!intersection_is_endnode && line_intersects) {
-              const Vector firstvector = BPS[0]->node->getPosition() - Intersection;
-              const Vector secondvector = BPS[1]->node->getPosition() - Intersection;
-              if (firstvector.ScalarProduct(secondvector) >= 0)
-                line_intersects = false;
-            }
-            no_line_intersects &= !line_intersects;
-          }
-
-          if (no_line_intersects) {
-            // calculate the volume
-            const double tmp = line->CalculateConvexity();
-            const double gain =
-              CalculateVolumeofGeneralTetraeder(
-                  BPS[0]->node->getPosition(),
-                  BPS[1]->node->getPosition(),
-                  BPS[2]->node->getPosition(),
-                  BPS[3]->node->getPosition());
-
-            GainMap.insert(std::make_pair(tmp, std::make_pair(line,gain) ));
-            LOG(2, "DEBUG: Adding concave line " << *line << " with gain of "
-                << gain << ".");
-          } else {
-            // if 2 or 3 don't
-            LOG(2, "DEBUG: We don't added concave line " << *line
-                << " as other line intersects with adjacent triangles.");
-          }
-        }
-      }
       LineRunner = LineAdvance;
     }
-    // flip line with most gain
-    if (!GainMap.empty()) {
-      line = GainMap.begin()->second.first;
-      const double tmp = GainMap.begin()->second.second;
-      volume += tmp;
-
-//      GainMap.clear();
-
-      // and flip the line
-      LOG(1, "INFO: Flipping current most concave line " << *line << " with gain of "
-          << tmp << ".");
-      TesselStruct->FlipBaseline(line);
-      Concavity = true;
-    }
-  } while ((Concavity)); // && (run < 100)
+    run++;
+  } while (Concavity);
+  //CalculateConcavityPerBoundaryPoint(TesselStruct);
+  //StoreTrianglesinFile(mol, filename, "-third");
+
+  // third step: IsConvexRectangle
+//  LineRunner = TesselStruct->LinesOnBoundary.begin();
+//  LineAdvance = LineRunner;  // we need an advanced line, as the LineRunner might get removed
+//  while (LineRunner != TesselStruct->LinesOnBoundary.end()) {
+//    LineAdvance++;
+//    line = LineRunner->second;
+//    LOG(1, "INFO: Current line is " << *line << ".");
+//    //if (LineAdvance != TesselStruct->LinesOnBoundary.end())
+//      //LOG(1, "INFO: Next line will be " << *(LineAdvance->second) << ".");
+//    if (!line->CheckConvexityCriterion(out)) {
+//      LOG(1, "INFO: ... line " << *line << " is concave, flipping it.");
+//
+//      // take highest of both lines
+//      point = TesselStruct->IsConvexRectangle(line);
+//      if (point != NULL)
+//        volume += TesselStruct->RemovePointFromTesselatedSurface(point);
+//    }
+//    LineRunner = LineAdvance;
+//  }
 
   CalculateConcavityPerBoundaryPoint(TesselStruct);
@@ -738 +615 @@
 
   // end
-  LOG(0, "RESULT: Added volume in convexization is " << volume << ".");
+  LOG(0, "Volume is " << volume << ".");
   return volume;
 };
@@ -889 +766 @@
 //  status = CheckListOfBaselines(TesselStruct);
 
-//  cout << "before correction" << endl;
+  cout << "before correction" << endl;
 
   // store before correction
@@ -902 +779 @@
   // write final envelope
   CalculateConcavityPerBoundaryPoint(TesselStruct);
-//  cout << "after correction" << endl;
+  cout << "after correction" << endl;
   StoreTrianglesinFile(mol, TesselStruct, filename, "");