Changeset 30ce7d for src

Timestamp:

Oct 19, 2014, 5:13:11 PM (10 years ago)

Author:

Frederik Heber <heber@…>

Branches:

Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable

Children:

6da9e9

Parents:

77eec5

git-author:

Frederik Heber <heber@…> (10/09/14 07:37:06)

git-committer:

Frederik Heber <heber@…> (10/19/14 17:13:11)

Message:

Extracted RemovePointSurroundedByPolygon() out of Tesselation::RemovePointFromTesselatedSurface()

• also added new functions Tesselation::RemoveFullConcavePointFromTesselatedSurface() and Tesselation::CheckAllConcaveInPolygon().

Location:

src/Tesselation

Files:

• tesselation.cpp (modified) (5 diffs)
• tesselation.hpp (modified) (1 diff)

src/Tesselation/tesselation.cpp

@@ -3330 +3330 @@
 }
 
+double Tesselation::RemovePointSurroundedByPolygon(
+    TesselPointList *connectedPath,
+    BoundaryPointSet *point)
+{
+  double volume = 0.;
+  const Vector OldPoint = point->node->getPosition();
+
+  TesselPoint *oldNode = point->node;
+  // remove present triangles for this connectedPath
+  unsigned int count = 0;
+  for (TesselPointList::const_iterator iter = connectedPath->begin();
+      iter != connectedPath->end(); ++iter)
+    LOG(4, "DEBUG: Node in connectedPath is " << **iter);
+
+  {
+    TesselPointList::iterator FirstNode, SecondNode;
+    SecondNode = connectedPath->begin();
+    FirstNode = SecondNode++;
+    for (;FirstNode != connectedPath->end(); ++SecondNode, ++FirstNode) {
+      LOG(3, "DEBUG: MiddleNode is " << **FirstNode << ".");
+      if (SecondNode == connectedPath->end())
+        SecondNode = connectedPath->begin();
+      TesselPoint *TriangleCandidates[3];
+      TriangleCandidates[0] = *FirstNode;
+      TriangleCandidates[1] = *SecondNode;
+      TriangleCandidates[2] = oldNode;
+      BoundaryTriangleSet *triangle = GetPresentTriangle(TriangleCandidates);
+      ASSERT( triangle != NULL,
+          "Tesselation::RemovePointSurroundedByPolygon() - triangle to points "
+          +toString((**SecondNode))+", "+toString((**FirstNode))+" and "
+          +toString(*oldNode)+" does not exist.");
+      LOG(3, "DEBUG: Removing triangle " << *triangle << ".");
+      RemoveTesselationTriangle(triangle);
+      ++count;
+    }
+    LOG(2, "INFO: " << count << " triangles were removed.");
+  }
+
+  // re-create all triangles by going through connected points list
+  LineList NewLines;
+  typedef std::vector<double> angles_t;
+  angles_t angles;
+  count = 0;
+  for (; !connectedPath->empty();) {
+    // search middle node with widest angle to next neighbours
+    TesselPointList::iterator StartNode, MiddleNode, EndNode;
+    for (MiddleNode = connectedPath->begin(); MiddleNode != connectedPath->end(); MiddleNode++) {
+      LOG(3, "INFO: MiddleNode is " << **MiddleNode << ".");
+      // construct vectors to next and previous neighbour
+      StartNode = MiddleNode;
+      if (StartNode == connectedPath->begin())
+        StartNode = connectedPath->end();
+      StartNode--;
+      //LOG(3, "INFO: StartNode is " << **StartNode << ".");
+      const Vector Point = ((*StartNode)->getPosition()) - ((*MiddleNode)->getPosition());
+      EndNode = MiddleNode;
+      EndNode++;
+      if (EndNode == connectedPath->end())
+        EndNode = connectedPath->begin();
+      //LOG(3, "INFO: EndNode is " << **StartNode << ".");
+      const Vector Reference = ((*EndNode)->getPosition()) - ((*MiddleNode)->getPosition());
+      Vector OrthogonalVector = ((*MiddleNode)->getPosition()) - OldPoint;
+      OrthogonalVector.MakeNormalTo(Reference);
+      angles.push_back(GetAngle(Point, Reference, OrthogonalVector));
+    }
+    ASSERT( !angles.empty(),
+        "Tesselation::RemovePointSurroundedByPolygon() - angles empty");
+    const angles_t::const_iterator maxiter = std::max_element(angles.begin(), angles.end());
+    angles_t::const_iterator miniter = angles.begin();
+    // distinguish between convex and nonconvex polygon
+    if (*maxiter > M_PI) {
+      // connectedPath is not convex: The idea is to fill any kinks pointing
+      // inside into the connectedPath close to this concave spot first, making
+      // it eventually become convex.
+      // Hence, use adjacent (and convex) fill-in point
+      miniter = (maxiter == angles.begin()) ? angles.end()-1 : maxiter-1;
+      if (*miniter > M_PI) {
+        miniter = (maxiter+1 == angles.end()) ? angles.begin() : maxiter+1;
+        if (*miniter > M_PI) {
+          miniter = std::min_element(angles.begin(), angles.end());
+        }
+      }
+    } else {
+      // is convex
+      miniter = std::min_element(angles.begin(), angles.end(), CloserToPiHalf());
+    }
+    MiddleNode = connectedPath->begin();
+    std::advance(MiddleNode, std::distance(const_cast<const angles_t &>(angles).begin(), miniter));
+
+    ASSERT(MiddleNode != connectedPath->end(),
+        "Tesselation::RemovePointSurroundedByPolygon() - Could not find a smallest angle!");
+    angles.clear();
+    StartNode = MiddleNode;
+    EndNode = MiddleNode;
+    if (StartNode == connectedPath->begin())
+      StartNode = connectedPath->end();
+    StartNode--;
+    EndNode++;
+    if (EndNode == connectedPath->end())
+      EndNode = connectedPath->begin();
+    LOG(2, "INFO: StartNode is " << **StartNode << ".");
+    LOG(2, "INFO: MiddleNode is " << **MiddleNode << ".");
+    LOG(2, "INFO: EndNode is " << **EndNode << ".");
+    LOG(1, "INFO: Attempting to create triangle " << (*StartNode)->getName() << ", " << (*MiddleNode)->getName() << " and " << (*EndNode)->getName() << ".");
+    TesselPoint *TriangleCandidates[3];
+    TriangleCandidates[0] = *StartNode;
+    TriangleCandidates[1] = *MiddleNode;
+    TriangleCandidates[2] = *EndNode;
+    BoundaryTriangleSet *triangle = GetPresentTriangle(TriangleCandidates);
+    if (triangle != NULL) {
+      const Vector center = 1./3. * ((*StartNode)->getPosition()
+              + (*MiddleNode)->getPosition()
+              + (*EndNode)->getPosition());
+      const Vector NormalVector = OldPoint - center;
+      // check orientation of normal vector (that points inside)
+      ASSERT( triangle->NormalVector.ScalarProduct(NormalVector) > std::numeric_limits<double>::epsilon()*1e2,
+          "Tesselation::RemovePointSurroundedByPolygon() - New triangle with same orientation already present as "
+          +toString(*triangle)+"!");
+    }
+
+    LOG(3, "DEBUG: Adding new triangle points.");
+    AddTesselationPoint(*StartNode, 0);
+    AddTesselationPoint(*MiddleNode, 1);
+    AddTesselationPoint(*EndNode, 2);
+    LOG(3, "DEBUG: Adding new triangle lines.");
+    AddTesselationLine(NULL, NULL, TPS[0], TPS[1], 0);
+    // line between start and end must be new (except for very last triangle)
+    if (AddTesselationLine(NULL, NULL, TPS[0], TPS[2], 1))
+      NewLines.push_back(BLS[1]);
+    AddTesselationLine(NULL, NULL, TPS[1], TPS[2], 2);
+    BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
+    const Vector center = 1./3. * ((*StartNode)->getPosition()
+            + (*MiddleNode)->getPosition()
+            + (*EndNode)->getPosition());
+    const Vector NormalVector = OldPoint - center;
+    BTS->GetNormalVector(NormalVector);
+    AddTesselationTriangle();
+    // calculate volume summand as a general tetraeder
+    volume += CalculateVolumeofGeneralTetraeder(
+        TPS[0]->node->getPosition(),
+        TPS[1]->node->getPosition(),
+        TPS[2]->node->getPosition(),
+        OldPoint);
+    // advance number
+    ++count;
+
+    // prepare nodes for next triangle
+    LOG(3, "DEBUG: Removing " << **MiddleNode << " from closed path.");
+    connectedPath->remove(*MiddleNode); // remove the middle node (it is surrounded by triangles)
+    LOG(3, "DEBUG: Remaining points: " << connectedPath->size() << ".");
+    ASSERT(connectedPath->size() >= 2,
+        "Tesselation::RemovePointSurroundedByPolygon() - There are only two endpoints left!");
+    if (connectedPath->size() == 2) { // we are done
+      connectedPath->remove(*StartNode); // remove the start node
+      connectedPath->remove(*EndNode); // remove the end node
+    }
+  }
+  LOG(1, "INFO: " << count << " triangles were created.");
+
+  return volume;
+}
+
 /** Removes a boundary point from the envelope while keeping it closed.
  * We remove the old triangles connected to the point and re-create new triangles to close the surface following this ansatz:
@@ -3343 +3505 @@
 double Tesselation::RemovePointFromTesselatedSurface(class BoundaryPointSet *point)
 {
-  class BoundaryLineSet *line = NULL;
-  class BoundaryTriangleSet *triangle = NULL;
-  Vector OldPoint; //, NormalVector;
   double volume = 0;
-  int count = 0;
 
   if (point == NULL) {
@@ -3355 +3513 @@
     LOG(4, "DEBUG: Removing point " << *point << " from tesselated boundary ...");
 
-  // copy old location for the volume
-  OldPoint = (point->node->getPosition());
-
   // get list of connected points
   if (point->lines.empty()) {
@@ -3365 +3520 @@
 
   list<TesselPointList *> *ListOfClosedPaths = GetClosedPathsOfConnectedPoints(point->node);
-  TesselPointList *connectedPath = NULL;
-
-  // gather all triangles
-  for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++)
-    count += LineRunner->second->triangles.size();
-  TriangleMap Candidates;
-  for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
-    line = LineRunner->second;
-    for (TriangleMap::iterator TriangleRunner = line->triangles.begin(); TriangleRunner != line->triangles.end(); TriangleRunner++) {
-      triangle = TriangleRunner->second;
-      Candidates.insert(TrianglePair(triangle->Nr, triangle));
-    }
-  }
-
-  // remove all triangles
-  count = 0;
-//  NormalVector.Zero();
-  for (TriangleMap::iterator Runner = Candidates.begin(); Runner != Candidates.end(); Runner++) {
-    LOG(3, "DEBUG: Removing triangle " << *(Runner->second) << ".");
-//    NormalVector -= Runner->second->NormalVector; // has to point inward
-    RemoveTesselationTriangle(Runner->second);
-    count++;
-  }
-  LOG(2, "INFO: " << count << " triangles were removed.");
-
+  list<TesselPointList *>::iterator ListAdvance = ListOfClosedPaths->begin();
+  list<TesselPointList *>::iterator ListRunner = ListAdvance;
+//  TriangleMap::iterator NumberRunner = Candidates.begin();
+  Vector Point, Reference, OrthogonalVector;
+  for (; ListRunner != ListOfClosedPaths->end(); ListRunner = ListAdvance) { // go through all closed paths
+    if (ListAdvance != ListOfClosedPaths->end())
+      ListAdvance++;
+
+    TesselPointList *connectedPath = *ListRunner;
+
+    volume += RemovePointSurroundedByPolygon(connectedPath, point);
+
+    ListOfClosedPaths->remove(connectedPath);
+    delete (connectedPath);
+  }
+  delete (ListOfClosedPaths);
+
+  LOG(1, "INFO: Removed volume is " << volume << ".");
+
+  return volume;
+}
+;
+
+bool Tesselation::CheckAllConcaveInPolygon(
+    const TesselPointList *connectedPath,
+    const BoundaryPointSet *point
+    )
+{
+  // check whether lines in this path to point to remove are all concave
+  bool all_lines_concave = true;
+  if (connectedPath->size() >= 2) {
+    TesselPointList::const_iterator StartNode, MiddleNode, EndNode;
+    // have nearest neighbors to a middle node to know adjacent triangles
+    for (MiddleNode = connectedPath->begin();
+        all_lines_concave && (MiddleNode != connectedPath->end());
+        MiddleNode++) {
+      LOG(3, "INFO: MiddleNode is " << **MiddleNode << ".");
+      EndNode = MiddleNode;
+      if (EndNode == connectedPath->begin())
+        EndNode = connectedPath->end();
+      --EndNode;
+      StartNode = MiddleNode;
+      ++StartNode;
+      if (StartNode == connectedPath->end())
+        StartNode = connectedPath->begin();
+
+      AddTesselationPoint(point->node, 0);
+      AddTesselationPoint(*MiddleNode, 1);
+
+      ASSERT( point->lines.find((*MiddleNode)->getNr()) != point->lines.end(),
+          "Tesselation::CheckAllConcaveInPolygon() - MiddleNode "
+          +toString(**MiddleNode)+" not present in "
+          +toString(*point)+"'s lines.");
+
+      // get line between Node and point and check
+      std::pair<LineMap::const_iterator, LineMap::const_iterator> FindPair =
+          point->lines.equal_range((*MiddleNode)->getNr());
+      LineMap::const_iterator FindLine = FindPair.first;
+      for (; FindLine != FindPair.second; ++FindLine) {
+        // line  has got two triangles, check whether they resemble those
+        // with start and endnode
+        const BoundaryLineSet *currentline = FindLine->second;
+        unsigned int matching_triangles = 0;
+        for (TriangleMap::const_iterator triangleiter = currentline->triangles.begin();
+            triangleiter != currentline->triangles.end(); ++triangleiter) {
+          const BoundaryTriangleSet *triangle = triangleiter->second;
+          AddTesselationPoint(*StartNode, 2);
+          if (triangle->IsPresentTupel(TPS))
+            ++matching_triangles;
+          AddTesselationPoint(*EndNode, 2);
+          if (triangle->IsPresentTupel(TPS))
+            ++matching_triangles;
+        }
+        if (matching_triangles == 2)
+          break;
+      }
+      if (FindLine != FindPair.second) {
+        LOG(3, "INFO: Current line of point " << *point << " is " << *FindLine->second << ".");
+        all_lines_concave &= !FindLine->second->CheckConvexityCriterion();
+      }
+    }
+  } else {
+    // check the single line
+    if (connectedPath->empty())
+      return false;
+    LineMap::const_iterator FindLine = point->lines.find((*connectedPath->begin())->getNr());
+    ASSERT( FindLine != point->lines.end(),
+        "Tesselation::CheckAllConcaveInPolygon() - point "
+        +toString((*connectedPath->begin())->getNr())+" not present in "
+        +toString(*point)+"'s lines.");
+    return !FindLine->second->CheckConvexityCriterion();
+  }
+
+  return all_lines_concave;
+}
+
+/** Removes a boundary point from the envelope while keeping it closed.
+ * We remove the old triangles connected to the point and re-create new triangles to close the surface following this ansatz:
+ *  -# a closed path(s) of boundary points surrounding the point to be removed is constructed
+ *  -# on each closed path, we pick three adjacent points, create a triangle with them and subtract the middle point from the path
+ *  -# we advance two points (i.e. the next triangle will start at the ending point of the last triangle) and continue as before
+ *  -# the surface is closed, when the path is empty
+ * Thereby, we (hopefully) make sure that the removed points remains beneath the surface (this is checked via IsInnerPoint eventually).
+ * \param *out output stream for debugging
+ * \param *point point to be removed
+ * \return volume added to the volume inside the tesselated surface by the removal
+ */
+double Tesselation::RemoveFullConcavePointFromTesselatedSurface(class BoundaryPointSet *point)
+{
+  double volume = 0;
+
+  if (point == NULL) {
+    ELOG(1, "Cannot remove the point " << point << ", it's NULL!");
+    return 0.;
+  } else
+    LOG(4, "DEBUG: Removing point " << *point << " from tesselated boundary ...");
+
+  // get list of connected points
+  if (point->lines.empty()) {
+    ELOG(1, "Cannot remove the point " << *point << ", it's connected to no lines!");
+    return 0.;
+  }
+
+  list<TesselPointList *> *ListOfClosedPaths = GetClosedPathsOfConnectedPoints(point->node);
   list<TesselPointList *>::iterator ListAdvance = ListOfClosedPaths->begin();
   list<TesselPointList *>::iterator ListRunner = ListAdvance;
@@ -3395 +3649 @@
   TesselPointList::iterator StartNode, MiddleNode, EndNode;
   Vector Point, Reference, OrthogonalVector;
-  if (count > 2) { // less than three triangles, then nothing will be created
-    class TesselPoint *TriangleCandidates[3];
-    count = 0;
-    for (; ListRunner != ListOfClosedPaths->end(); ListRunner = ListAdvance) { // go through all closed paths
-      if (ListAdvance != ListOfClosedPaths->end())
-        ListAdvance++;
-
-      connectedPath = *ListRunner;
-      // re-create all triangles by going through connected points list
-      LineList NewLines;
-      typedef std::vector<double> angles_t;
-      angles_t angles;
-      for (; !connectedPath->empty();) {
-        // search middle node with widest angle to next neighbours
-        EndNode = connectedPath->end();
-        for (MiddleNode = connectedPath->begin(); MiddleNode != connectedPath->end(); MiddleNode++) {
-          LOG(3, "INFO: MiddleNode is " << **MiddleNode << ".");
-          // construct vectors to next and previous neighbour
-          StartNode = MiddleNode;
-          if (StartNode == connectedPath->begin())
-            StartNode = connectedPath->end();
-          StartNode--;
-          //LOG(3, "INFO: StartNode is " << **StartNode << ".");
-          Point = ((*StartNode)->getPosition()) - ((*MiddleNode)->getPosition());
-          StartNode = MiddleNode;
-          StartNode++;
-          if (StartNode == connectedPath->end())
-            StartNode = connectedPath->begin();
-          //LOG(3, "INFO: EndNode is " << **StartNode << ".");
-          Reference = ((*StartNode)->getPosition()) - ((*MiddleNode)->getPosition());
-          OrthogonalVector = ((*MiddleNode)->getPosition()) - OldPoint;
-          OrthogonalVector.MakeNormalTo(Reference);
-          angles.push_back(GetAngle(Point, Reference, OrthogonalVector));
-        }
-        const angles_t::const_iterator maxiter = std::max_element(angles.begin(), angles.end());
-        angles_t::const_iterator miniter = angles.begin();
-        // distinguish between convex and nonconvex polygon
-        if (*maxiter > M_PI) {
-          // connectedPath is not convex
-          // use adjacent (and convx) fill-in point
-          miniter = (maxiter == angles.begin()) ? angles.end()-1 : maxiter-1;
-          if (*miniter > M_PI) {
-            miniter = (maxiter+1 == angles.end()) ? angles.begin() : maxiter+1;
-            if (*(++miniter) > M_PI) {
-              miniter = std::min_element(angles.begin(), angles.end());
-            }
-          }
-        } else {
-          // is convex
-          miniter = std::min_element(angles.begin(), angles.end(), CloserToPiHalf());
-        }
-        MiddleNode = connectedPath->begin();
-        std::advance(MiddleNode, std::distance(const_cast<const angles_t &>(angles).begin(), miniter));
-        angles.clear();
-
-        ASSERT(MiddleNode != connectedPath->end(), 
-            "Tesselation::RemovePointFromTesselatedSurface() - Could not find a smallest angle!");
-        StartNode = MiddleNode;
-        EndNode = MiddleNode;
-        if (StartNode == connectedPath->begin())
-          StartNode = connectedPath->end();
-        StartNode--;
-        EndNode++;
-        if (EndNode == connectedPath->end())
-          EndNode = connectedPath->begin();
-        LOG(2, "INFO: StartNode is " << **StartNode << ".");
-        LOG(2, "INFO: MiddleNode is " << **MiddleNode << ".");
-        LOG(2, "INFO: EndNode is " << **EndNode << ".");
-        LOG(1, "INFO: Attempting to create triangle " << (*StartNode)->getName() << ", " << (*MiddleNode)->getName() << " and " << (*EndNode)->getName() << ".");
-        TriangleCandidates[0] = *StartNode;
-        TriangleCandidates[1] = *MiddleNode;
-        TriangleCandidates[2] = *EndNode;
-        triangle = GetPresentTriangle(TriangleCandidates);
-        if (triangle != NULL) {
-          const Vector center = 1./3. * ((*StartNode)->getPosition()
-                  + (*MiddleNode)->getPosition()
-                  + (*EndNode)->getPosition());
-          const Vector NormalVector = OldPoint - center;
-          // check orientation of normal vector (that points inside)
-          ASSERT( triangle->NormalVector.ScalarProduct(NormalVector) > std::numeric_limits<double>::epsilon()*1e2,
-              "Tesselation::RemovePointFromTesselatedSurface() - New triangle with same orientation already present as "
-              +toString(*triangle)+"!");
-        }
-        if (0) {
-          StartNode++;
-          MiddleNode++;
-          EndNode++;
-          if (StartNode == connectedPath->end())
-            StartNode = connectedPath->begin();
-          if (MiddleNode == connectedPath->end())
-            MiddleNode = connectedPath->begin();
-          if (EndNode == connectedPath->end())
-            EndNode = connectedPath->begin();
-          continue;
-        }
-        LOG(3, "DEBUG: Adding new triangle points.");
-        AddTesselationPoint(*StartNode, 0);
-        AddTesselationPoint(*MiddleNode, 1);
-        AddTesselationPoint(*EndNode, 2);
-        LOG(3, "DEBUG: Adding new triangle lines.");
-        AddTesselationLine(NULL, NULL, TPS[0], TPS[1], 0);
-        // line between start and end must be new (except for very last triangle)
-        if (AddTesselationLine(NULL, NULL, TPS[0], TPS[2], 1))
-          NewLines.push_back(BLS[1]);
-        AddTesselationLine(NULL, NULL, TPS[1], TPS[2], 2);
-        BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);
-        const Vector center = 1./3. * ((*StartNode)->getPosition()
-                + (*MiddleNode)->getPosition()
-                + (*EndNode)->getPosition());
-        const Vector NormalVector = OldPoint - center;
-        BTS->GetNormalVector(NormalVector);
-        AddTesselationTriangle();
-        // calculate volume summand as a general tetraeder
-        volume += CalculateVolumeofGeneralTetraeder(TPS[0]->node->getPosition(), TPS[1]->node->getPosition(), TPS[2]->node->getPosition(), OldPoint);
-        // advance number
-        count++;
-
-        // prepare nodes for next triangle
-        LOG(3, "DEBUG: Removing " << **MiddleNode << " from closed path.");
-        connectedPath->remove(*MiddleNode); // remove the middle node (it is surrounded by triangles)
-        LOG(3, "DEBUG: Remaining points: " << connectedPath->size() << ".");
-        ASSERT(connectedPath->size() >= 2, "Tesselation::RemovePointFromTesselatedSurface() - There are only two endpoints left!");
-        if (connectedPath->size() == 2) { // we are done
-          connectedPath->remove(*StartNode); // remove the start node
-          connectedPath->remove(*EndNode); // remove the end node
-          break;
-        } else {
-          StartNode = EndNode;
-          MiddleNode = StartNode;
-          MiddleNode++;
-          if (MiddleNode == connectedPath->end())
-            MiddleNode = connectedPath->begin();
-          EndNode = MiddleNode;
-          EndNode++;
-          if (EndNode == connectedPath->end())
-            EndNode = connectedPath->begin();
-        }
-      }
-
-      ListOfClosedPaths->remove(connectedPath);
-      delete (connectedPath);
-    }
-    LOG(1, "INFO: " << count << " triangles were created.");
-  } else {
-    while (!ListOfClosedPaths->empty()) {
-      ListRunner = ListOfClosedPaths->begin();
-      connectedPath = *ListRunner;
-      ListOfClosedPaths->remove(connectedPath);
-      delete (connectedPath);
-    }
-    LOG(3, "DEBUG: No need to create any triangles.");
+  for (; ListRunner != ListOfClosedPaths->end(); ListRunner = ListAdvance) { // go through all closed paths
+    if (ListAdvance != ListOfClosedPaths->end())
+      ListAdvance++;
+
+    TesselPointList *connectedPath = *ListRunner;
+
+    if (CheckAllConcaveInPolygon(connectedPath, point)) {
+      LOG(1, "INFO: ... point " << *point << " cannot be on convex envelope, all lines concave.");
+      volume += RemovePointSurroundedByPolygon(connectedPath, point);
+    }
+
+    ListOfClosedPaths->remove(connectedPath);
+    delete (connectedPath);
   }
   delete (ListOfClosedPaths);
 
-  LOG(1, "INFO: Removed volume is " << volume << ".");
+  if (volume > 0.)
+    LOG(1, "INFO: Removed volume is " << volume << ".");
 
   return volume;

src/Tesselation/tesselation.hpp

@@ -103 +103 @@
     void GuessStartingTriangle();
     bool InsertStraddlingPoints(IPointCloud & cloud, const LinkedCell_deprecated *LC);
+    double RemovePointSurroundedByPolygon(
+        TesselPointList *connectedPath,
+        BoundaryPointSet *point);
+    bool CheckAllConcaveInPolygon(
+        const TesselPointList *connectedPath,
+        const BoundaryPointSet *point
+        );
+    double RemoveFullConcavePointFromTesselatedSurface(class BoundaryPointSet *point);
     double RemovePointFromTesselatedSurface(class BoundaryPointSet *point);
     class BoundaryLineSet * FlipBaseline(class BoundaryLineSet *Base);