Changeset 75d42a for src/Tesselation/tesselation.cpp

Timestamp:

Oct 19, 2014, 5:13:10 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:

dfe054

Parents:

c2f306

git-author:

Frederik Heber <heber@…> (09/25/14 18:23:50)

git-committer:

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

Message:

FIX: RemovePointFromTesselatedSurface() needs to pick smallest angle when polygon is not convex.

• AddTesselationLine() now returns whether new or present line was used.
• FIX: Last inner line in closed polygon has been added by one but last triangle already, we added it twice.

File:

• src/Tesselation/tesselation.cpp (modified) (8 diffs)

src/Tesselation/tesselation.cpp

@@ -35 +35 @@
 #include "CodePatterns/MemDebug.hpp"
 
+#include <algorithm>
+#include <boost/foreach.hpp>
 #include <fstream>
 #include <iomanip>
+#include <iterator>
 #include <sstream>
 
@@ -731 +734 @@
  * @param *b second endpoint
  * @param n index of Tesselation::BLS giving the line with both endpoints
- */
-void Tesselation::AddTesselationLine(const Vector * const OptCenter, const BoundaryPointSet * const candidate, class BoundaryPointSet *a, class BoundaryPointSet *b, const int n)
+ * @return true - inserted a new line, false - present line used
+ */
+bool Tesselation::AddTesselationLine(const Vector * const OptCenter, const BoundaryPointSet * const candidate, class BoundaryPointSet *a, class BoundaryPointSet *b, const int n)
 {
   bool insertNewLine = true;
@@ -784 +788 @@
     AddExistingTesselationTriangleLine(WinningLine, n);
   }
+
+  return insertNewLine;
 }
 ;
@@ -3198 +3204 @@
 ;
 
+struct CloserToPiHalf
+{
+  bool operator()(double angle, double smallestangle) {
+    return (fabs(angle - M_PI/2.) < fabs(smallestangle - M_PI/2.));
+  }
+};
+
 /** 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:
@@ -3262 +3275 @@
 //  TriangleMap::iterator NumberRunner = Candidates.begin();
   TesselPointList::iterator StartNode, MiddleNode, EndNode;
-  double angle;
-  double smallestangle;
   Vector Point, Reference, OrthogonalVector;
   if (count > 2) { // less than three triangles, then nothing will be created
@@ -3275 +3286 @@
       // 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();
-        smallestangle = 0.;
         for (MiddleNode = connectedPath->begin(); MiddleNode != connectedPath->end(); MiddleNode++) {
           LOG(1, "INFO: MiddleNode is " << **MiddleNode << ".");
@@ -3296 +3308 @@
           OrthogonalVector = ((*MiddleNode)->getPosition()) - OldPoint;
           OrthogonalVector.MakeNormalTo(Reference);
-          angle = GetAngle(Point, Reference, OrthogonalVector);
-          //if (angle < M_PI)  // no wrong-sided triangles, please?
-          if (fabs(angle - M_PI) < fabs(smallestangle - M_PI)) { // get straightest angle (i.e. construct those triangles with smallest area first)
-            smallestangle = angle;
-            EndNode = MiddleNode;
-          }
+          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
+          miniter = std::min_element(angles.begin(), angles.end());
+        } else {
+          // is convex
+          miniter = std::min_element(angles.begin(), angles.end(), CloserToPiHalf());
+        }
+        EndNode = connectedPath->begin();
+        std::advance(EndNode,
+            std::distance(const_cast<const angles_t &>(angles).begin(), miniter));
+        angles.clear();
+
         MiddleNode = EndNode;
-        if (MiddleNode == connectedPath->end()) {
-          ELOG(0, "CRITICAL: Could not find a smallest angle!");
-          performCriticalExit();
-        }
+        ASSERT( MiddleNode == connectedPath->end(),
+          "Tesselation::RemovePointFromTesselatedSurface() - Could not find a smallest angle!");
         StartNode = MiddleNode;
         if (StartNode == connectedPath->begin())
@@ -3347 +3367 @@
         LOG(3, "Adding new triangle lines.");
         AddTesselationLine(NULL, NULL, TPS[0], TPS[1], 0);
-        AddTesselationLine(NULL, NULL, TPS[0], TPS[2], 1);
-        NewLines.push_back(BLS[1]);
+        // 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);