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

File:

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

src/Tesselation/boundary.cpp

@@ -513 +513 @@
  * \return volume difference between the non- and the created convex envelope
  */
-double ConvexizeNonconvexEnvelope(class Tesselation *&TesselStruct, const molecule * const mol, const char * const filename)
+double ConvexizeNonconvexEnvelope(
+    Tesselation *&TesselStruct,
+    const molecule * const mol,
+    const char * const filename,
+    bool DebugOutputEveryStep)
 {
         //Info FunctionInfo(__func__);
@@ -535 +539 @@
   }
 
-  // First step: RemovePointFromTesselatedSurface
+  LOG(1, "INFO: Making tesselated surface with " << TesselStruct->TrianglesOnBoundaryCount
+      << " convex ...");
+
+  // first purge all degenerate triangles
+  TesselStruct->RemoveDegeneratedTriangles();
+
   do {
     Concavity = false;
-    sprintf(dummy, "-first-%d", run);
-    //CalculateConcavityPerBoundaryPoint(TesselStruct);
-    StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, dummy);
-
+
+    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
     PointRunner = TesselStruct->PointsOnBoundary.begin();
     PointAdvance = PointRunner; // we need an advanced point, as the PointRunner might get removed
@@ -547 +561 @@
       PointAdvance++;
       point = PointRunner->second;
-      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);
+      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())
+          break;
+      }
+      // remove the point if needed
+      if (LineRunner != point->lines.end()) {
+        const double tmp = TesselStruct->RemoveFullConcavePointFromTesselatedSurface(point);
+        if (tmp > 0.) {
+          volume += tmp;
           Concavity = true;
-          break;
         }
       }
@@ -564 +581 @@
     }
 
-    sprintf(dummy, "-second-%d", run);
-    //CalculateConcavityPerBoundaryPoint(TesselStruct);
-    StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, dummy);
-
-    // second step: PickFarthestofTwoBaselines
+    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;
-      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;
+      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;
     }
-    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;
-//  }
+    // 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)
 
   CalculateConcavityPerBoundaryPoint(TesselStruct);
@@ -615 +738 @@
 
   // end
-  LOG(0, "Volume is " << volume << ".");
+  LOG(0, "RESULT: Added volume in convexization is " << volume << ".");
   return volume;
 };
@@ -766 +889 @@
 //  status = CheckListOfBaselines(TesselStruct);
 
-  cout << "before correction" << endl;
+//  cout << "before correction" << endl;
 
   // store before correction
@@ -779 +902 @@
   // write final envelope
   CalculateConcavityPerBoundaryPoint(TesselStruct);
-  cout << "after correction" << endl;
+//  cout << "after correction" << endl;
   StoreTrianglesinFile(mol, TesselStruct, filename, "");