Changes in molecuilder/src/boundary.cpp [d6b011:78dac6]

File:

• molecuilder/src/boundary.cpp (modified) (27 diffs)

molecuilder/src/boundary.cpp

@@ -118 +118 @@
  * \param *mol molecule structure with atom positions
  */
-void WriteVrmlFile(ofstream *out, ofstream *vrmlfile, class Tesselation *Tess, class molecule *mol)
+void write_vrml_file(ofstream *out, ofstream *vrmlfile, class Tesselation *Tess, class molecule *mol)
 {
   atom *Walker = mol->start;
@@ -172 +172 @@
  * \param *mol molecule structure with atom positions
  */
-void WriteRaster3dFile(ofstream *out, ofstream *rasterfile, class Tesselation *Tess, class molecule *mol)
+void write_raster3d_file(ofstream *out, ofstream *rasterfile, class Tesselation *Tess, class molecule *mol)
 {
   atom *Walker = mol->start;
@@ -227 +227 @@
  * \param N arbitrary number to differentiate various zones in the tecplot format
  */
-void WriteTecplotFile(ofstream *out, ofstream *tecplot, class Tesselation *TesselStruct, class molecule *mol, int N)
+void write_tecplot_file(ofstream *out, ofstream *tecplot, class Tesselation *TesselStruct, class molecule *mol, int N)
 {
   if ((tecplot != NULL) && (TesselStruct != NULL)) {
@@ -308 +308 @@
 
       //*out << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl;
-      if (ProjectedVector.Projection(&AngleReferenceNormalVector) > 0) {
+      if (ProjectedVector.ScalarProduct(&AngleReferenceNormalVector) > 0) {
         angle = 2. * M_PI - angle;
       }
@@ -440 +440 @@
  * \return *TesselStruct is filled with convex boundary and tesselation is stored under \a *filename.
  */
-void FindConvexBorder(ofstream *out, molecule* mol, class LinkedCell *LCList, const char *filename)
+void Find_convex_border(ofstream *out, molecule* mol, class LinkedCell *LCList, const char *filename)
 {
   bool BoundaryFreeFlag = false;
   Boundaries *BoundaryPoints = NULL;
 
-  cout << Verbose(1) << "Begin of FindConvexBorder" << endl;
+  cout << Verbose(1) << "Begin of find_convex_border" << endl;
 
   if (mol->TesselStruct != NULL) // free if allocated
@@ -509 +509 @@
       OutputName.append(TecplotSuffix);
       ofstream *tecplot = new ofstream(OutputName.c_str());
-      WriteTecplotFile(out, tecplot, mol->TesselStruct, mol, 0);
+      write_tecplot_file(out, tecplot, mol->TesselStruct, mol, 0);
       tecplot->close();
       delete(tecplot);
@@ -518 +518 @@
       OutputName.append(Raster3DSuffix);
       ofstream *rasterplot = new ofstream(OutputName.c_str());
-      WriteRaster3dFile(out, rasterplot, mol->TesselStruct, mol);
+      write_raster3d_file(out, rasterplot, mol->TesselStruct, mol);
       rasterplot->close();
       delete(rasterplot);
@@ -556 +556 @@
       OutputName.append(TecplotSuffix);
       ofstream *tecplot = new ofstream(OutputName.c_str());
-      WriteTecplotFile(out, tecplot, mol->TesselStruct, mol, 0);
+      write_tecplot_file(out, tecplot, mol->TesselStruct, mol, 0);
       tecplot->close();
       delete(tecplot);
@@ -564 +564 @@
       OutputName.append(Raster3DSuffix);
       ofstream *rasterplot = new ofstream(OutputName.c_str());
-      WriteRaster3dFile(out, rasterplot, mol->TesselStruct, mol);
+      write_raster3d_file(out, rasterplot, mol->TesselStruct, mol);
       rasterplot->close();
       delete(rasterplot);
@@ -575 +575 @@
     delete[] (BoundaryPoints);
 
-  cout << Verbose(1) << "End of FindConvexBorder" << endl;
+  cout << Verbose(1) << "End of find_convex_border" << endl;
 };
 
@@ -620 +620 @@
   }
 
+  //CalculateConcavityPerBoundaryPoint(out, TesselStruct);
+  StoreTrianglesinFile(out, mol, filename, "-first");
+
   // First step: RemovePointFromTesselatedSurface
-  PointRunner = TesselStruct->PointsOnBoundary.begin();
-  PointAdvance = PointRunner; // we need an advanced point, as the PointRunner might get removed
-  while (PointRunner != TesselStruct->PointsOnBoundary.end()) {
-    PointAdvance++;
-    point = PointRunner->second;
-    *out << Verbose(1) << "INFO: Current point is " << *point << "." << endl;
-    Concavity = true;
-    for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
+  do {
+    Concavity = false;
+    PointRunner = TesselStruct->PointsOnBoundary.begin();
+    PointAdvance = PointRunner; // we need an advanced point, as the PointRunner might get removed
+    while (PointRunner != TesselStruct->PointsOnBoundary.end()) {
+      PointAdvance++;
+      point = PointRunner->second;
+      *out << Verbose(1) << "INFO: Current point is " << *point << "." << endl;
+      for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
+        line = LineRunner->second;
+        *out << Verbose(2) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;
+      }
+      if (!line->CheckConvexityCriterion(out)) {
+        *out << Verbose(1) << "... point " << *point << " cannot be on convex envelope." << endl;
+        // remove the point
+        Concavity = true;
+        TesselStruct->RemovePointFromTesselatedSurface(out, point);
+      }
+      PointRunner = PointAdvance;
+    }
+
+    //CalculateConcavityPerBoundaryPoint(out, TesselStruct);
+    //StoreTrianglesinFile(out, mol, filename, "-second");
+
+    // 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;
-      *out << Verbose(2) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;
-      Concavity = Concavity && (!line->CheckConvexityCriterion(out));
-    }
-    if (Concavity) {
-      *out << Verbose(1) << "... point " << *point << " cannot be on convex envelope." << endl;
-      // remove the point
-      TesselStruct->RemovePointFromTesselatedSurface(out, point);
-    }
-    PointRunner = PointAdvance;
-  }
-
-
-//  // print all lines
-//  LineRunner = TesselStruct->LinesOnBoundary.begin();
-//  LineAdvance = LineRunner;  // we need an advanced line, as the LineRunner might get removed
-//  *out << Verbose(1) << "Printing all boundary lines for debugging." << endl;
-//  while (LineRunner != TesselStruct->LinesOnBoundary.end()) {
-//    LineAdvance++;
-//    line = LineRunner->second;
-//    *out << Verbose(2) << "INFO: Current line is " << *line << "." << endl;
-//    if (LineAdvance != TesselStruct->LinesOnBoundary.end())
-//      *out << Verbose(2) << "INFO: Next line will be " << *(LineAdvance->second) << "." << endl;
-//    LineRunner = LineAdvance;
-//  }
-//
-//  // print all triangles
-//  TriangleRunner = TesselStruct->TrianglesOnBoundary.begin();
-//  TriangleAdvance = TriangleRunner;  // we need an advanced line, as the LineRunner might get removed
-//  *out << Verbose(1) << "Printing all boundary triangles for debugging." << endl;
-//  while (TriangleRunner != TesselStruct->TrianglesOnBoundary.end()) {
-//    TriangleAdvance++;
-//    *out << Verbose(2) << "INFO: Current triangle is " << *(TriangleRunner->second) << "." << endl;
-//    if (TriangleAdvance != TesselStruct->TrianglesOnBoundary.end())
-//      *out << Verbose(2) << "INFO: Next triangle will be " << *(TriangleAdvance->second) << "." << endl;
-//    TriangleRunner = TriangleAdvance;
-//  }
-
-  // second step: PickFarthestofTwoBaselines
+      *out << Verbose(1) << "INFO: Picking farthest baseline for line is " << *line << "." << endl;
+      // take highest of both lines
+      if (TesselStruct->IsConvexRectangle(out, line) == NULL) {
+        TesselStruct->PickFarthestofTwoBaselines(out, line);
+        Concavity = true;
+      }
+      LineRunner = LineAdvance;
+    }
+  } while (Concavity);
+  CalculateConcavityPerBoundaryPoint(out, TesselStruct);
+  StoreTrianglesinFile(out, mol, filename, "-third");
+
+  // third step: IsConvexRectangle
   LineRunner = TesselStruct->LinesOnBoundary.begin();
   LineAdvance = LineRunner;  // we need an advanced line, as the LineRunner might get removed
@@ -673 +672 @@
     LineAdvance++;
     line = LineRunner->second;
-    *out << Verbose(1) << "INFO: Picking farthest baseline for line is " << *line << "." << endl;
-    if (LineAdvance != TesselStruct->LinesOnBoundary.end())
-    // take highest of both lines
-    TesselStruct->PickFarthestofTwoBaselines(out, line);
+    *out << Verbose(1) << "INFO: Current line is " << *line << "." << endl;
+    //if (LineAdvance != TesselStruct->LinesOnBoundary.end())
+      //*out << Verbose(1) << "INFO: Next line will be " << *(LineAdvance->second) << "." << endl;
+    if (!line->CheckConvexityCriterion(out)) {
+      *out << Verbose(1) << "... line " << *line << " is concave, flipping it." << endl;
+
+      // take highest of both lines
+      point = TesselStruct->IsConvexRectangle(out, line);
+      if (point != NULL)
+        TesselStruct->RemovePointFromTesselatedSurface(out, point);
+    }
     LineRunner = LineAdvance;
   }
 
+  CalculateConcavityPerBoundaryPoint(out, TesselStruct);
+  StoreTrianglesinFile(out, mol, filename, "-fourth");
+
+  // end
+  *out << Verbose(0) << "End of ConvexizeNonconvexEnvelope" << endl;
+  return volume;
+};
+
+/** Calculates the concavity for each of the BoundaryPointSet's in a Tesselation.
+ * Sets BoundaryPointSet::value equal to the number of connected lines that are not convex.
+ * \param *out output stream for debugging
+ * \param *TesselStruct pointer to Tesselation structure
+ */
+void CalculateConcavityPerBoundaryPoint(ofstream *out, class Tesselation *TesselStruct)
+{
+  class BoundaryPointSet *point = NULL;
+  class BoundaryLineSet *line = NULL;
   // calculate remaining concavity
-  for (PointRunner = TesselStruct->PointsOnBoundary.begin(); PointRunner != TesselStruct->PointsOnBoundary.end(); PointRunner++) {
+  for (PointMap::iterator PointRunner = TesselStruct->PointsOnBoundary.begin(); PointRunner != TesselStruct->PointsOnBoundary.end(); PointRunner++) {
     point = PointRunner->second;
     *out << Verbose(1) << "INFO: Current point is " << *point << "." << endl;
@@ -692 +715 @@
     }
   }
-
+};
+
+/** Stores triangles to file.
+ * \param *out output stream for debugging
+ * \param *mol molecule with atoms and bonds
+ * \param *filename prefix of filename
+ * \param *extraSuffix intermediate suffix
+ */
+void StoreTrianglesinFile(ofstream *out, molecule *mol, const char *filename, const char *extraSuffix)
+{
   // 4. Store triangles in tecplot file
   if (filename != NULL) {
     if (DoTecplotOutput) {
       string OutputName(filename);
-      OutputName.append("_intermed");
+      OutputName.append(extraSuffix);
       OutputName.append(TecplotSuffix);
       ofstream *tecplot = new ofstream(OutputName.c_str());
-      WriteTecplotFile(out, tecplot, mol->TesselStruct, mol, 0);
+      write_tecplot_file(out, tecplot, mol->TesselStruct, mol, 0);
       tecplot->close();
       delete(tecplot);
@@ -706 +738 @@
     if (DoRaster3DOutput) {
       string OutputName(filename);
-      OutputName.append("_intermed");
+      OutputName.append(extraSuffix);
       OutputName.append(Raster3DSuffix);
       ofstream *rasterplot = new ofstream(OutputName.c_str());
-      WriteRaster3dFile(out, rasterplot, mol->TesselStruct, mol);
+      write_raster3d_file(out, rasterplot, mol->TesselStruct, mol);
       rasterplot->close();
       delete(rasterplot);
     }
   }
-
-  // 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;
-    *out << Verbose(1) << "INFO: Current line is " << *line << "." << endl;
-    if (LineAdvance != TesselStruct->LinesOnBoundary.end())
-      *out << Verbose(1) << "INFO: Next line will be " << *(LineAdvance->second) << "." << endl;
-    if (!line->CheckConvexityCriterion(out)) {
-      *out << Verbose(1) << "... line " << *line << " is concave, flipping it." << endl;
-
-      // take highest of both lines
-      point = TesselStruct->IsConvexRectangle(out, line);
-      if (point != NULL)
-        TesselStruct->RemovePointFromTesselatedSurface(out, point);
-    }
-    LineRunner = LineAdvance;
-  }
-
-  // calculate remaining concavity
-  for (PointRunner = TesselStruct->PointsOnBoundary.begin(); PointRunner != TesselStruct->PointsOnBoundary.end(); PointRunner++) {
-    point = PointRunner->second;
-    *out << Verbose(1) << "INFO: Current point is " << *point << "." << endl;
-    point->value = 0;
-    for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
-      line = LineRunner->second;
-      *out << Verbose(2) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;
-      if (!line->CheckConvexityCriterion(out))
-        point->value += 1;
-    }
-  }
-
-  // 4. Store triangles in tecplot file
-  if (filename != NULL) {
-    if (DoTecplotOutput) {
-      string OutputName(filename);
-      OutputName.append(TecplotSuffix);
-      ofstream *tecplot = new ofstream(OutputName.c_str());
-      WriteTecplotFile(out, tecplot, mol->TesselStruct, mol, 0);
-      tecplot->close();
-      delete(tecplot);
-    }
-    if (DoRaster3DOutput) {
-      string OutputName(filename);
-      OutputName.append(Raster3DSuffix);
-      ofstream *rasterplot = new ofstream(OutputName.c_str());
-      WriteRaster3dFile(out, rasterplot, mol->TesselStruct, mol);
-      rasterplot->close();
-      delete(rasterplot);
-    }
-  }
-
-  // end
-  *out << Verbose(0) << "End of ConvexizeNonconvexEnvelope" << endl;
-  return volume;
 };
 
@@ -804 +779 @@
       G = sqrt(((a + b + c) * (a + b + c) - 2 * (a * a + b * b + c * c)) / 16.); // area of tesselated triangle
       x.MakeNormalVector(runner->second->endpoints[0]->node->node, runner->second->endpoints[1]->node->node, runner->second->endpoints[2]->node->node);
-      x.Scale(runner->second->endpoints[1]->node->node->Projection(&x));
+      x.Scale(runner->second->endpoints[1]->node->node->ScalarProduct(&x));
       h = x.Norm(); // distance of CoG to triangle
       PyramidVolume = (1. / 3.) * G * h; // this formula holds for _all_ pyramids (independent of n-edge base or (not) centered peak)
@@ -841 +816 @@
   class Tesselation *TesselStruct = NULL;
   LinkedCell LCList(mol, 10.);
-  FindConvexBorder(out, mol, &LCList, NULL);
+  Find_convex_border(out, mol, &LCList, NULL);
   double clustervolume;
   if (ClusterVolume == 0)
@@ -982 +957 @@
     if ((*ListRunner)->TesselStruct == NULL) {
       *out << Verbose(1) << "Pre-creating tesselation for molecule " << *ListRunner << "." << endl;
-      FindNonConvexBorder((ofstream *)&cout, (*ListRunner), LCList[i], NULL, 5.);
+      Find_non_convex_border((ofstream *)&cout, (*ListRunner), LCList[i], NULL, 5.);
     }
     i++;
@@ -1107 +1082 @@
  * \para RADIUS radius of the virtual sphere
  */
-void FindNonConvexBorder(ofstream *out, molecule* mol, class LinkedCell *LCList, const char *filename, const double RADIUS)
+void Find_non_convex_border(ofstream *out, molecule* mol, class LinkedCell *LCList, const char *filename, const double RADIUS)
 {
   int N = 0;
@@ -1126 +1101 @@
   LineMap::iterator baseline;
   LineMap::iterator testline;
-  *out << Verbose(0) << "Begin of FindNonConvexBorder\n";
+  *out << Verbose(0) << "Begin of Find_non_convex_border\n";
   bool flag = false;  // marks whether we went once through all baselines without finding any without two triangles
   bool failflag = false;
@@ -1135 +1110 @@
   }
 
-  mol->TesselStruct->FindStartingTriangle(out, RADIUS, LCList);
+  mol->TesselStruct->Find_starting_triangle(out, RADIUS, LCList);
 
   baseline = mol->TesselStruct->LinesOnBoundary.begin();
@@ -1144 +1119 @@
   while ((baseline != mol->TesselStruct->LinesOnBoundary.end()) || (flag)) {
     if (baseline->second->triangles.size() == 1) {
-      failflag = mol->TesselStruct->FindNextSuitableTriangle(out, *(baseline->second), *(((baseline->second->triangles.begin()))->second), RADIUS, N, LCList); //the line is there, so there is a triangle, but only one.
+      failflag = mol->TesselStruct->Find_next_suitable_triangle(out, *(baseline->second), *(((baseline->second->triangles.begin()))->second), RADIUS, N, LCList); //the line is there, so there is a triangle, but only one.
       flag = flag || failflag;
       if (!failflag)
-        cerr << "WARNING: FindNextSuitableTriangle failed." << endl;
+        cerr << "WARNING: Find_next_suitable_triangle failed." << endl;
       // write temporary envelope
       if ((DoSingleStepOutput && (mol->TesselStruct->TrianglesOnBoundaryCount % SingleStepWidth == 0))) { // if we have a new triangle and want to output each new triangle configuration
@@ -1163 +1138 @@
             *out << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
             tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
-            WriteTecplotFile(out, tempstream, mol->TesselStruct, mol, TriangleFilesWritten);
+            write_tecplot_file(out, tempstream, mol->TesselStruct, mol, TriangleFilesWritten);
             tempstream->close();
             tempstream->flush();
@@ -1177 +1152 @@
             *out << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";
             tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc);
-            WriteRaster3dFile(out, tempstream, mol->TesselStruct, mol);
+            write_raster3d_file(out, tempstream, mol->TesselStruct, mol);
     //        // include the current position of the virtual sphere in the temporary raster3d file
     //        // make the circumsphere's center absolute again
@@ -1216 +1191 @@
   }
 
-  // Purges surplus triangles.
-  mol->TesselStruct->RemoveDegeneratedTriangles();
-
   // write final envelope
   if (filename != 0) {
@@ -1226 +1198 @@
       OutputName.append(TecplotSuffix);
       ofstream *tecplot = new ofstream(OutputName.c_str());
-      WriteTecplotFile(out, tecplot, mol->TesselStruct, mol, -1);
+      write_tecplot_file(out, tecplot, mol->TesselStruct, mol, -1);
       tecplot->close();
       delete(tecplot);
@@ -1234 +1206 @@
       OutputName.append(Raster3DSuffix);
       ofstream *raster = new ofstream(OutputName.c_str());
-      WriteRaster3dFile(out, raster, mol->TesselStruct, mol);
+      write_raster3d_file(out, raster, mol->TesselStruct, mol);
       raster->close();
       delete(raster);
@@ -1285 +1257 @@
 //    << "for atom " << a << " (inside)." << endl;
 
+
   if (freeLC)
     delete(LCList);
-  *out << Verbose(0) << "End of FindNonConvexBorder\n";
+  *out << Verbose(0) << "End of Find_non_convex_border\n";
 };