Changeset d289c6 for src/Tesselation/boundary.cpp

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:

81ffd8

Parents:

6da9e9

git-author:

Frederik Heber <heber@…> (10/04/14 14:34:41)

git-committer:

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

Message:

Rewrote ConvexizeNonconvexEnvelope().

• now works for some sphere radii on amylose (e.g. 40 and 50).
• we first remove all points that are only attached to concave lines.
• also, we always flip the baseline (i.e. add a general tetraeder to the non-convex volume) for concave lines when the tetraeder does not intersect with present triangles.

File:

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

src/Tesselation/boundary.cpp

@@ -538 +538 @@
       << " convex ...");
 
-  // First step: RemovePointFromTesselatedSurface
   do {
     Concavity = false;
-    sprintf(dummy, "-first-%d", run);
+
+    sprintf(dummy, "-%d", run++);
     //CalculateConcavityPerBoundaryPoint(TesselStruct);
     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
@@ -551 +552 @@
       point = PointRunner->second;
       LOG(2, "INFO: Current point is " << *point << ".");
+      bool all_lines_concave = true;
       for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
         line = LineRunner->second;
         LOG(3, "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;
-        }
+        all_lines_concave &= !line->CheckConvexityCriterion();
+      }
+      if (all_lines_concave) {
+        // remove the point if needed
+        LOG(1, "... point " << *point << " cannot be on convex envelope, all lines concave.");
+        volume += TesselStruct->RemovePointFromTesselatedSurface(point);
+        Concavity = true;
       }
       PointRunner = PointAdvance;
     }
 
-    sprintf(dummy, "-second-%d", run);
+    sprintf(dummy, "-%d", run++);
     //CalculateConcavityPerBoundaryPoint(TesselStruct);
     StoreTrianglesinFile(mol, (const Tesselation *&)TesselStruct, filename, dummy);
 
-    // second step: PickFarthestofTwoBaselines
+    // 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, BoundaryLineSet *> 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 =
+                CalculateVolumeofGeneralTetraeder(
+                    BPS[0]->node->getPosition(),
+                    BPS[1]->node->getPosition(),
+                    BPS[2]->node->getPosition(),
+                    BPS[3]->node->getPosition());
+            GainMap.insert(std::make_pair(tmp, line));
+            LOG(2, "DEBUG: Adding concave line " << *line << " with gain of "
+                << tmp << ".");
+          } 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.rbegin()->second;
+      const double tmp = GainMap.rbegin()->first;
+      volume += tmp;
+//      GainMap.clear();
+
+      // and flip the line
+      LOG(1, "INFO: Flipping concave line " << *line << " with gain of "
+          << tmp << ".");
+      TesselStruct->FlipBaseline(line);
+      Concavity = true;
+    }
+  } while ((Concavity)); // && (run < 100)
 
   CalculateConcavityPerBoundaryPoint(TesselStruct);
@@ -618 +718 @@
 
   // end
-  LOG(0, "Volume is " << volume << ".");
+  LOG(0, "RESULT: Added volume in convexization is " << volume << ".");
   return volume;
 };