Changeset b8d215 for src/Tesselation
- Timestamp:
- Oct 19, 2014, 5:13:10 PM (10 years ago)
- 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:
- 6d5280
- Parents:
- dfe054
- git-author:
- Frederik Heber <heber@…> (09/25/14 18:26:03)
- git-committer:
- Frederik Heber <heber@…> (10/19/14 17:13:10)
- Location:
- src/Tesselation
- Files:
-
- 8 edited
Legend:
- Unmodified
- Added
- Removed
-
src/Tesselation/BoundaryLineSet.cpp
rdfe054 rb8d215 309 309 ostream & operator <<(ostream &ost, const BoundaryLineSet &a) 310 310 { 311 ost << "[" << a.Nr << "|" << a.endpoints[0]->node->getName() << " at " << a.endpoints[0]->node->getPosition() << "," << a.endpoints[1]->node->getName() << " at " << a.endpoints[1]->node->getPosition() << "]"; 311 ost << "[" << a.Nr << "|" << a.endpoints[0]->node->getName() << "," << a.endpoints[1]->node->getName() << "]"; 312 //ost << "[" << a.Nr << "|" << a.endpoints[0]->node->getName() << " at " << a.endpoints[0]->node->getPosition() << "," << a.endpoints[1]->node->getName() << " at " << a.endpoints[1]->node->getPosition() << "]"; 312 313 return ost; 313 314 } -
src/Tesselation/BoundaryPointSet.cpp
rdfe054 rb8d215 133 133 ostream & operator <<(ostream &ost, const BoundaryPointSet &a) 134 134 { 135 ost << "[" << a.Nr << "|" << a.getName() << " at " << a.getPosition() << "]";135 ost << "[" << a.Nr << "|" << a.getName() << "]"; // " at " << a.getPosition() << "]"; 136 136 return ost; 137 137 } -
src/Tesselation/BoundaryTriangleSet.cpp
rdfe054 rb8d215 94 94 // set endpoints 95 95 int Counter = 0; 96 LOG(4, "DEBUG: New triangle " << Nr << " with end points :");96 LOG(4, "DEBUG: New triangle " << Nr << " with end points ... and lines:"); 97 97 for (PointMap::iterator runner = OrderMap.begin(); runner != OrderMap.end(); runner++) { 98 98 endpoints[Counter] = runner->second; 99 LOG(4, "DEBUG: " << *endpoints[Counter] );99 LOG(4, "DEBUG: " << *endpoints[Counter] << "\t\t" << *lines[Counter]); 100 100 Counter++; 101 101 } … … 246 246 247 247 // 1. get intersection with plane 248 LOG( 3, "DEBUG: Looking for closest point of triangle " << *this << " to " << x << ".");249 LOG( 3, "DEBUG: endpoints are " << endpoints[0]->node->getPosition() << ","248 LOG(4, "DEBUG: Looking for closest point of triangle " << *this << " to " << x << "."); 249 LOG(5, "DEBUG: endpoints are " << endpoints[0]->node->getPosition() << "," 250 250 << endpoints[1]->node->getPosition() << ", and " << endpoints[2]->node->getPosition() << "."); 251 251 try { … … 256 256 } 257 257 Vector InPlane(ClosestPoint); // points from plane intersection to straight-down point 258 LOG( 4, "DEBUG: Closest point on triangle plane is " << ClosestPoint << ".");258 LOG(5, "DEBUG: Closest point on triangle plane is " << ClosestPoint << "."); 259 259 260 260 // 2. Calculate in plane part of line (x, intersection) … … 269 269 CrossPoint[i] = l.getClosestPoint(InPlane); 270 270 // NOTE: direction of line is normalized, hence s must not necessarily be in [0,1] for the baseline 271 LOG( 4, "DEBUG: Closest point on line from " << (endpoints[(i+1)%3]->node->getPosition())271 LOG(5, "DEBUG: Closest point on line from " << (endpoints[(i+1)%3]->node->getPosition()) 272 272 << " to " << (endpoints[i%3]->node->getPosition()) << " is " << CrossPoint[i] << "."); 273 273 CrossPoint[i] -= (endpoints[(i+1)%3]->node->getPosition()); // cross point was returned as absolute vector 274 274 const double s = CrossPoint[i].ScalarProduct(Direction)/Direction.NormSquared(); 275 LOG( 5, "DEBUG: Factor s is " << s << ".");275 LOG(6, "DEBUG: Factor s is " << s << "."); 276 276 if ((s >= -MYEPSILON) && ((s-1.) <= MYEPSILON)) { 277 277 CrossPoint[i] += (endpoints[(i+1)%3]->node->getPosition()); // make cross point absolute again 278 LOG( 5, "DEBUG: Crosspoint is " << CrossPoint[i] << ", intersecting BoundaryLine between "278 LOG(6, "DEBUG: Crosspoint is " << CrossPoint[i] << ", intersecting BoundaryLine between " 279 279 << endpoints[i % 3]->node->getPosition() << " and " 280 280 << endpoints[(i + 1) % 3]->node->getPosition() << "."); … … 285 285 else 286 286 CrossPoint[i] = (endpoints[i%3]->node->getPosition()); 287 LOG( 5, "DEBUG: Crosspoint is " << CrossPoint[i] << ", intersecting outside of BoundaryLine between "287 LOG(6, "DEBUG: Crosspoint is " << CrossPoint[i] << ", intersecting outside of BoundaryLine between " 288 288 << endpoints[i % 3]->node->getPosition() << " and " 289 289 << endpoints[(i + 1) % 3]->node->getPosition() << "."); … … 312 312 ShortestDistance = InPlane.DistanceSquared(x); 313 313 } 314 LOG( 3, "DEBUG: Closest Point is " << ClosestPoint << " with shortest squared distance is " << ShortestDistance << ".");314 LOG(4, "DEBUG: Closest Point is " << ClosestPoint << " with shortest squared distance is " << ShortestDistance << "."); 315 315 316 316 return ShortestDistance; … … 367 367 { 368 368 //Info FunctionInfo(__func__); 369 LOG(1, "INFO: Checking " << Points[0] << "," << Points[1] << "," << Points[2] << " against " << endpoints[0] << "," << endpoints[1] << "," << endpoints[2] << "."); 370 return (((endpoints[0] == Points[0]) || (endpoints[0] == Points[1]) || (endpoints[0] == Points[2])) && ((endpoints[1] == Points[0]) || (endpoints[1] == Points[1]) || (endpoints[1] == Points[2])) && ((endpoints[2] == Points[0]) || (endpoints[2] == Points[1]) || (endpoints[2] == Points[2]) 369 LOG(5, "DEBUG: Checking " << *Points[0] << "," << *Points[1] << "," << *Points[2] 370 << " against " << *this); //*endpoints[0] << "," << *endpoints[1] << "," << *endpoints[2] << "."); 371 return (((endpoints[0] == Points[0]) || (endpoints[0] == Points[1]) || (endpoints[0] == Points[2])) 372 && ((endpoints[1] == Points[0]) || (endpoints[1] == Points[1]) || (endpoints[1] == Points[2])) 373 && ((endpoints[2] == Points[0]) || (endpoints[2] == Points[1]) || (endpoints[2] == Points[2]) 371 374 372 375 )); -
src/Tesselation/CandidateForTesselation.cpp
rdfe054 rb8d215 108 108 109 109 if (!pointlist.empty()) 110 LOG(3, "DEBUG: Checking whether sphere contains candidate list and baseline " << *BaseLine->endpoints[0] << "<->" << *BaseLine->endpoints[1] << " only ...");110 LOG(3, "DEBUG: Checking validity whether sphere contains candidate list and baseline " << *BaseLine->endpoints[0] << "<->" << *BaseLine->endpoints[1] << " only ..."); 111 111 else 112 LOG(3, "DEBUG: Checking whether sphere with no candidates contains baseline " << *BaseLine->endpoints[0] << "<->" << *BaseLine->endpoints[1] << " only ...");112 LOG(3, "DEBUG: Checking validity whether sphere with no candidates contains baseline " << *BaseLine->endpoints[0] << "<->" << *BaseLine->endpoints[1] << " only ..."); 113 113 // check baseline for OptCenter and OtherOptCenter being on sphere's surface 114 114 for (list<const Vector *>::const_iterator VRunner = VectorList.begin(); VRunner != VectorList.end(); ++VRunner) { … … 131 131 return false; 132 132 } else { 133 LOG( 3, "DEBUG: Candidate " << *Walker << " is inside by " << distance << ".");133 LOG(4, "DEBUG: Candidate " << *Walker << " is inside by " << distance << "."); 134 134 } 135 135 } 136 136 } 137 137 138 LOG( 2, "DEBUG: Checkingwhether sphere contains no others points ...");138 LOG(3, "DEBUG: Checking validity whether sphere contains no others points ..."); 139 139 bool flag = true; 140 140 for (list<const Vector *>::const_iterator VRunner = VectorList.begin(); VRunner != VectorList.end(); ++VRunner) { … … 143 143 144 144 { 145 LOG( 3, "DEBUG: The following atoms are inside sphere at " << (*VRunner) << ":");145 LOG(4, "DEBUG: The following atoms are inside sphere at " << (*VRunner) << ":"); 146 146 for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner) 147 LOG( 3, "DEBUG: " << *(*Runner) << " with distance " << (*Runner)->distance(*(*VRunner)) << ".");147 LOG(4, "DEBUG: " << *(*Runner) << " with distance " << (*Runner)->distance(*(*VRunner)) << "."); 148 148 } 149 149 150 150 // remove baseline's endpoints and candidates 151 151 for (int i = 0; i < 2; i++) { 152 LOG( 3, "DEBUG: removing baseline tesselpoint " << *BaseLine->endpoints[i]->node << ".");152 LOG(5, "DEBUG: removing baseline tesselpoint " << *BaseLine->endpoints[i]->node << "."); 153 153 ListofPoints->remove(BaseLine->endpoints[i]->node); 154 154 } 155 155 for (TesselPointList::const_iterator Runner = pointlist.begin(); Runner != pointlist.end(); ++Runner) { 156 LOG( 3, "DEBUG: removing candidate tesselpoint " << *(*Runner) << ".");156 LOG(5, "DEBUG: removing candidate tesselpoint " << *(*Runner) << "."); 157 157 ListofPoints->remove(*Runner); 158 158 } -
src/Tesselation/boundary.cpp
rdfe054 rb8d215 535 535 } 536 536 537 LOG(1, "INFO: Making tesselated surface with " << TesselStruct->TrianglesOnBoundaryCount 538 << " convex ..."); 539 537 540 // First step: RemovePointFromTesselatedSurface 538 541 do { … … 547 550 PointAdvance++; 548 551 point = PointRunner->second; 549 LOG( 1, "INFO: Current point is " << *point << ".");552 LOG(2, "INFO: Current point is " << *point << "."); 550 553 for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) { 551 554 line = LineRunner->second; 552 LOG( 1, "INFO: Current line of point " << *point << " is " << *line << ".");555 LOG(3, "INFO: Current line of point " << *point << " is " << *line << "."); 553 556 if (!line->CheckConvexityCriterion()) { 554 557 // remove the point if needed … … 766 769 // status = CheckListOfBaselines(TesselStruct); 767 770 768 cout << "before correction" << endl;771 // cout << "before correction" << endl; 769 772 770 773 // store before correction … … 779 782 // write final envelope 780 783 CalculateConcavityPerBoundaryPoint(TesselStruct); 781 cout << "after correction" << endl;784 // cout << "after correction" << endl; 782 785 StoreTrianglesinFile(mol, TesselStruct, filename, ""); 783 786 -
src/Tesselation/tesselation.cpp
rdfe054 rb8d215 69 69 class molecule; 70 70 71 const char *TecplotSuffix =".dat";72 const char *Raster3DSuffix =".r3d";73 const char *VRMLSUffix =".wrl";74 75 const double ParallelEpsilon =1e-3;71 const char *TecplotSuffix = ".dat"; 72 const char *Raster3DSuffix = ".r3d"; 73 const char *VRMLSUffix = ".wrl"; 74 75 const double ParallelEpsilon = 1e-3; 76 76 const double Tesselation::HULLEPSILON = 1e-9; 77 77 … … 79 79 */ 80 80 Tesselation::Tesselation() : 81 PointsOnBoundaryCount(0), 82 LinesOnBoundaryCount(0), 83 TrianglesOnBoundaryCount(0), 84 LastTriangle(NULL), 85 TriangleFilesWritten(0), 86 InternalPointer(PointsOnBoundary.begin()) 81 PointsOnBoundaryCount(0), LinesOnBoundaryCount(0), TrianglesOnBoundaryCount(0), LastTriangle(NULL), TriangleFilesWritten(0), InternalPointer(PointsOnBoundary.begin()) 87 82 { 88 83 //Info FunctionInfo(__func__); … … 116 111 { 117 112 // create linkedcell 118 LinkedCell_deprecated *LinkedList = new LinkedCell_deprecated(cloud, 2.*SPHERERADIUS);119 120 121 122 123 124 for (size_t i =0;i<3;++i, cloud.GoToNext())113 LinkedCell_deprecated *LinkedList = new LinkedCell_deprecated(cloud, 2. * SPHERERADIUS); 114 115 // check for at least three points 116 { 117 bool ThreePointsFound = true; 118 cloud.GoToFirst(); 119 for (size_t i = 0; i < 3; ++i, cloud.GoToNext()) 125 120 ThreePointsFound &= (!cloud.IsEnd()); 126 121 cloud.GoToFirst(); … … 129 124 return; 130 125 } 131 132 133 126 } 127 128 // find a starting triangle 134 129 FindStartingTriangle(SPHERERADIUS, LinkedList); 135 130 … … 145 140 //the line is there, so there is a triangle, but only one. 146 141 const bool TesselationFailFlag = FindNextSuitableTriangle(*baseline, *T, SPHERERADIUS, LinkedList); 147 ASSERT( TesselationFailFlag, 148 "Tesselation::operator() - no suitable candidate triangle found."); 142 ASSERT(TesselationFailFlag, "Tesselation::operator() - no suitable candidate triangle found."); 149 143 } 150 144 } 151 145 152 146 // 2b. search for smallest ShortestAngle among all candidates 153 double ShortestAngle = 4. *M_PI;147 double ShortestAngle = 4. * M_PI; 154 148 for (CandidateMap::iterator Runner = OpenLines.begin(); Runner != OpenLines.end(); Runner++) { 155 149 if (Runner->second->ShortestAngle < ShortestAngle) { … … 158 152 } 159 153 } 160 if ((ShortestAngle == 4. *M_PI) || (baseline->pointlist.empty()))154 if ((ShortestAngle == 4. * M_PI) || (baseline->pointlist.empty())) 161 155 OneLoopWithoutSuccessFlag = false; 162 156 else { … … 180 174 181 175 // 6a. Every triangle forms a pyramid with the center of gravity as its peak, sum up the volumes 182 for (TriangleMap::const_iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) 183 { // go through every triangle, calculate volume of its pyramid with CoG as peak 184 x = runner->second->getEndpoint(0) - runner->second->getEndpoint(1); 185 const double G = runner->second->getArea(); 186 x = runner->second->getPlane().getNormal(); 187 x.Scale(runner->second->getEndpoint(1).ScalarProduct(x)); 188 const double h = x.Norm(); // distance of CoG to triangle 189 const double PyramidVolume = (1. / 3.) * G * h; // this formula holds for _all_ pyramids (independent of n-edge base or (not) centered peak) 190 LOG(1, "INFO: Area of triangle is " << setprecision(10) << G << " " 191 << (IsAngstroem ? "angstrom" : "atomiclength") << "^2, height is " 192 << h << " and the volume is " << PyramidVolume << " " 193 << (IsAngstroem ? "angstrom" : "atomiclength") << "^3."); 194 volume += PyramidVolume; 195 } 196 LOG(0, "RESULT: The summed volume is " << setprecision(6) 197 << volume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3."); 176 for (TriangleMap::const_iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) { // go through every triangle, calculate volume of its pyramid with CoG as peak 177 x = runner->second->getEndpoint(0) - runner->second->getEndpoint(1); 178 const double G = runner->second->getArea(); 179 x = runner->second->getPlane().getNormal(); 180 x.Scale(runner->second->getEndpoint(1).ScalarProduct(x)); 181 const double h = x.Norm(); // distance of CoG to triangle 182 const double PyramidVolume = (1. / 3.) * G * h; // this formula holds for _all_ pyramids (independent of n-edge base or (not) centered peak) 183 LOG(1, "INFO: Area of triangle is " << setprecision(10) << G << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^2, height is " << h << " and the volume is " << PyramidVolume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3."); 184 volume += PyramidVolume; 185 } 186 LOG(0, "RESULT: The summed volume is " << setprecision(6) << volume << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3."); 198 187 199 188 return volume; … … 212 201 213 202 // 6a. Every triangle forms a pyramid with the center of gravity as its peak, sum up the volumes 214 for (TriangleMap::const_iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) 215 { // go through every triangle, calculate volume of its pyramid with CoG as peak 216 const double area = runner->second->getArea(); 217 LOG(1, "INFO: Area of triangle is " << setprecision(10) << area << " " 218 << (IsAngstroem ? "angstrom" : "atomiclength") << "^2."); 219 surfacearea += area; 220 } 221 LOG(0, "RESULT: The summed surface area is " << setprecision(6) 222 << surfacearea << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3."); 203 for (TriangleMap::const_iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) { // go through every triangle, calculate volume of its pyramid with CoG as peak 204 const double area = runner->second->getArea(); 205 LOG(1, "INFO: Area of triangle is " << setprecision(10) << area << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^2."); 206 surfacearea += area; 207 } 208 LOG(0, "RESULT: The summed surface area is " << setprecision(6) << surfacearea << " " << (IsAngstroem ? "angstrom" : "atomiclength") << "^3."); 223 209 224 210 return surfacearea; 225 211 } 226 227 212 228 213 /** Gueses first starting triangle of the convex envelope. … … 256 241 tmp = B->second->node->DistanceSquared(C->second->node->getPosition()); 257 242 distance += tmp * tmp; 258 DistanceMMap.insert(DistanceMultiMapPair(distance, pair<PointMap::iterator, PointMap::iterator> 243 DistanceMMap.insert(DistanceMultiMapPair(distance, pair<PointMap::iterator, PointMap::iterator>(B, C))); 259 244 } 260 245 } … … 276 261 // 2. next, we have to check whether all points reside on only one side of the triangle 277 262 // 3. construct plane vector 278 PlaneVector = Plane(A->second->node->getPosition(), 279 baseline->second.first->second->node->getPosition(), 280 baseline->second.second->second->node->getPosition()).getNormal(); 263 PlaneVector = Plane(A->second->node->getPosition(), baseline->second.first->second->node->getPosition(), baseline->second.second->second->node->getPosition()).getNormal(); 281 264 LOG(2, "Plane vector of candidate triangle is " << PlaneVector); 282 265 // 4. loop over all points … … 393 376 // prepare some auxiliary vectors 394 377 Vector BaseLineCenter, BaseLine; 395 BaseLineCenter = 0.5 * ((baseline->second->endpoints[0]->node->getPosition()) + 396 (baseline->second->endpoints[1]->node->getPosition())); 378 BaseLineCenter = 0.5 * ((baseline->second->endpoints[0]->node->getPosition()) + (baseline->second->endpoints[1]->node->getPosition())); 397 379 BaseLine = (baseline->second->endpoints[0]->node->getPosition()) - (baseline->second->endpoints[1]->node->getPosition()); 398 380 … … 412 394 // vector in propagation direction (out of triangle) 413 395 // project center vector onto triangle plane (points from intersection plane-NormalVector to plane-CenterVector intersection) 414 PropagationVector = Plane(BaseLine, NormalVector, 0).getNormal();396 PropagationVector = Plane(BaseLine, NormalVector, 0).getNormal(); 415 397 TempVector = CenterVector - (baseline->second->endpoints[0]->node->getPosition()); // TempVector is vector on triangle plane pointing from one baseline egde towards center! 416 398 //LOG(0, "Projection of propagation onto temp: " << PropagationVector.Projection(&TempVector) << "."); … … 465 447 // in case NOT both were found, create virtually this triangle, get its normal vector, calculate angle 466 448 flag = true; 467 VirtualNormalVector = Plane((baseline->second->endpoints[0]->node->getPosition()), 468 (baseline->second->endpoints[1]->node->getPosition()), 469 (target->second->node->getPosition())).getNormal(); 470 TempVector = (1./3.) * ((baseline->second->endpoints[0]->node->getPosition()) + 471 (baseline->second->endpoints[1]->node->getPosition()) + 472 (target->second->node->getPosition())); 449 VirtualNormalVector = Plane((baseline->second->endpoints[0]->node->getPosition()), (baseline->second->endpoints[1]->node->getPosition()), (target->second->node->getPosition())).getNormal(); 450 TempVector = (1. / 3.) * ((baseline->second->endpoints[0]->node->getPosition()) + (baseline->second->endpoints[1]->node->getPosition()) + (target->second->node->getPosition())); 473 451 TempVector -= (*Center); 474 452 // make it always point outward … … 482 460 winner = target; 483 461 LOG(5, "DEBUG: New winner " << *winner->second->node << " due to smaller angle between normal vectors."); 484 } else if (fabs(SmallestAngle - TempAngle) < MYEPSILON) { // check the angle to propagation, both possible targets are in one plane! (their normals have same angle) 485 // hence, check the angles to some normal direction from our base line but in this common plane of both targets... 486 helper = (target->second->node->getPosition()) - BaseLineCenter; 487 helper.ProjectOntoPlane(BaseLine); 488 // ...the one with the smaller angle is the better candidate 489 TempVector = (target->second->node->getPosition()) - BaseLineCenter; 490 TempVector.ProjectOntoPlane(VirtualNormalVector); 491 TempAngle = TempVector.Angle(helper); 492 TempVector = (winner->second->node->getPosition()) - BaseLineCenter; 493 TempVector.ProjectOntoPlane(VirtualNormalVector); 494 if (TempAngle < TempVector.Angle(helper)) { 495 TempAngle = NormalVector.Angle(VirtualNormalVector); 496 SmallestAngle = TempAngle; 497 winner = target; 498 LOG(5, "DEBUG: New winner " << *winner->second->node << " due to smaller angle " << TempAngle << " to propagation direction."); 462 } else 463 if (fabs(SmallestAngle - TempAngle) < MYEPSILON) { // check the angle to propagation, both possible targets are in one plane! (their normals have same angle) 464 // hence, check the angles to some normal direction from our base line but in this common plane of both targets... 465 helper = (target->second->node->getPosition()) - BaseLineCenter; 466 helper.ProjectOntoPlane(BaseLine); 467 // ...the one with the smaller angle is the better candidate 468 TempVector = (target->second->node->getPosition()) - BaseLineCenter; 469 TempVector.ProjectOntoPlane(VirtualNormalVector); 470 TempAngle = TempVector.Angle(helper); 471 TempVector = (winner->second->node->getPosition()) - BaseLineCenter; 472 TempVector.ProjectOntoPlane(VirtualNormalVector); 473 if (TempAngle < TempVector.Angle(helper)) { 474 TempAngle = NormalVector.Angle(VirtualNormalVector); 475 SmallestAngle = TempAngle; 476 winner = target; 477 LOG(5, "DEBUG: New winner " << *winner->second->node << " due to smaller angle " << TempAngle << " to propagation direction."); 478 } else 479 LOG(5, "DEBUG: Keeping old winner " << *winner->second->node << " due to smaller angle to propagation direction."); 499 480 } else 500 LOG(5, "DEBUG: Keeping old winner " << *winner->second->node << " due to smaller angle to propagation direction."); 501 } else 502 LOG(5, "DEBUG: Keeping old winner " << *winner->second->node << " due to smaller angle between normal vectors."); 481 LOG(5, "DEBUG: Keeping old winner " << *winner->second->node << " due to smaller angle between normal vectors."); 503 482 } 504 483 } // end of loop over all boundary points … … 567 546 568 547 cloud.GoToFirst(); 569 PointCloudAdaptor< 548 PointCloudAdaptor<Tesselation, MapValueIterator<Tesselation::iterator> > newcloud(this, cloud.GetName()); 570 549 BoundaryPoints = new LinkedCell_deprecated(newcloud, 5.); 571 550 while (!cloud.IsEnd()) { // we only have to go once through all points, as boundary can become only bigger … … 752 731 if (FindLine->second->triangles.size() == 1) { 753 732 CandidateMap::iterator Finder = OpenLines.find(FindLine->second); 754 ASSERT( Finder != OpenLines.end(), 755 "Tesselation::AddTesselationLine() - "+toString(*FindLine->second) 756 +" is not a new line and not in OpenLines."); 733 ASSERT(Finder != OpenLines.end(), "Tesselation::AddTesselationLine() - " + toString(*FindLine->second) + " is not a new line and not in OpenLines."); 757 734 if (Finder->second != NULL) { 758 735 if (!Finder->second->pointlist.empty()) … … 765 742 // get open line 766 743 for (TesselPointList::const_iterator CandidateChecker = Finder->second->pointlist.begin(); CandidateChecker != Finder->second->pointlist.end(); ++CandidateChecker) { 767 if ((*(CandidateChecker) == candidate->node) && (OptCenter == NULL || OptCenter->DistanceSquared(Finder->second->OptCenter) < MYEPSILON 744 if ((*(CandidateChecker) == candidate->node) && (OptCenter == NULL || OptCenter->DistanceSquared(Finder->second->OptCenter) < MYEPSILON)) { // stop searching if candidate matches 768 745 LOG(4, "ACCEPT: Candidate " << *(*CandidateChecker) << " has the right center " << Finder->second->OptCenter << "."); 769 746 insertNewLine = false; … … 814 791 // also add to open lines 815 792 CandidateForTesselation *CFT = new CandidateForTesselation(BLS[n]); 816 OpenLines.insert(pair<BoundaryLineSet *, CandidateForTesselation *> 793 OpenLines.insert(pair<BoundaryLineSet *, CandidateForTesselation *>(BLS[n], CFT)); 817 794 } 818 795 ; … … 905 882 for (TriangleMap::iterator TriangleRunner = triangle->lines[i]->triangles.begin(); TriangleRunner != triangle->lines[i]->triangles.end(); TriangleRunner++) 906 883 output << "\t[" << (TriangleRunner->second)->Nr << "|" << *((TriangleRunner->second)->endpoints[0]) << ", " << *((TriangleRunner->second)->endpoints[1]) << ", " << *((TriangleRunner->second)->endpoints[2]) << "] \t"; 907 OpenLines.insert(pair<BoundaryLineSet *, CandidateForTesselation *> 884 OpenLines.insert(pair<BoundaryLineSet *, CandidateForTesselation *>(triangle->lines[i], NULL)); 908 885 } 909 886 LOG(3, "DEBUG: " << output.str()); … … 964 941 LOG(4, "DEBUG: " << *line->endpoints[i] << " has no more lines it's attached to, erasing."); 965 942 RemoveTesselationPoint(line->endpoints[i]); 966 } else if (DoLog(0)) { 967 std::stringstream output; 968 output << "DEBUG: " << *line->endpoints[i] << " has still lines it's attached to: "; 969 for (LineMap::iterator LineRunner = line->endpoints[i]->lines.begin(); LineRunner != line->endpoints[i]->lines.end(); LineRunner++) 970 output << "[" << *(LineRunner->second) << "] \t"; 971 LOG(4, output.str()); 972 } 943 } else 944 if (DoLog(0)) { 945 std::stringstream output; 946 output << "DEBUG: " << *line->endpoints[i] << " has still lines it's attached to: "; 947 for (LineMap::iterator LineRunner = line->endpoints[i]->lines.begin(); LineRunner != line->endpoints[i]->lines.end(); LineRunner++) 948 output << "[" << *(LineRunner->second) << "] \t"; 949 LOG(4, output.str()); 950 } 973 951 line->endpoints[i] = NULL; // free'd or not: disconnect 974 952 } else … … 1012 990 //Info FunctionInfo(__func__); 1013 991 1014 LOG(3, "DEBUG: Checking whether sphere contains no others points ...");992 LOG(3, "DEBUG: Checking degeneracy by whether sphere contains no others points ..."); 1015 993 bool flag = true; 1016 994 … … 1019 997 TesselPointList *ListofPoints = LC->GetPointsInsideSphere(RADIUS, &CandidateLine.OtherOptCenter); 1020 998 1021 LOG(3, "DEBUG: The following atoms are inside sphere at " << CandidateLine.OtherOptCenter << ":"); 999 LOG(3, "DEBUG: CheckDegeneracy: There are " << ListofPoints->size() << " points inside the sphere."); 1000 LOG(4, "DEBUG: The following atoms are inside sphere at " << CandidateLine.OtherOptCenter << ":"); 1022 1001 for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner) 1023 LOG( 3, "DEBUG: " << *(*Runner) << " with distance " << (*Runner)->distance(CandidateLine.OtherOptCenter) << ".");1002 LOG(4, "DEBUG: " << *(*Runner) << " with distance " << (*Runner)->distance(CandidateLine.OtherOptCenter) << "."); 1024 1003 1025 1004 // remove triangles's endpoints … … 1035 1014 LOG(3, "DEBUG: CheckDegeneracy: There are still " << ListofPoints->size() << " points inside the sphere."); 1036 1015 flag = false; 1037 LOG( 3, "DEBUG: External atoms inside of sphere at " << CandidateLine.OtherOptCenter << ":");1016 LOG(4, "DEBUG: External atoms inside of sphere at " << CandidateLine.OtherOptCenter << ":"); 1038 1017 for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner) 1039 LOG( 3, "DEBUG: " << *(*Runner) << " with distance " << (*Runner)->distance(CandidateLine.OtherOptCenter) << ".");1018 LOG(4, "DEBUG: " << *(*Runner) << " with distance " << (*Runner)->distance(CandidateLine.OtherOptCenter) << "."); 1040 1019 } 1041 1020 delete ListofPoints; … … 1182 1161 1183 1162 // 1. searching topmost point with respect to each axis 1163 LOG(2, "INFO: Searching for topmost pointer from each axis."); 1184 1164 for (int i = 0; i < NDIM; i++) { // each axis 1185 1165 LC->n[i] = LC->N[i] - 1; // current axis is topmost cell … … 1203 1183 } 1204 1184 1205 if (DoLog( 1)) {1185 if (DoLog(3)) { 1206 1186 std::stringstream output; 1207 1187 output << "Found maximum coordinates: "; … … 1217 1197 BaseLine = new BoundaryLineSet(); 1218 1198 BaseLine->endpoints[0] = new BoundaryPointSet(MaxPoint[k]); 1219 LOG(2, " DEBUG: Coordinates of start node at " << *BaseLine->endpoints[0]->node << ".");1199 LOG(2, "INFO: Coordinates of start node at " << *BaseLine->endpoints[0]->node << "."); 1220 1200 1221 1201 double ShortestAngle; … … 1230 1210 } 1231 1211 BaseLine->endpoints[1] = new BoundaryPointSet(Temporary); 1232 LOG( 1, "INFO: Second node is at " << *Temporary << ".");1212 LOG(2, "INFO: Second node is at " << *Temporary << "."); 1233 1213 1234 1214 // construct center of circle 1235 1215 CircleCenter = 0.5 * ((BaseLine->endpoints[0]->node->getPosition()) + (BaseLine->endpoints[1]->node->getPosition())); 1236 LOG( 1, "INFO: CircleCenter is at " << CircleCenter << ".");1216 LOG(4, "DEBUG: CircleCenter is at " << CircleCenter << "."); 1237 1217 1238 1218 // construct normal vector of circle 1239 1219 CirclePlaneNormal = (BaseLine->endpoints[0]->node->getPosition()) - (BaseLine->endpoints[1]->node->getPosition()); 1240 LOG( 1, "INFO: CirclePlaneNormal is at " << CirclePlaneNormal << ".");1220 LOG(4, "DEBUG: CirclePlaneNormal is at " << CirclePlaneNormal << "."); 1241 1221 1242 1222 double radius = CirclePlaneNormal.NormSquared(); … … 1245 1225 NormalVector.ProjectOntoPlane(CirclePlaneNormal); 1246 1226 NormalVector.Normalize(); 1247 LOG( 1, "INFO: NormalVector is at " << NormalVector << ".");1227 LOG(4, "DEBUG: NormalVector is at " << NormalVector << "."); 1248 1228 ShortestAngle = 2. * M_PI; // This will indicate the quadrant. 1249 1229 1250 SphereCenter = (CircleRadius * NormalVector) + 1230 SphereCenter = (CircleRadius * NormalVector) + CircleCenter; 1251 1231 // Now, NormalVector and SphereCenter are two orthonormalized vectors in the plane defined by CirclePlaneNormal (not normalized) 1252 1232 1253 1233 // look in one direction of baseline for initial candidate 1254 1234 try { 1255 SearchDirection = Plane(CirclePlaneNormal, NormalVector,0).getNormal(); // whether we look "left" first or "right" first is not important ... 1256 } catch(LinearAlgebraException) { 1257 ELOG(1, "Vectors are linear dependent: " 1258 << CirclePlaneNormal << ", " << NormalVector << "."); 1235 SearchDirection = Plane(CirclePlaneNormal, NormalVector, 0).getNormal(); // whether we look "left" first or "right" first is not important ... 1236 } catch (LinearAlgebraException) { 1237 ELOG(1, "Vectors are linear dependent: " << CirclePlaneNormal << ", " << NormalVector << "."); 1259 1238 delete BaseLine; 1260 1239 continue; … … 1262 1241 1263 1242 // adding point 1 and point 2 and add the line between them 1264 LOG( 2, "DEBUG: Found second point is at " << *BaseLine->endpoints[1]->node << ".");1243 LOG(3, "DEBUG: Found second point is at " << *BaseLine->endpoints[1]->node << "."); 1265 1244 1266 1245 //LOG(1, "INFO: OldSphereCenter is at " << helper << "."); … … 1271 1250 for (TesselPointList::iterator it = OptCandidates.pointlist.begin(); it != OptCandidates.pointlist.end(); it++) 1272 1251 output << *(*it); 1273 LOG( 2, "DEBUG: List of third Points is: " << output.str());1252 LOG(3, "DEBUG: List of third Points is: " << output.str()); 1274 1253 } 1275 1254 if (!OptCandidates.pointlist.empty()) { … … 1433 1412 // return result; 1434 1413 //}; 1435 1436 1414 /** This function finds a triangle to a line, adjacent to an existing one. 1437 1415 * @param out output stream for debugging … … 1465 1443 1466 1444 // construct center of circle 1467 CircleCenter = 0.5 * ((CandidateLine.BaseLine->endpoints[0]->node->getPosition()) + 1468 (CandidateLine.BaseLine->endpoints[1]->node->getPosition())); 1445 CircleCenter = 0.5 * ((CandidateLine.BaseLine->endpoints[0]->node->getPosition()) + (CandidateLine.BaseLine->endpoints[1]->node->getPosition())); 1469 1446 1470 1447 // construct normal vector of circle 1471 CirclePlaneNormal = (CandidateLine.BaseLine->endpoints[0]->node->getPosition()) - 1472 (CandidateLine.BaseLine->endpoints[1]->node->getPosition()); 1448 CirclePlaneNormal = (CandidateLine.BaseLine->endpoints[0]->node->getPosition()) - (CandidateLine.BaseLine->endpoints[1]->node->getPosition()); 1473 1449 1474 1450 // calculate squared radius of circle … … 1480 1456 CircleRadius = RADIUS * RADIUS - radius / 4.; 1481 1457 CirclePlaneNormal.Normalize(); 1482 LOG( 3, "DEBUG: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << ".");1483 1484 LOG( 3, "DEBUG: OldSphereCenter is at " << T.SphereCenter << ".");1458 LOG(4, "DEBUG: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "."); 1459 1460 LOG(4, "DEBUG: OldSphereCenter is at " << T.SphereCenter << "."); 1485 1461 1486 1462 // construct SearchDirection and an "outward pointer" 1487 SearchDirection = Plane(RelativeSphereCenter, CirclePlaneNormal, 0).getNormal();1463 SearchDirection = Plane(RelativeSphereCenter, CirclePlaneNormal, 0).getNormal(); 1488 1464 helper = CircleCenter - (ThirdPoint->node->getPosition()); 1489 if (helper.ScalarProduct(SearchDirection) < -HULLEPSILON) // ohoh, SearchDirection points inwards!1465 if (helper.ScalarProduct(SearchDirection) < -HULLEPSILON) // ohoh, SearchDirection points inwards! 1490 1466 SearchDirection.Scale(-1.); 1491 LOG( 3, "DEBUG: SearchDirection is " << SearchDirection << ".");1467 LOG(4, "DEBUG: SearchDirection is " << SearchDirection << "."); 1492 1468 if (fabs(RelativeSphereCenter.ScalarProduct(SearchDirection)) > HULLEPSILON) { 1493 1469 // rotated the wrong way! … … 1499 1475 1500 1476 } else { 1501 LOG( 3, "DEBUG: Circumcircle for base line " << *CandidateLine.BaseLine << " and base triangle " << T << " is too big!");1477 LOG(4, "DEBUG: Circumcircle for base line " << *CandidateLine.BaseLine << " and base triangle " << T << " is too big!"); 1502 1478 } 1503 1479 … … 1532 1508 baseline = Runner->second; 1533 1509 if (baseline->pointlist.empty()) { 1534 ASSERT((baseline->BaseLine->triangles.size() == 1), "Open line without exactly one attached triangle");1510 ASSERT((baseline->BaseLine->triangles.size() == 1), "Open line without exactly one attached triangle"); 1535 1511 T = (((baseline->BaseLine->triangles.begin()))->second); 1536 1512 LOG(4, "DEBUG: Finding best candidate for open line " << *baseline->BaseLine << " of triangle " << *T); … … 1565 1541 TesselPointList *connectedClosestPoints = GetCircleOfSetOfPoints(&SetOfNeighbours, TurningPoint, CandidateLine.BaseLine->endpoints[1]->node->getPosition()); 1566 1542 1567 {1543 if (DoLog(3)) { 1568 1544 std::stringstream output; 1569 1545 for (TesselPointList::iterator TesselRunner = connectedClosestPoints->begin(); TesselRunner != connectedClosestPoints->end(); ++TesselRunner) 1570 1546 output << **TesselRunner; 1571 LOG(3, "DEBUG: List of Candidates for Turning Point " << *TurningPoint << ":" );1547 LOG(3, "DEBUG: List of Candidates for Turning Point " << *TurningPoint << ":" << output.str()); 1572 1548 } 1573 1549 … … 1610 1586 } 1611 1587 delete (connectedClosestPoints); 1612 }; 1588 } 1589 ; 1613 1590 1614 1591 /** for polygons (multiple candidates for a baseline) sets internal edges to the correct next candidate. … … 1630 1607 else { 1631 1608 LOG(4, "DEBUG: line " << *(FindLine->second) << " is open with no candidate, setting to next Sprinter" << (*Sprinter)); 1632 Finder->second->T = BTS; 1609 Finder->second->T = BTS; // is last triangle 1633 1610 Finder->second->pointlist.push_back(Sprinter); 1634 1611 Finder->second->ShortestAngle = 0.; … … 1637 1614 } 1638 1615 } 1639 }; 1616 } 1617 ; 1640 1618 1641 1619 /** If a given \a *triangle is degenerated, this adds both sides. … … 1673 1651 // give some verbose output about the whole procedure 1674 1652 if (CandidateLine.T != NULL) 1675 LOG(2, " DEBUG: --> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector << ", from " << *CandidateLine.T << " and angle " << CandidateLine.ShortestAngle << ".");1653 LOG(2, "INFO: --> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector << ", from " << *CandidateLine.T << " and angle " << CandidateLine.ShortestAngle << "."); 1676 1654 else 1677 LOG(2, " DEBUG: --> New starting triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " and no top triangle.");1655 LOG(2, "INFO: --> New starting triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " and no top triangle."); 1678 1656 triangle = BTS; 1679 1657 … … 1895 1873 } 1896 1874 for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) { 1897 LOG(4, "DEBUG: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << ".");1875 LOG(4, "DEBUG: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << "."); 1898 1876 BaseLineNormal += (runner->second->NormalVector); 1899 1877 } … … 1935 1913 } 1936 1914 for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) { 1937 LOG( 1, "INFO: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << ".");1915 LOG(5, "INFO: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << "."); 1938 1916 BaseLineNormal += (runner->second->NormalVector); 1939 1917 } … … 1998 1976 m = 0; 1999 1977 // first obtain all triangle to delete ... (otherwise we pull the carpet (Base) from under the for-loop's feet) 2000 list 1978 list<BoundaryTriangleSet *> TrianglesOfBase; 2001 1979 for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); ++runner) 2002 1980 TrianglesOfBase.push_back(runner->second); 2003 1981 // .. then delete each triangle (which deletes the line as well) 2004 for (list 1982 for (list<BoundaryTriangleSet *>::iterator runner = TrianglesOfBase.begin(); !TrianglesOfBase.empty(); runner = TrianglesOfBase.begin()) { 2005 1983 LOG(3, "DEBUG: Deleting triangle " << *(*runner) << "."); 2006 1984 OldTriangleNrs[m++] = (*runner)->Nr; … … 2182 2160 TesselPoint *Candidate = NULL; 2183 2161 2184 LOG( 3, "DEBUG: NormalVector of BaseTriangle is " << NormalVector << ".");2162 LOG(4, "DEBUG: NormalVector of BaseTriangle is " << NormalVector << "."); 2185 2163 2186 2164 // copy old center … … 2190 2168 2191 2169 // construct center of circle 2192 CircleCenter = 0.5 * ((CandidateLine.BaseLine->endpoints[0]->node->getPosition()) + 2193 (CandidateLine.BaseLine->endpoints[1]->node->getPosition())); 2170 CircleCenter = 0.5 * ((CandidateLine.BaseLine->endpoints[0]->node->getPosition()) + (CandidateLine.BaseLine->endpoints[1]->node->getPosition())); 2194 2171 2195 2172 // construct normal vector of circle 2196 CirclePlaneNormal = (CandidateLine.BaseLine->endpoints[0]->node->getPosition()) - 2197 (CandidateLine.BaseLine->endpoints[1]->node->getPosition()); 2173 CirclePlaneNormal = (CandidateLine.BaseLine->endpoints[0]->node->getPosition()) - (CandidateLine.BaseLine->endpoints[1]->node->getPosition()); 2198 2174 2199 2175 RelativeOldSphereCenter = OldSphereCenter - CircleCenter; … … 2204 2180 CircleRadius = RADIUS * RADIUS - radius; 2205 2181 CirclePlaneNormal.Normalize(); 2206 LOG( 3, "DEBUG: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << ".");2182 LOG(4, "DEBUG: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "."); 2207 2183 2208 2184 // test whether old center is on the band's plane … … 2213 2189 radius = RelativeOldSphereCenter.NormSquared(); 2214 2190 if (fabs(radius - CircleRadius) < HULLEPSILON) { 2215 LOG( 3, "DEBUG: RelativeOldSphereCenter is at " << RelativeOldSphereCenter << ".");2191 LOG(4, "DEBUG: RelativeOldSphereCenter is at " << RelativeOldSphereCenter << "."); 2216 2192 2217 2193 // check SearchDirection 2218 LOG( 3, "DEBUG: SearchDirection is " << SearchDirection << ".");2194 LOG(4, "DEBUG: SearchDirection is " << SearchDirection << "."); 2219 2195 if (fabs(RelativeOldSphereCenter.ScalarProduct(SearchDirection)) > HULLEPSILON) { // rotated the wrong way! 2220 2196 ELOG(1, "SearchDirection and RelativeOldSphereCenter are not orthogonal!"); … … 2257 2233 // find center on the plane 2258 2234 GetCenterofCircumcircle(NewPlaneCenter, CandidateLine.BaseLine->endpoints[0]->node->getPosition(), CandidateLine.BaseLine->endpoints[1]->node->getPosition(), Candidate->getPosition()); 2259 LOG( 3, "DEBUG: NewPlaneCenter is " << NewPlaneCenter << ".");2235 LOG(5, "DEBUG: NewPlaneCenter is " << NewPlaneCenter << "."); 2260 2236 2261 2237 try { 2262 NewNormalVector = Plane((CandidateLine.BaseLine->endpoints[0]->node->getPosition()), 2263 (CandidateLine.BaseLine->endpoints[1]->node->getPosition()), 2264 (Candidate->getPosition())).getNormal(); 2265 LOG(3, "DEBUG: NewNormalVector is " << NewNormalVector << "."); 2238 NewNormalVector = Plane((CandidateLine.BaseLine->endpoints[0]->node->getPosition()), (CandidateLine.BaseLine->endpoints[1]->node->getPosition()), (Candidate->getPosition())).getNormal(); 2239 LOG(5, "DEBUG: NewNormalVector is " << NewNormalVector << "."); 2266 2240 radius = CandidateLine.BaseLine->endpoints[0]->node->DistanceSquared(NewPlaneCenter); 2267 LOG( 3, "DEBUG: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << ".");2268 LOG( 3, "DEBUG: SearchDirection is " << SearchDirection << ".");2269 LOG( 3, "DEBUG: Radius of CircumCenterCircle is " << radius << ".");2241 LOG(5, "DEBUG: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "."); 2242 LOG(5, "DEBUG: SearchDirection is " << SearchDirection << "."); 2243 LOG(5, "DEBUG: Radius of CircumCenterCircle is " << radius << "."); 2270 2244 if (radius < RADIUS * RADIUS) { 2271 2245 otherradius = CandidateLine.BaseLine->endpoints[1]->node->DistanceSquared(NewPlaneCenter); … … 2273 2247 // construct both new centers 2274 2248 NewSphereCenter = NewPlaneCenter; 2275 OtherNewSphereCenter = NewPlaneCenter;2249 OtherNewSphereCenter = NewPlaneCenter; 2276 2250 helper = NewNormalVector; 2277 2251 helper.Scale(sqrt(RADIUS * RADIUS - radius)); 2278 LOG( 4, "DEBUG: Distance of NewPlaneCenter " << NewPlaneCenter << " to either NewSphereCenter is " << helper.Norm() << " of vector " << helper << " with sphere radius " << RADIUS << ".");2252 LOG(5, "DEBUG: Distance of NewPlaneCenter " << NewPlaneCenter << " to either NewSphereCenter is " << helper.Norm() << " of vector " << helper << " with sphere radius " << RADIUS << "."); 2279 2253 NewSphereCenter += helper; 2280 LOG( 4, "DEBUG: NewSphereCenter is at " << NewSphereCenter << ".");2254 LOG(5, "DEBUG: NewSphereCenter is at " << NewSphereCenter << "."); 2281 2255 // OtherNewSphereCenter is created by the same vector just in the other direction 2282 2256 helper.Scale(-1.); 2283 2257 OtherNewSphereCenter += helper; 2284 LOG( 4, "DEBUG: OtherNewSphereCenter is at " << OtherNewSphereCenter << ".");2258 LOG(5, "DEBUG: OtherNewSphereCenter is at " << OtherNewSphereCenter << "."); 2285 2259 alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, NormalVector, SearchDirection, HULLEPSILON); 2286 2260 Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, NormalVector, SearchDirection, HULLEPSILON); … … 2303 2277 if ((CandidateLine.ShortestAngle - HULLEPSILON) < alpha) { 2304 2278 CandidateLine.pointlist.push_back(Candidate); 2305 LOG( 2, "ACCEPT: We have found an equally good candidate: " << *(Candidate) << " with " << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << ".");2279 LOG(3, "ACCEPT: We have found an equally good candidate: " << *(Candidate) << " with " << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << "."); 2306 2280 } else { 2307 2281 // remove all candidates from the list and then the list itself 2308 2282 CandidateLine.pointlist.clear(); 2309 2283 CandidateLine.pointlist.push_back(Candidate); 2310 LOG( 2, "ACCEPT: We have found a better candidate: " << *(Candidate) << " with " << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << ".");2284 LOG(3, "ACCEPT: We have found a better candidate: " << *(Candidate) << " with " << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << "."); 2311 2285 } 2312 2286 CandidateLine.ShortestAngle = alpha; 2313 LOG( 2, "DEBUG: There are " << CandidateLine.pointlist.size() << " candidates in the list now.");2287 LOG(4, "DEBUG: There are " << CandidateLine.pointlist.size() << " candidates in the list now."); 2314 2288 } else { 2315 2289 if ((Candidate != NULL) && (CandidateLine.pointlist.begin() != CandidateLine.pointlist.end())) { 2316 LOG( 3, "REJECT: Old candidate " << *(*CandidateLine.pointlist.begin()) << " with " << CandidateLine.ShortestAngle << " is better than new one " << *Candidate << " with " << alpha << " .");2290 LOG(4, "REJECT: Old candidate " << *(*CandidateLine.pointlist.begin()) << " with " << CandidateLine.ShortestAngle << " is better than new one " << *Candidate << " with " << alpha << " ."); 2317 2291 } else { 2318 LOG( 3, "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected.");2292 LOG(4, "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected."); 2319 2293 } 2320 2294 } 2321 2295 } else { 2322 ELOG(0, "REJECT: Distance to center of circumcircle is not the same from each corner of the triangle: " << fabs(radius - otherradius));2296 ASSERT(0, std::string("FindThirdPointForTesselation() - ") + std::string("REJECT: Distance to center of circumcircle is not the same from each corner of the triangle: ") + toString(fabs(radius - otherradius))); 2323 2297 } 2324 2298 } else { 2325 LOG( 3, "REJECT: NewSphereCenter " << NewSphereCenter << " for " << *Candidate << " is too far away: " << radius << ".");2299 LOG(4, "REJECT: NewSphereCenter " << NewSphereCenter << " for " << *Candidate << " is too far away: " << radius << "."); 2326 2300 } 2327 } 2328 catch (LinearDependenceException &excp){ 2329 LOG(3, boost::diagnostic_information(excp)); 2330 LOG(3, "REJECT: Three points from " << *CandidateLine.BaseLine << " and Candidate " << *Candidate << " are linear-dependent."); 2301 } catch (LinearDependenceException &excp) { 2302 LOG(4, boost::diagnostic_information(excp)); 2303 LOG(4, "REJECT: Three points from " << *CandidateLine.BaseLine << " and Candidate " << *Candidate << " are linear-dependent."); 2331 2304 } 2332 2305 } else { 2333 2306 if (ThirdPoint != NULL) { 2334 LOG( 3, "REJECT: Base triangle " << *CandidateLine.BaseLine << " and " << *ThirdPoint << " contains Candidate " << *Candidate << ".");2307 LOG(4, "REJECT: Base triangle " << *CandidateLine.BaseLine << " and " << *ThirdPoint << " contains Candidate " << *Candidate << "."); 2335 2308 } else { 2336 LOG( 3, "REJECT: Base triangle " << *CandidateLine.BaseLine << " contains Candidate " << *Candidate << ".");2309 LOG(4, "REJECT: Base triangle " << *CandidateLine.BaseLine << " contains Candidate " << *Candidate << "."); 2337 2310 } 2338 2311 } … … 2345 2318 } else { 2346 2319 if (ThirdPoint != NULL) 2347 LOG( 3, "Circumcircle for base line " << *CandidateLine.BaseLine << " and third node " << *ThirdPoint << " is too big!");2320 LOG(4, "DEBUG: Circumcircle for base line " << *CandidateLine.BaseLine << " and third node " << *ThirdPoint << " is too big!"); 2348 2321 else 2349 LOG( 3, "Circumcircle for base line " << *CandidateLine.BaseLine << " is too big!");2350 } 2351 2352 LOG( 2, "DEBUG: Sorting candidate list ...");2322 LOG(4, "DEBUG: Circumcircle for base line " << *CandidateLine.BaseLine << " is too big!"); 2323 } 2324 2325 LOG(5, "DEBUG: Sorting candidate list ..."); 2353 2326 if (CandidateLine.pointlist.size() > 1) { 2354 2327 CandidateLine.pointlist.unique(); … … 2356 2329 } 2357 2330 2358 ASSERT(CandidateLine.pointlist.empty() || (CandidateLine.CheckValidity(RADIUS, LC)), 2359 std::string("Tesselation::FindThirdPointForTesselation()") 2360 +std::string("There were other points contained in the rolling sphere as well!")); 2331 ASSERT(CandidateLine.pointlist.empty() || (CandidateLine.CheckValidity(RADIUS, LC)), std::string("Tesselation::FindThirdPointForTesselation()") + std::string("There were other points contained in the rolling sphere as well!")); 2361 2332 } 2362 2333 ; … … 2370 2341 { 2371 2342 //Info FunctionInfo(__func__); 2372 const BoundaryLineSet * lines[2] = { line1, line2};2343 const BoundaryLineSet * lines[2] = {line1, line2}; 2373 2344 class BoundaryPointSet *node = NULL; 2374 2345 PointMap OrderMap; … … 2377 2348 // for both lines 2378 2349 for (int j = 0; j < 2; j++) { // for both endpoints 2379 OrderTest = OrderMap.insert(pair<int, class BoundaryPointSet *> 2350 OrderTest = OrderMap.insert(pair<int, class BoundaryPointSet *>(lines[i]->endpoints[j]->Nr, lines[i]->endpoints[j])); 2380 2351 if (!OrderTest.second) { // if insertion fails, we have common endpoint 2381 2352 node = OrderTest.first->second; … … 2426 2397 // we should make sure that both triangles end up in the same entry 2427 2398 // in the distance multimap. Hence, we round to 6 digit precision. 2428 const double distance = 2429 1e-6*floor(FindPoint->second->node->DistanceSquared(x)*1e+6); 2399 const double distance = 1e-6 * floor(FindPoint->second->node->DistanceSquared(x) * 1e+6); 2430 2400 points->insert(DistanceToPointPair(distance, FindPoint->second)); 2431 2401 LOG(3, "DEBUG: Putting " << *FindPoint->second << " into the list."); … … 2472 2442 for (LineMap::iterator LineRunner = Runner->second->lines.begin(); LineRunner != Runner->second->lines.end(); LineRunner++) { 2473 2443 // calculate closest point on line to desired point 2474 helper = 0.5 * (((LineRunner->second)->endpoints[0]->node->getPosition()) + 2475 ((LineRunner->second)->endpoints[1]->node->getPosition())); 2444 helper = 0.5 * (((LineRunner->second)->endpoints[0]->node->getPosition()) + ((LineRunner->second)->endpoints[1]->node->getPosition())); 2476 2445 Center = (x) - helper; 2477 BaseLine = ((LineRunner->second)->endpoints[0]->node->getPosition()) - 2478 ((LineRunner->second)->endpoints[1]->node->getPosition()); 2446 BaseLine = ((LineRunner->second)->endpoints[0]->node->getPosition()) - ((LineRunner->second)->endpoints[1]->node->getPosition()); 2479 2447 Center.ProjectOntoPlane(BaseLine); 2480 2448 const double distance = Center.NormSquared(); … … 2533 2501 for (LineMap::iterator LineRunner = Runner->second->lines.begin(); LineRunner != Runner->second->lines.end(); LineRunner++) { 2534 2502 2535 BaseLine = ((LineRunner->second)->endpoints[0]->node->getPosition()) - 2536 ((LineRunner->second)->endpoints[1]->node->getPosition()); 2503 BaseLine = ((LineRunner->second)->endpoints[0]->node->getPosition()) - ((LineRunner->second)->endpoints[1]->node->getPosition()); 2537 2504 const double lengthBase = BaseLine.NormSquared(); 2538 2505 … … 2550 2517 MinDistance = lengthEnd; 2551 2518 LOG(1, "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[0]->node << " is closer with " << lengthEnd << "."); 2552 } else if (fabs(lengthEnd - MinDistance) < MYEPSILON) { // additional best candidate 2553 ClosestLines.insert(LineRunner->second); 2554 LOG(1, "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[1]->node << " is equally good with " << lengthEnd << "."); 2555 } else { // line is worse 2556 LOG(1, "REJECT: Line " << *LineRunner->second << " to either endpoints is further away than present closest line candidate: " << lengthEndA << ", " << lengthEndB << ", and distance is longer than baseline:" << lengthBase << "."); 2557 } 2519 } else 2520 if (fabs(lengthEnd - MinDistance) < MYEPSILON) { // additional best candidate 2521 ClosestLines.insert(LineRunner->second); 2522 LOG(1, "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[1]->node << " is equally good with " << lengthEnd << "."); 2523 } else { // line is worse 2524 LOG(1, "REJECT: Line " << *LineRunner->second << " to either endpoints is further away than present closest line candidate: " << lengthEndA << ", " << lengthEndB << ", and distance is longer than baseline:" << lengthBase << "."); 2525 } 2558 2526 } else { // intersection is closer, calculate 2559 2527 // calculate closest point on line to desired point … … 2681 2649 double distance = 0.; 2682 2650 2683 if (triangle == NULL) { // is boundary point or only point in point cloud?2651 if (triangle == NULL) { // is boundary point or only point in point cloud? 2684 2652 LOG(1, "No triangle given!"); 2685 2653 return -1.; … … 2790 2758 takePoint = true; 2791 2759 current = findLines->second->endpoints[1]->node; 2792 } else if (findLines->second->endpoints[1]->Nr == Point->getNr()) { 2793 takePoint = true; 2794 current = findLines->second->endpoints[0]->node; 2795 } 2760 } else 2761 if (findLines->second->endpoints[1]->Nr == Point->getNr()) { 2762 takePoint = true; 2763 current = findLines->second->endpoints[0]->node; 2764 } 2796 2765 2797 2766 if (takePoint) { … … 2833 2802 Vector OrthogonalVector; 2834 2803 Vector helper; 2835 const TesselPoint * const TrianglePoints[3] = { Point, NULL, NULL};2804 const TesselPoint * const TrianglePoints[3] = {Point, NULL, NULL}; 2836 2805 TriangleList *triangles = NULL; 2837 2806 … … 2844 2813 // calculate central point 2845 2814 triangles = FindTriangles(TrianglePoints); 2846 ASSERT((triangles == NULL) || (triangles->empty()), 2847 std::string("Tesselation::GetCircleOfConnectedTriangles()") 2848 +std::string("Could not find any triangles for point "+toString(*Point)+".")); 2815 ASSERT((triangles == NULL) || (triangles->empty()), std::string("Tesselation::GetCircleOfConnectedTriangles()") + std::string("Could not find any triangles for point " + toString(*Point) + ".")); 2849 2816 for (TriangleList::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++) 2850 2817 PlaneNormal += (*Runner)->NormalVector; … … 2860 2827 AngleZero = ((*SetOfNeighbours->begin())->getPosition()) - (Point->getPosition()); 2861 2828 AngleZero.ProjectOntoPlane(PlaneNormal); 2862 ASSERT (AngleZero.NormSquared() > MYEPSILON, 2863 std::string("Tesselation::GetCircleOfConnectedTriangles() - ") 2864 +std::string("AngleZero is 0 even with alternative reference.") 2865 +std::string("The algorithm has to be changed here!")); 2829 ASSERT(AngleZero.NormSquared() > MYEPSILON, std::string("Tesselation::GetCircleOfConnectedTriangles() - ") + std::string("AngleZero is 0 even with alternative reference.") + std::string("The algorithm has to be changed here!")); 2866 2830 } 2867 2831 LOG(4, "DEBUG: Reference vector on this plane representing angle 0 is " << AngleZero << "."); 2868 2832 if (AngleZero.NormSquared() > MYEPSILON) 2869 OrthogonalVector = Plane(PlaneNormal, AngleZero, 0).getNormal();2833 OrthogonalVector = Plane(PlaneNormal, AngleZero, 0).getNormal(); 2870 2834 else 2871 2835 OrthogonalVector.MakeNormalTo(PlaneNormal); … … 2878 2842 double angle = GetAngle(helper, AngleZero, OrthogonalVector); 2879 2843 LOG(4, "DEBUG" << angle << " for point " << **listRunner << "."); 2880 anglesOfPoints.insert(pair<double, TesselPoint*> 2844 anglesOfPoints.insert(pair<double, TesselPoint*>(angle, (*listRunner))); 2881 2845 } 2882 2846 … … 2933 2897 int counter = 0; 2934 2898 while (TesselC != SetOfNeighbours->end()) { 2935 helper = Plane(((*TesselA)->getPosition()), 2936 ((*TesselB)->getPosition()), 2937 ((*TesselC)->getPosition())).getNormal(); 2899 helper = Plane(((*TesselA)->getPosition()), ((*TesselB)->getPosition()), ((*TesselC)->getPosition())).getNormal(); 2938 2900 LOG(5, "DEBUG: Making normal vector out of " << *(*TesselA) << ", " << *(*TesselB) << " and " << *(*TesselC) << ":" << helper); 2939 2901 counter++; … … 2944 2906 } 2945 2907 //LOG(0, "Summed vectors " << center << "; number of points " << connectedPoints.size() << "; scale factor " << counter); 2946 PlaneNormal.Scale(1.0 / (double) 2908 PlaneNormal.Scale(1.0 / (double)counter); 2947 2909 // LOG(1, "INFO: Calculated center of all circle points is " << center << "."); 2948 2910 // … … 2958 2920 AngleZero.ProjectOntoPlane(PlaneNormal); 2959 2921 } 2960 if ((Reference.IsZero()) || (AngleZero.NormSquared() < MYEPSILON 2922 if ((Reference.IsZero()) || (AngleZero.NormSquared() < MYEPSILON)) { 2961 2923 LOG(4, "DEBUG: Using alternatively " << (*SetOfNeighbours->begin())->getPosition() << " as angle 0 referencer."); 2962 2924 AngleZero = ((*SetOfNeighbours->begin())->getPosition()) - (Point->getPosition()); 2963 2925 AngleZero.ProjectOntoPlane(PlaneNormal); 2964 ASSERT(AngleZero.NormSquared() > MYEPSILON, 2965 std::string("Tesselation::GetCircleOfSetOfPoints() - ") 2966 +std::string("AngleZero is 0 even with alternative reference.") 2967 +std::string("The algorithm has to be changed here!")); 2926 ASSERT(AngleZero.NormSquared() > MYEPSILON, std::string("Tesselation::GetCircleOfSetOfPoints() - ") + std::string("AngleZero is 0 even with alternative reference.") + std::string("The algorithm has to be changed here!")); 2968 2927 } 2969 2928 LOG(4, "DEBUG: Reference vector on this plane representing angle 0 is " << AngleZero << "."); 2970 2929 if (AngleZero.NormSquared() > MYEPSILON) 2971 OrthogonalVector = Plane(PlaneNormal, AngleZero, 0).getNormal();2930 OrthogonalVector = Plane(PlaneNormal, AngleZero, 0).getNormal(); 2972 2931 else 2973 2932 OrthogonalVector.MakeNormalTo(PlaneNormal); … … 2983 2942 angle = 2. * M_PI - angle; 2984 2943 LOG(4, "DEBUG: Calculated angle between " << helper << " and " << AngleZero << " is " << angle << " for point " << **listRunner << "."); 2985 InserterTest = anglesOfPoints.insert(pair<double, TesselPoint*> (angle, (*listRunner))); 2986 ASSERT(InserterTest.second, 2987 std::string("Tesselation::GetCircleOfSetOfPoints() - ") 2988 +std::string("got two atoms with same angle "+toString(*((InserterTest.first)->second))) 2989 +std::string(" and "+toString((*listRunner)))); 2944 InserterTest = anglesOfPoints.insert(pair<double, TesselPoint*>(angle, (*listRunner))); 2945 ASSERT(InserterTest.second, std::string("Tesselation::GetCircleOfSetOfPoints() - ") + std::string("got two atoms with same angle " + toString(*((InserterTest.first)->second))) + std::string(" and " + toString((*listRunner)))); 2990 2946 } 2991 2947 … … 3007 2963 //Info FunctionInfo(__func__); 3008 2964 map<double, TesselPoint*> anglesOfPoints; 3009 list<TesselPointList *> *ListOfPaths = new list<TesselPointList *> 2965 list<TesselPointList *> *ListOfPaths = new list<TesselPointList *>; 3010 2966 TesselPointList *connectedPath = NULL; 3011 2967 Vector center; … … 3033 2989 map<class BoundaryTriangleSet *, bool>::iterator TriangleRunner; 3034 2990 for (LineMap::iterator Runner = ReferencePoint->lines.begin(); Runner != ReferencePoint->lines.end(); Runner++) { 3035 TouchedLine.insert(pair<class BoundaryLineSet *, bool> (Runner->second, false)); 3036 for (TriangleMap::iterator Sprinter = Runner->second->triangles.begin(); Sprinter != Runner->second->triangles.end(); Sprinter++) 3037 TouchedTriangle.insert(pair<class BoundaryTriangleSet *, bool> (Sprinter->second, false)); 2991 LOG(4, "DEBUG: Adding " << *Runner->second << " to TouchedLine map."); 2992 TouchedLine.insert(pair<class BoundaryLineSet *, bool>(Runner->second, false)); 2993 for (TriangleMap::iterator Sprinter = Runner->second->triangles.begin(); Sprinter != Runner->second->triangles.end(); Sprinter++) { 2994 LOG(4, "DEBUG: Adding " << *Sprinter->second << " to TouchedTriangle map."); 2995 TouchedTriangle.insert(pair<class BoundaryTriangleSet *, bool>(Sprinter->second, false)); 2996 } 3038 2997 } 3039 2998 if (!ReferencePoint->lines.empty()) { … … 3042 3001 if (LineRunner == TouchedLine.end()) { 3043 3002 ELOG(1, "I could not find " << *runner->second << " in the touched list."); 3044 } else if (!LineRunner->second) { 3045 LineRunner->second = true; 3046 connectedPath = new TesselPointList; 3047 triangle = NULL; 3048 CurrentLine = runner->second; 3049 StartLine = CurrentLine; 3050 CurrentPoint = CurrentLine->GetOtherEndpoint(ReferencePoint); 3051 LOG(1, "INFO: Beginning path retrieval at " << *CurrentPoint << " of line " << *CurrentLine << "."); 3052 do { 3053 // push current one 3054 LOG(1, "INFO: Putting " << *CurrentPoint << " at end of path."); 3055 connectedPath->push_back(CurrentPoint->node); 3056 3057 // find next triangle 3058 for (TriangleMap::iterator Runner = CurrentLine->triangles.begin(); Runner != CurrentLine->triangles.end(); Runner++) { 3059 LOG(1, "INFO: Inspecting triangle " << *Runner->second << "."); 3060 if ((Runner->second != triangle)) { // look for first triangle not equal to old one 3061 triangle = Runner->second; 3062 TriangleRunner = TouchedTriangle.find(triangle); 3063 if (TriangleRunner != TouchedTriangle.end()) { 3064 if (!TriangleRunner->second) { 3065 TriangleRunner->second = true; 3066 LOG(1, "INFO: Connecting triangle is " << *triangle << "."); 3067 break; 3003 } else 3004 if (!LineRunner->second) { 3005 LineRunner->second = true; 3006 connectedPath = new TesselPointList; 3007 triangle = NULL; 3008 CurrentLine = runner->second; 3009 StartLine = CurrentLine; 3010 CurrentPoint = CurrentLine->GetOtherEndpoint(ReferencePoint); 3011 LOG(3, "INFO: Beginning path retrieval at " << *CurrentPoint << " of line " << *CurrentLine << "."); 3012 do { 3013 // push current one 3014 LOG(3, "INFO: Putting " << *CurrentPoint << " at end of path."); 3015 connectedPath->push_back(CurrentPoint->node); 3016 3017 // find next triangle 3018 for (TriangleMap::iterator Runner = CurrentLine->triangles.begin(); Runner != CurrentLine->triangles.end(); Runner++) { 3019 LOG(4, "DEBUG: Inspecting triangle " << *Runner->second << "."); 3020 if ((Runner->second != triangle)) { // look for first triangle not equal to old one 3021 triangle = Runner->second; 3022 TriangleRunner = TouchedTriangle.find(triangle); 3023 if (TriangleRunner != TouchedTriangle.end()) { 3024 if (!TriangleRunner->second) { 3025 TriangleRunner->second = true; 3026 LOG(4, "DEBUG: Connecting triangle is " << *triangle << "."); 3027 break; 3028 } else { 3029 LOG(4, "DEBUG: Skipping " << *triangle << ", as we have already visited it."); 3030 triangle = NULL; 3031 } 3068 3032 } else { 3069 LOG(1, "INFO: Skipping " << *triangle << ", as we have already visited it.");3033 ELOG(1, "I could not find " << *triangle << " in the touched list."); 3070 3034 triangle = NULL; 3071 3035 } 3072 3036 } else { 3073 ELOG(1, "I could not find " << *triangle << " in the touched list.");3037 // as we have stumbled upon the same triangle, we don't need the check anymore 3074 3038 triangle = NULL; 3075 3039 } 3076 } else {3077 // as we have stumbled upon the same triangle, we don't need the check anymore3078 triangle = NULL;3079 3040 } 3041 if (triangle == NULL) 3042 break; 3043 // find next line 3044 for (int i = 0; i < 3; i++) { 3045 if ((triangle->lines[i] != CurrentLine) && (triangle->lines[i]->ContainsBoundaryPoint(ReferencePoint))) { // not the current line and still containing Point 3046 CurrentLine = triangle->lines[i]; 3047 LOG(3, "INFO: Connecting line is " << *CurrentLine << "."); 3048 break; 3049 } 3050 } 3051 LineRunner = TouchedLine.find(CurrentLine); 3052 if (LineRunner == TouchedLine.end()) 3053 ELOG(1, "I could not find " << *CurrentLine << " in the touched list."); 3054 else 3055 LineRunner->second = true; 3056 // find next point 3057 CurrentPoint = CurrentLine->GetOtherEndpoint(ReferencePoint); 3058 3059 } while (CurrentLine != StartLine); 3060 // last point is missing, as it's on start line 3061 if (StartLine->GetOtherEndpoint(ReferencePoint)->node != connectedPath->back()) { 3062 LOG(3, "INFO: Putting " << *CurrentPoint << " at end of path to close it."); 3063 connectedPath->push_back(StartLine->GetOtherEndpoint(ReferencePoint)->node); 3080 3064 } 3081 if (triangle == NULL) 3082 break; 3083 // find next line 3084 for (int i = 0; i < 3; i++) { 3085 if ((triangle->lines[i] != CurrentLine) && (triangle->lines[i]->ContainsBoundaryPoint(ReferencePoint))) { // not the current line and still containing Point 3086 CurrentLine = triangle->lines[i]; 3087 LOG(1, "INFO: Connecting line is " << *CurrentLine << "."); 3088 break; 3089 } 3090 } 3091 LineRunner = TouchedLine.find(CurrentLine); 3092 if (LineRunner == TouchedLine.end()) 3093 ELOG(1, "I could not find " << *CurrentLine << " in the touched list."); 3094 else 3095 LineRunner->second = true; 3096 // find next point 3097 CurrentPoint = CurrentLine->GetOtherEndpoint(ReferencePoint); 3098 3099 } while (CurrentLine != StartLine); 3100 // last point is missing, as it's on start line 3101 LOG(1, "INFO: Putting " << *CurrentPoint << " at end of path."); 3102 if (StartLine->GetOtherEndpoint(ReferencePoint)->node != connectedPath->back()) 3103 connectedPath->push_back(StartLine->GetOtherEndpoint(ReferencePoint)->node); 3104 3105 ListOfPaths->push_back(connectedPath); 3106 } else { 3107 LOG(1, "INFO: Skipping " << *runner->second << ", as we have already visited it."); 3108 } 3065 3066 ListOfPaths->push_back(connectedPath); 3067 } else { 3068 LOG(3, "DEBUG: Skipping " << *runner->second << ", as we have already visited it."); 3069 } 3109 3070 } 3110 3071 } else { … … 3125 3086 //Info FunctionInfo(__func__); 3126 3087 list<TesselPointList *> *ListofPaths = GetPathsOfConnectedPoints(Point); 3127 list<TesselPointList *> *ListofClosedPaths = new list<TesselPointList *> 3088 list<TesselPointList *> *ListofClosedPaths = new list<TesselPointList *>; 3128 3089 TesselPointList *connectedPath = NULL; 3129 3090 TesselPointList *newPath = NULL; … … 3135 3096 connectedPath = *ListRunner; 3136 3097 3137 LOG(1, "INFO: Current path is " << connectedPath << "."); 3098 if (DoLog(2)) { 3099 std::stringstream output; 3100 output << "INFO: Current path is "; 3101 BOOST_FOREACH(const TesselPoint * const item, *connectedPath) { 3102 output << *item << " "; 3103 } 3104 LOG(1, output.str()); 3105 } 3138 3106 3139 3107 // go through list, look for reappearance of starting Point and count … … 3144 3112 if ((*CircleRunner == *CircleStart) && (CircleRunner != CircleStart)) { // is not the very first point 3145 3113 // we have a closed circle from Marker to new Marker 3146 if (DoLog( 1)) {3114 if (DoLog(3)) { 3147 3115 std::stringstream output; 3148 output << count + 1 << ". closed path consists of: "; 3149 for (TesselPointList::iterator CircleSprinter = Marker; 3150 CircleSprinter != CircleRunner; 3151 CircleSprinter++) 3116 output << "DEBUG: " << count + 1 << ". closed path consists of: "; 3117 for (TesselPointList::iterator CircleSprinter = Marker; CircleSprinter != CircleRunner; CircleSprinter++) 3152 3118 output << (**CircleSprinter) << " <-> "; 3153 3119 LOG(1, output.str()); … … 3165 3131 } 3166 3132 } 3167 LOG( 1, "INFO: " << count << " closed additional path(s) have been created.");3133 LOG(2, "DEBUG: " << count << " closed additional path(s) have been created."); 3168 3134 3169 3135 // delete list of paths … … 3206 3172 struct CloserToPiHalf 3207 3173 { 3208 bool operator()(double angle, double smallestangle) { 3209 return (fabs(angle - M_PI/2.) < fabs(smallestangle - M_PI/2.)); 3174 bool operator()(double angle, double smallestangle) 3175 { 3176 return (fabs(angle - M_PI / 2.) < fabs(smallestangle - M_PI / 2.)); 3210 3177 } 3211 3178 }; … … 3264 3231 NormalVector.Zero(); 3265 3232 for (TriangleMap::iterator Runner = Candidates.begin(); Runner != Candidates.end(); Runner++) { 3266 LOG( 1, "INFO: Removing triangle " << *(Runner->second) << ".");3233 LOG(3, "DEBUG: Removing triangle " << *(Runner->second) << "."); 3267 3234 NormalVector -= Runner->second->NormalVector; // has to point inward 3268 3235 RemoveTesselationTriangle(Runner->second); 3269 3236 count++; 3270 3237 } 3271 LOG( 1,count << " triangles were removed.");3238 LOG(2, "INFO: " << count << " triangles were removed."); 3272 3239 3273 3240 list<TesselPointList *>::iterator ListAdvance = ListOfClosedPaths->begin(); … … 3292 3259 EndNode = connectedPath->end(); 3293 3260 for (MiddleNode = connectedPath->begin(); MiddleNode != connectedPath->end(); MiddleNode++) { 3294 LOG( 1, "INFO: MiddleNode is " << **MiddleNode << ".");3261 LOG(3, "INFO: MiddleNode is " << **MiddleNode << "."); 3295 3262 // construct vectors to next and previous neighbour 3296 3263 StartNode = MiddleNode; … … 3308 3275 OrthogonalVector = ((*MiddleNode)->getPosition()) - OldPoint; 3309 3276 OrthogonalVector.MakeNormalTo(Reference); 3310 angles.push_back( GetAngle(Point, Reference, OrthogonalVector));3277 angles.push_back(GetAngle(Point, Reference, OrthogonalVector)); 3311 3278 } 3312 3279 const angles_t::const_iterator maxiter = std::max_element(angles.begin(), angles.end()); … … 3321 3288 } 3322 3289 EndNode = connectedPath->begin(); 3323 std::advance(EndNode, 3324 std::distance(const_cast<const angles_t &>(angles).begin(), miniter)); 3290 std::advance(EndNode, std::distance(const_cast<const angles_t &>(angles).begin(), miniter)); 3325 3291 angles.clear(); 3326 3292 3327 3293 MiddleNode = EndNode; 3328 ASSERT( MiddleNode == connectedPath->end(), 3329 "Tesselation::RemovePointFromTesselatedSurface() - Could not find a smallest angle!"); 3294 ASSERT(MiddleNode == connectedPath->end(), "Tesselation::RemovePointFromTesselatedSurface() - Could not find a smallest angle!"); 3330 3295 StartNode = MiddleNode; 3331 3296 if (StartNode == connectedPath->begin()) … … 3361 3326 continue; 3362 3327 } 3363 LOG(3, " Adding new triangle points.");3328 LOG(3, "DEBUG: Adding new triangle points."); 3364 3329 AddTesselationPoint(*StartNode, 0); 3365 3330 AddTesselationPoint(*MiddleNode, 1); 3366 3331 AddTesselationPoint(*EndNode, 2); 3367 LOG(3, " Adding new triangle lines.");3332 LOG(3, "DEBUG: Adding new triangle lines."); 3368 3333 AddTesselationLine(NULL, NULL, TPS[0], TPS[1], 0); 3369 3334 // line between start and end must be new (except for very last triangle) … … 3381 3346 // prepare nodes for next triangle 3382 3347 StartNode = EndNode; 3383 LOG( 2, "Removing " << **MiddleNode << " from closed path, remaining points: " << connectedPath->size() << ".");3348 LOG(3, "DEBUG: Removing " << **MiddleNode << " from closed path, remaining points: " << connectedPath->size() << "."); 3384 3349 connectedPath->remove(*MiddleNode); // remove the middle node (it is surrounded by triangles) 3385 ASSERT(connectedPath->size() >= 2, 3386 "Tesselation::RemovePointFromTesselatedSurface() - There are only two endpoints left!"); 3350 ASSERT(connectedPath->size() >= 2, "Tesselation::RemovePointFromTesselatedSurface() - There are only two endpoints left!"); 3387 3351 if (connectedPath->size() == 2) { // we are done 3388 3352 connectedPath->remove(*StartNode); // remove the start node … … 3401 3365 } 3402 3366 // maximize the inner lines (we preferentially created lines with a huge angle, which is for the tesselation not wanted though useful for the closing) 3367 LOG(3, "INFO: Flipping inner lines to maximize volume"); 3403 3368 if (NewLines.size() > 1) { 3404 3369 LineList::iterator Candidate; … … 3416 3381 if (maxgain != 0) { 3417 3382 volume += maxgain; 3418 LOG(1, "Flipping baseline with highest volume" << **Candidate << "."); 3383 LOG(3, "DEBUG: Flipping baseline with highest volume gain of " 3384 << maxgain << ": " << **Candidate << "."); 3419 3385 OtherBase = FlipBaseline(*Candidate); 3420 3386 NewLines.erase(Candidate); … … 3483 3449 if (TrianglePoints[j] != NULL) { 3484 3450 for (FindLine = TrianglePoints[i]->lines.find(TrianglePoints[j]->node->getNr()); // is a multimap! 3485 (FindLine != TrianglePoints[i]->lines.end()) && (FindLine->first == TrianglePoints[j]->node->getNr()); FindLine++) {3451 (FindLine != TrianglePoints[i]->lines.end()) && (FindLine->first == TrianglePoints[j]->node->getNr()); FindLine++) { 3486 3452 for (FindTriangle = FindLine->second->triangles.begin(); FindTriangle != FindLine->second->triangles.end(); FindTriangle++) { 3487 3453 if (FindTriangle->second->IsPresentTupel(TrianglePoints)) { … … 3504 3470 break; 3505 3471 for (FindLine = TrianglePoints[(i + 1) % 3]->lines.find(TrianglePoints[(i + 2) % 3]->node->getNr()); // is a multimap! 3506 (FindLine != TrianglePoints[(i + 1) % 3]->lines.end()) && (FindLine->first == TrianglePoints[(i + 2) % 3]->node->getNr()); FindLine++) {3472 (FindLine != TrianglePoints[(i + 1) % 3]->lines.end()) && (FindLine->first == TrianglePoints[(i + 2) % 3]->node->getNr()); FindLine++) { 3507 3473 for (FindTriangle = FindLine->second->triangles.begin(); FindTriangle != FindLine->second->triangles.end(); FindTriangle++) { 3508 3474 if (FindTriangle->second->IsPresentTupel(TrianglePoints)) { … … 3533 3499 } 3534 3500 default: 3535 ASSERT(0, 3536 "Tesselation::FindTriangles() - Number of wildcards is greater than 3, cannot happen!"); 3501 ASSERT(0, "Tesselation::FindTriangles() - Number of wildcards is greater than 3, cannot happen!"); 3537 3502 break; 3538 3503 } … … 3560 3525 if (a->endpoints[lowerNra] < b->endpoints[lowerNrb]) 3561 3526 return true; 3562 else if (a->endpoints[lowerNra] > b->endpoints[lowerNrb]) 3563 return false; 3564 else { // both lower-numbered endpoints are the same ... 3565 if (a->endpoints[(lowerNra + 1) % 2] < b->endpoints[(lowerNrb + 1) % 2]) 3566 return true; 3567 else if (a->endpoints[(lowerNra + 1) % 2] > b->endpoints[(lowerNrb + 1) % 2]) 3527 else 3528 if (a->endpoints[lowerNra] > b->endpoints[lowerNrb]) 3568 3529 return false; 3569 } 3530 else { // both lower-numbered endpoints are the same ... 3531 if (a->endpoints[(lowerNra + 1) % 2] < b->endpoints[(lowerNrb + 1) % 2]) 3532 return true; 3533 else 3534 if (a->endpoints[(lowerNra + 1) % 2] > b->endpoints[(lowerNrb + 1) % 2]) 3535 return false; 3536 } 3570 3537 return false; 3571 3538 } … … 3588 3555 3589 3556 // sanity check 3590 if (LinesOnBoundary.empty()) {3591 ELOG(2, "FindAllDegeneratedTriangles() was called without any tesselation structure.");3592 return DegeneratedLines;3593 }3594 LineMap::iterator LineRunner1;3595 pair<UniqueLines::iterator, bool> tester;3557 if (LinesOnBoundary.empty()) { 3558 ELOG(2, "FindAllDegeneratedTriangles() was called without any tesselation structure."); 3559 return DegeneratedLines; 3560 } 3561 LineMap::iterator LineRunner1; 3562 pair<UniqueLines::iterator, bool> tester; 3596 3563 for (LineRunner1 = LinesOnBoundary.begin(); LineRunner1 != LinesOnBoundary.end(); ++LineRunner1) { 3597 3564 tester = AllLines.insert(LineRunner1->second); … … 3610 3577 const LineMap::const_iterator Line2 = LinesOnBoundary.find((*it).second); 3611 3578 if (Line1 != LinesOnBoundary.end() && Line2 != LinesOnBoundary.end()) 3612 3579 LOG(3, "DEBUG: " << *Line1->second << " => " << *Line2->second); 3613 3580 else 3614 3581 ELOG(1, "Either " << (*it).first << " or " << (*it).second << " are not in LinesOnBoundary!"); 3615 3582 } 3616 3583 … … 3644 3611 for (TriangleRunner2 = line2->triangles.begin(); TriangleRunner2 != line2->triangles.end(); ++TriangleRunner2) { 3645 3612 if ((TriangleRunner1->second != TriangleRunner2->second) && (TriangleRunner1->second->IsPresentTupel(TriangleRunner2->second))) { 3646 DegeneratedTriangles->insert(pair<int, int> 3647 DegeneratedTriangles->insert(pair<int, int> 3613 DegeneratedTriangles->insert(pair<int, int>(TriangleRunner1->second->Nr, TriangleRunner2->second->Nr)); 3614 DegeneratedTriangles->insert(pair<int, int>(TriangleRunner2->second->Nr, TriangleRunner1->second->Nr)); 3648 3615 } 3649 3616 } … … 3713 3680 // OtherpartnerTriangle = TriangleRunner->second; 3714 3681 // } 3715 /// interchanges their lines so that triangle->lines[i] == partnerTriangle->lines[j]3716 // the line of triangle receives the degenerated ones3717 triangle->lines[i]->triangles.erase(Othertriangle->Nr);3682 /// interchanges their lines so that triangle->lines[i] == partnerTriangle->lines[j] 3683 // the line of triangle receives the degenerated ones 3684 triangle->lines[i]->triangles.erase(Othertriangle->Nr); 3718 3685 triangle->lines[i]->triangles.insert(TrianglePair(partnerTriangle->Nr, partnerTriangle)); 3719 3686 for (int k = 0; k < 3; k++) … … 3729 3696 3730 3697 // erase the pair 3731 count += (int) 3698 count += (int)DegeneratedTriangles->erase(triangle->Nr); 3732 3699 LOG(4, "DEBUG: RemoveDegeneratedTriangles() removes triangle " << *triangle << "."); 3733 3700 RemoveTesselationTriangle(triangle); 3734 count += (int) 3701 count += (int)DegeneratedTriangles->erase(partnerTriangle->Nr); 3735 3702 LOG(4, "DEBUG: RemoveDegeneratedTriangles() removes triangle " << *partnerTriangle << "."); 3736 3703 RemoveTesselationTriangle(partnerTriangle); … … 3780 3747 class BoundaryLineSet *BestLine = NULL; 3781 3748 for (LineMap::iterator Runner = NearestBoundaryPoint->lines.begin(); Runner != NearestBoundaryPoint->lines.end(); Runner++) { 3782 BaseLine = (Runner->second->endpoints[0]->node->getPosition()) - 3783 (Runner->second->endpoints[1]->node->getPosition()); 3784 CenterToPoint = 0.5 * ((Runner->second->endpoints[0]->node->getPosition()) + 3785 (Runner->second->endpoints[1]->node->getPosition())); 3749 BaseLine = (Runner->second->endpoints[0]->node->getPosition()) - (Runner->second->endpoints[1]->node->getPosition()); 3750 CenterToPoint = 0.5 * ((Runner->second->endpoints[0]->node->getPosition()) + (Runner->second->endpoints[1]->node->getPosition())); 3786 3751 CenterToPoint -= (point->getPosition()); 3787 3752 angle = CenterToPoint.Angle(BaseLine); 3788 if (fabs(angle - M_PI /2.) < fabs(BestAngle - M_PI/2.)) {3753 if (fabs(angle - M_PI / 2.) < fabs(BestAngle - M_PI / 2.)) { 3789 3754 BestAngle = angle; 3790 3755 BestLine = Runner->second; … … 3835 3800 for (int i = 0; i < 3; i++) { // look for the same line as BestLine (only it's its degenerated companion) 3836 3801 if ((BTS->lines[i]->ContainsBoundaryPoint(BestLine->endpoints[0])) && (BTS->lines[i]->ContainsBoundaryPoint(BestLine->endpoints[1]))) { 3837 ASSERT(BestLine != BTS->lines[i], 3838 std::string("Tesselation::AddBoundaryPointByDegeneratedTriangle() - ") 3839 +std::string("BestLine is same as found line, something's wrong here!")); 3840 BTS->lines[i]->triangles.insert(pair<int, class BoundaryTriangleSet *> (TempTriangle->Nr, TempTriangle)); 3802 ASSERT(BestLine != BTS->lines[i], std::string("Tesselation::AddBoundaryPointByDegeneratedTriangle() - ") + std::string("BestLine is same as found line, something's wrong here!")); 3803 BTS->lines[i]->triangles.insert(pair<int, class BoundaryTriangleSet *>(TempTriangle->Nr, TempTriangle)); 3841 3804 TempTriangle->lines[nr] = BTS->lines[i]; 3842 3805 break; … … 3860 3823 if (LastTriangle != NULL) { 3861 3824 stringstream sstr; 3862 sstr << "-" << TrianglesOnBoundary.size() << "-" << LastTriangle->getEndpointName(0) << "_" << LastTriangle->getEndpointName(1) << "_" << LastTriangle->getEndpointName(2);3825 sstr << "-" << TrianglesOnBoundary.size() << "-" << LastTriangle->getEndpointName(0) << "_" << LastTriangle->getEndpointName(1) << "_" << LastTriangle->getEndpointName(2); 3863 3826 NumberName = sstr.str(); 3864 3827 if (DoTecplotOutput) { … … 3902 3865 if (s1->endpoints.size() < s2->endpoints.size()) 3903 3866 return true; 3904 else if (s1->endpoints.size() > s2->endpoints.size()) 3905 return false; 3906 else { // equality of number of endpoints 3907 PointSet::const_iterator Walker1 = s1->endpoints.begin(); 3908 PointSet::const_iterator Walker2 = s2->endpoints.begin(); 3909 while ((Walker1 != s1->endpoints.end()) || (Walker2 != s2->endpoints.end())) { 3910 if ((*Walker1)->Nr < (*Walker2)->Nr) 3911 return true; 3912 else if ((*Walker1)->Nr > (*Walker2)->Nr) 3913 return false; 3914 Walker1++; 3915 Walker2++; 3867 else 3868 if (s1->endpoints.size() > s2->endpoints.size()) 3869 return false; 3870 else { // equality of number of endpoints 3871 PointSet::const_iterator Walker1 = s1->endpoints.begin(); 3872 PointSet::const_iterator Walker2 = s2->endpoints.begin(); 3873 while ((Walker1 != s1->endpoints.end()) || (Walker2 != s2->endpoints.end())) { 3874 if ((*Walker1)->Nr < (*Walker2)->Nr) 3875 return true; 3876 else 3877 if ((*Walker1)->Nr > (*Walker2)->Nr) 3878 return false; 3879 Walker1++; 3880 Walker2++; 3881 } 3882 return false; 3916 3883 } 3917 return false;3918 }3919 3884 } 3920 3885 }; … … 3941 3906 for (TriangleMap::const_iterator TriangleRunner = (LineRunner->second)->triangles.begin(); TriangleRunner != (LineRunner->second)->triangles.end(); TriangleRunner++) { 3942 3907 if (DegeneratedTriangles->find(TriangleRunner->second->Nr) == DegeneratedTriangles->end()) { 3943 TriangleInsertionTester = TriangleVectors.insert(pair<int, Vector *> 3908 TriangleInsertionTester = TriangleVectors.insert(pair<int, Vector *>((TriangleRunner->second)->Nr, &((TriangleRunner->second)->NormalVector))); 3944 3909 if (TriangleInsertionTester.second) 3945 3910 LOG(5, "DEBUG: Adding triangle " << *(TriangleRunner->second) << " to triangles to check-list."); … … 3954 3919 if (VectorWalker != VectorRunner) { // skip equals 3955 3920 const double SCP = VectorWalker->second->ScalarProduct(*VectorRunner->second); // ScalarProduct should result in -1. for degenerated triangles 3956 LOG(4, "DEBUG: Checking " << * VectorWalker->second << " against " << *VectorRunner->second<< ": " << SCP);3921 LOG(4, "DEBUG: Checking " << *(VectorWalker->second) << " against " << *(VectorRunner->second) << ": " << SCP); 3957 3922 if (fabs(SCP + 1.) < ParallelEpsilon) { 3958 3923 InsertionTester = EndpointCandidateList.insert((Runner->second)); … … 4036 4001 // The Problem is probably due to two degenerated polygons being connected by a bridging, non-degenerated polygon, as somehow one node has 4037 4002 // connections to either polygon ... 4038 ASSERT (T->size() % 2 == 0, 4039 std::string("Tesselation::CorrectAllDegeneratedPolygons() - ") 4040 +std::string(" degenerated polygon contains an odd number of triangles,") 4041 +std::string(" probably contains bridging non-degenerated ones, too!")); 4003 ASSERT(T->size() % 2 == 0, std::string("Tesselation::CorrectAllDegeneratedPolygons() - ") + std::string(" degenerated polygon contains an odd number of triangles,") + std::string(" probably contains bridging non-degenerated ones, too!")); 4042 4004 TriangleSet::iterator TriangleWalker = T->begin(); // is the inner iterator 4043 4005 /// 4a. Get NormalVector for one side (this is "front") -
src/Tesselation/tesselationhelpers.cpp
rdfe054 rb8d215 235 235 if (helper.ScalarProduct(SearchDirection) < -HULLEPSILON) // acos is not unique on [0, 2.*M_PI), hence extra check to decide between two half intervals 236 236 alpha = 2.*M_PI - alpha; 237 LOG( 3, "DEBUG: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << RelativeOldSphereCenter << " and resulting angle is " << alpha << ".");237 LOG(5, "DEBUG: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << RelativeOldSphereCenter << " and resulting angle is " << alpha << "."); 238 238 radius = helper.distance(RelativeOldSphereCenter); 239 239 helper.ProjectOntoPlane(NormalVector); 240 240 // check whether new center is somewhat away or at least right over the current baseline to prevent intersecting triangles 241 241 if ((radius > HULLEPSILON) || (helper.Norm() < HULLEPSILON)) { 242 LOG( 4, "DEBUG: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << ".");242 LOG(6, "DEBUG: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << "."); 243 243 return alpha; 244 244 } else { 245 LOG( 3, "DEBUG: NewSphereCenter " << RelativeNewSphereCenter << " is too close to RelativeOldSphereCenter" << RelativeOldSphereCenter << ".");245 LOG(5, "DEBUG: NewSphereCenter " << RelativeNewSphereCenter << " is too close to RelativeOldSphereCenter" << RelativeOldSphereCenter << "."); 246 246 return 2.*M_PI; 247 247 } … … 278 278 } 279 279 280 LOG( 1, "INFO: " << point << " has angle " << phi << " with respect to reference " << reference << ".");280 LOG(4, "INFO: " << point << " has angle " << phi << " with respect to reference " << reference << "."); 281 281 282 282 return phi; … … 545 545 Normal.VectorProduct(OtherBaseline); 546 546 Normal.Normalize(); 547 LOG( 1, "First direction is " << Baseline << ", second direction is " << OtherBaseline << ", normal of intersection plane is " << Normal << ".");547 LOG(3, "First direction is " << Baseline << ", second direction is " << OtherBaseline << ", normal of intersection plane is " << Normal << "."); 548 548 549 549 // project one offset point of OtherBase onto this plane (and add plane offset vector) … … 559 559 *Intersection = line1.getIntersection(line2); 560 560 Normal = (*Intersection) - (Base->endpoints[0]->node->getPosition()); 561 LOG( 1, "Found closest point on " << *Base << " at " << *Intersection << ", factor in line is " << fabs(Normal.ScalarProduct(Baseline)/Baseline.NormSquared()) << ".");561 LOG(3, "Found closest point on " << *Base << " at " << *Intersection << ", factor in line is " << fabs(Normal.ScalarProduct(Baseline)/Baseline.NormSquared()) << "."); 562 562 563 563 return Intersection; … … 751 751 *tecplot << endl; 752 752 // print connectivity 753 LOG( 1, "The following triangles were created:");753 LOG(4, "DEBUG: The following triangles were created:"); 754 754 for (TriangleMap::const_iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++) { 755 LOG( 1, " " << runner->second->endpoints[0]->node->getName() << "<->" << runner->second->endpoints[1]->node->getName() << "<->" << runner->second->endpoints[2]->node->getName());755 LOG(4, " " << runner->second->endpoints[0]->node->getName() << "<->" << runner->second->endpoints[1]->node->getName() << "<->" << runner->second->endpoints[2]->node->getName()); 756 756 *tecplot << LookupList[runner->second->endpoints[0]->node->getNr()] << " " << LookupList[runner->second->endpoints[1]->node->getNr()] << " " << LookupList[runner->second->endpoints[2]->node->getNr()] << endl; 757 757 } … … 778 778 for (PointMap::const_iterator PointRunner = TesselStruct->PointsOnBoundary.begin(); PointRunner != TesselStruct->PointsOnBoundary.end(); PointRunner++) { 779 779 point = PointRunner->second; 780 LOG( 1, "INFO: Current point is " << *point << ".");780 LOG(2, "INFO: Current point is " << *point << "."); 781 781 782 782 // calculate mean concavity over all connected line -
src/Tesselation/triangleintersectionlist.cpp
rdfe054 rb8d215 132 132 // if we have found one, check Scalarproduct between the vector 133 133 Vector TestVector = (Point) - (*(*runner).second); 134 LOG( 1, "INFO: Distance vector between point " << Point134 LOG(4, "DEBUG: Distance vector between point " << Point 135 135 << " and triangle intersection at " << (*(*runner).second) << " is " 136 136 << TestVector << "."); 137 137 if (fabs(TestVector.NormSquared()) < MYEPSILON) {// 138 LOG( 1, "ACCEPT: Point is on the intersected triangle.");138 LOG(3, "ACCEPT: Point is on the intersected triangle."); 139 139 return true; 140 140 } 141 141 const double sign = (*runner).first->NormalVector.ScalarProduct(TestVector); 142 LOG( 1, "INFO: Checking sign of SKP between Distance vector and triangle's NormalVector "142 LOG(4, "DEBUG: Checking sign of SKP between Distance vector and triangle's NormalVector " 143 143 << (*runner).first->NormalVector << ": " << sign << "."); 144 144 if (sign < 0) { 145 LOG( 1, "ACCEPT: Point lies on inner side of intersected triangle.");145 LOG(3, "ACCEPT: Point lies on inner side of intersected triangle."); 146 146 return true; 147 147 } else { 148 LOG( 1, "REJECT: Point lies on outer side of intersected triangle.");148 LOG(3, "REJECT: Point lies on outer side of intersected triangle."); 149 149 return false; 150 150 } … … 224 224 for (DistanceTriangleMap::const_iterator iter = DistanceRange.first; 225 225 iter != DistanceRange.second; ++iter) { 226 LOG( 3, "DEBUG: " << *(iter->second) << " is in the list with " << iter->first << ".");226 LOG(4, "DEBUG: " << *(iter->second) << " is in the list with " << iter->first << "."); 227 227 ++count; 228 228 } 229 LOG( 1, "INFO: There are " << count << " possible triangles at the smallest distance.");229 LOG(3, "INFO: There are " << count << " possible triangles at the smallest distance."); 230 230 } 231 231 // if there is more than one, check all of their normal vectors … … 233 233 // would be truely inside, all of the closest triangles would have to show 234 234 // their inner side 235 LOG( 1, "INFO: Looking for first SKP greater than zero ...");235 LOG(4, "DEBUG: Looking for first SKP greater than zero ..."); 236 236 for (DistanceTriangleMap::const_iterator iter = DistanceRange.first; 237 237 iter != DistanceRange.second; ++iter) { … … 241 241 // construct SKP 242 242 const double sign = (*iter).second->NormalVector.ScalarProduct(TestVector); 243 LOG( 1, "INFO: Checking SKP of closest triangle " << *(iter->second) << " = " << sign << ".");243 LOG(4, "DEBUG: Checking SKP of closest triangle " << *(iter->second) << " = " << sign << "."); 244 244 // if positive return 245 245 if (sign > 0)
Note:
See TracChangeset
for help on using the changeset viewer.