Changeset f67b6e
- Timestamp:
- Nov 25, 2009, 3:29:18 PM (15 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:
- 7273fc
- Parents:
- 8725ed
- Location:
- src
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
src/boundary.cpp
r8725ed rf67b6e 10 10 #include "element.hpp" 11 11 #include "helpers.hpp" 12 #include "info.hpp" 12 13 #include "linkedcell.hpp" 13 14 #include "log.hpp" … … 33 34 double *GetDiametersOfCluster(const Boundaries *BoundaryPtr, const molecule *mol, Tesselation *&TesselStruct, const bool IsAngstroem) 34 35 { 36 Info FunctionInfo(__func__); 35 37 // get points on boundary of NULL was given as parameter 36 38 bool BoundaryFreeFlag = false; … … 53 55 } else { 54 56 BoundaryPoints = BoundaryPtr; 55 Log() << Verbose( 1) << "Using given boundary points set." << endl;57 Log() << Verbose(0) << "Using given boundary points set." << endl; 56 58 } 57 59 // determine biggest "diameter" of cluster for each axis … … 67 69 //Log() << Verbose(1) << "Current component is " << component << ", Othercomponent is " << Othercomponent << "." << endl; 68 70 for (Boundaries::const_iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) { 69 //Log() << Verbose( 2) << "Current runner is " << *(runner->second.second) << "." << endl;71 //Log() << Verbose(1) << "Current runner is " << *(runner->second.second) << "." << endl; 70 72 // seek for the neighbours pair where the Othercomponent sign flips 71 73 Neighbour = runner; … … 82 84 DistanceVector.CopyVector(&runner->second.second->x); 83 85 DistanceVector.SubtractVector(&Neighbour->second.second->x); 84 //Log() << Verbose( 3) << "OldComponent is " << OldComponent << ", new one is " << DistanceVector.x[Othercomponent] << "." << endl;86 //Log() << Verbose(2) << "OldComponent is " << OldComponent << ", new one is " << DistanceVector.x[Othercomponent] << "." << endl; 85 87 } while ((runner != Neighbour) && (fabs(OldComponent / fabs( 86 88 OldComponent) - DistanceVector.x[Othercomponent] / fabs( … … 91 93 OtherNeighbour = BoundaryPoints[axis].end(); 92 94 OtherNeighbour--; 93 //Log() << Verbose( 2) << "The pair, where the sign of OtherComponent flips, is: " << *(Neighbour->second.second) << " and " << *(OtherNeighbour->second.second) << "." << endl;95 //Log() << Verbose(1) << "The pair, where the sign of OtherComponent flips, is: " << *(Neighbour->second.second) << " and " << *(OtherNeighbour->second.second) << "." << endl; 94 96 // now we have found the pair: Neighbour and OtherNeighbour 95 97 OtherVector.CopyVector(&runner->second.second->x); 96 98 OtherVector.SubtractVector(&OtherNeighbour->second.second->x); 97 //Log() << Verbose( 2) << "Distances to Neighbour and OtherNeighbour are " << DistanceVector.x[component] << " and " << OtherVector.x[component] << "." << endl;98 //Log() << Verbose( 2) << "OtherComponents to Neighbour and OtherNeighbour are " << DistanceVector.x[Othercomponent] << " and " << OtherVector.x[Othercomponent] << "." << endl;99 //Log() << Verbose(1) << "Distances to Neighbour and OtherNeighbour are " << DistanceVector.x[component] << " and " << OtherVector.x[component] << "." << endl; 100 //Log() << Verbose(1) << "OtherComponents to Neighbour and OtherNeighbour are " << DistanceVector.x[Othercomponent] << " and " << OtherVector.x[Othercomponent] << "." << endl; 99 101 // do linear interpolation between points (is exact) to extract exact intersection between Neighbour and OtherNeighbour 100 102 w1 = fabs(OtherVector.x[Othercomponent]); … … 103 105 * OtherVector.x[component]) / (w1 + w2)); 104 106 // mark if it has greater diameter 105 //Log() << Verbose( 2) << "Comparing current greatest " << GreatestDiameter[component] << " to new " << tmp << "." << endl;107 //Log() << Verbose(1) << "Comparing current greatest " << GreatestDiameter[component] << " to new " << tmp << "." << endl; 106 108 GreatestDiameter[component] = (GreatestDiameter[component] 107 109 > tmp) ? GreatestDiameter[component] : tmp; 108 110 } //else 109 //Log() << Verbose( 2) << "Saw no sign flip, probably top or bottom node." << endl;111 //Log() << Verbose(1) << "Saw no sign flip, probably top or bottom node." << endl; 110 112 } 111 113 } … … 135 137 Boundaries *GetBoundaryPoints(const molecule *mol, Tesselation *&TesselStruct) 136 138 { 139 Info FunctionInfo(__func__); 137 140 atom *Walker = NULL; 138 141 PointMap PointsOnBoundary; … … 149 152 double angle = 0.; 150 153 151 Log() << Verbose(1) << "Finding all boundary points." << endl;152 154 // 3a. Go through every axis 153 155 for (int axis = 0; axis < NDIM; axis++) { … … 176 178 angle = 0.; // otherwise it's a vector in Axis Direction and unimportant for boundary issues 177 179 178 //Log() << Verbose( 2) << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl;180 //Log() << Verbose(1) << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl; 179 181 if (ProjectedVector.ScalarProduct(&AngleReferenceNormalVector) > 0) { 180 182 angle = 2. * M_PI - angle; 181 183 } 182 Log() << Verbose( 2) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): " << ProjectedVector << endl;184 Log() << Verbose(1) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): " << ProjectedVector << endl; 183 185 BoundaryTestPair = BoundaryPoints[axis].insert(BoundariesPair(angle, DistancePair (radius, Walker))); 184 186 if (!BoundaryTestPair.second) { // same point exists, check first r, then distance of original vectors to center of gravity … … 210 212 // printing all inserted for debugging 211 213 // { 212 // Log() << Verbose( 2) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;214 // Log() << Verbose(1) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl; 213 215 // int i=0; 214 216 // for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) { 215 217 // if (runner != BoundaryPoints[axis].begin()) 216 // Log() << Verbose( 2) << ", " << i << ": " << *runner->second.second;218 // Log() << Verbose(0) << ", " << i << ": " << *runner->second.second; 217 219 // else 218 // Log() << Verbose( 2) << i << ": " << *runner->second.second;220 // Log() << Verbose(0) << i << ": " << *runner->second.second; 219 221 // i++; 220 222 // } 221 // Log() << Verbose( 2) << endl;223 // Log() << Verbose(0) << endl; 222 224 // } 223 225 // 3c. throw out points whose distance is less than the mean of left and right neighbours … … 249 251 SideA.SubtractVector(MolCenter); 250 252 SideA.ProjectOntoPlane(&AxisVector); 251 // Log() << Verbose( 0) << "SideA: " << SideA << endl;253 // Log() << Verbose(1) << "SideA: " << SideA << endl; 252 254 253 255 SideB.CopyVector(&right->second.second->x); 254 256 SideB.SubtractVector(MolCenter); 255 257 SideB.ProjectOntoPlane(&AxisVector); 256 // Log() << Verbose( 0) << "SideB: " << SideB << endl;258 // Log() << Verbose(1) << "SideB: " << SideB << endl; 257 259 258 260 SideC.CopyVector(&left->second.second->x); 259 261 SideC.SubtractVector(&right->second.second->x); 260 262 SideC.ProjectOntoPlane(&AxisVector); 261 // Log() << Verbose( 0) << "SideC: " << SideC << endl;263 // Log() << Verbose(1) << "SideC: " << SideC << endl; 262 264 263 265 SideH.CopyVector(&runner->second.second->x); 264 266 SideH.SubtractVector(MolCenter); 265 267 SideH.ProjectOntoPlane(&AxisVector); 266 // Log() << Verbose( 0) << "SideH: " << SideH << endl;268 // Log() << Verbose(1) << "SideH: " << SideH << endl; 267 269 268 270 // calculate each length … … 277 279 const double delta = SideC.Angle(&SideH); 278 280 const double MinDistance = a * sin(beta) / (sin(delta)) * (((alpha < M_PI / 2.) || (gamma < M_PI / 2.)) ? 1. : -1.); 279 //Log() << Verbose( 2) << " I calculated: a = " << a << ", h = " << h << ", beta(" << left->second.second->Name << "," << left->second.second->Name << "-" << right->second.second->Name << ") = " << beta << ", delta(" << left->second.second->Name << "," << runner->second.second->Name << ") = " << delta << ", Min = " << MinDistance << "." << endl;281 //Log() << Verbose(1) << " I calculated: a = " << a << ", h = " << h << ", beta(" << left->second.second->Name << "," << left->second.second->Name << "-" << right->second.second->Name << ") = " << beta << ", delta(" << left->second.second->Name << "," << runner->second.second->Name << ") = " << delta << ", Min = " << MinDistance << "." << endl; 280 282 Log() << Verbose(1) << "Checking CoG distance of runner " << *runner->second.second << " " << h << " against triangle's side length spanned by (" << *left->second.second << "," << *right->second.second << ") of " << MinDistance << "." << endl; 281 283 if ((fabs(h / fabs(h) - MinDistance / fabs(MinDistance)) < MYEPSILON) && ((h - MinDistance)) < -MYEPSILON) { … … 303 305 void FindConvexBorder(const molecule* mol, Tesselation *&TesselStruct, const LinkedCell *LCList, const char *filename) 304 306 { 307 Info FunctionInfo(__func__); 305 308 bool BoundaryFreeFlag = false; 306 309 Boundaries *BoundaryPoints = NULL; 307 308 Log() << Verbose(1) << "Begin of FindConvexBorder" << endl;309 310 310 311 if (TesselStruct != NULL) // free if allocated … … 317 318 BoundaryPoints = GetBoundaryPoints(mol, TesselStruct); 318 319 } else { 319 Log() << Verbose( 1) << "Using given boundary points set." << endl;320 Log() << Verbose(0) << "Using given boundary points set." << endl; 320 321 } 321 322 … … 323 324 for (int axis=0; axis < NDIM; axis++) 324 325 { 325 Log() << Verbose( 2) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;326 Log() << Verbose(1) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl; 326 327 int i=0; 327 328 for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) { 328 329 if (runner != BoundaryPoints[axis].begin()) 329 Log() << Verbose( 2) << ", " << i << ": " << *runner->second.second;330 Log() << Verbose(0) << ", " << i << ": " << *runner->second.second; 330 331 else 331 Log() << Verbose( 2) << i << ": " << *runner->second.second;332 Log() << Verbose(0) << i << ": " << *runner->second.second; 332 333 i++; 333 334 } 334 Log() << Verbose( 2) << endl;335 Log() << Verbose(0) << endl; 335 336 } 336 337 … … 341 342 eLog() << Verbose(2) << "Point " << *(runner->second.second) << " is already present!" << endl; 342 343 343 Log() << Verbose( 2) << "I found " << TesselStruct->PointsOnBoundaryCount << " points on the convex boundary." << endl;344 Log() << Verbose(0) << "I found " << TesselStruct->PointsOnBoundaryCount << " points on the convex boundary." << endl; 344 345 // now we have the whole set of edge points in the BoundaryList 345 346 … … 347 348 // Log() << Verbose(1) << "Listing PointsOnBoundary:"; 348 349 // for(PointMap::iterator runner = PointsOnBoundary.begin(); runner != PointsOnBoundary.end(); runner++) { 349 // Log() << Verbose( 1) << " " << *runner->second;350 // Log() << Verbose(0) << " " << *runner->second; 350 351 // } 351 // Log() << Verbose( 1) << endl;352 // Log() << Verbose(0) << endl; 352 353 353 354 // 3a. guess starting triangle … … 359 360 // 3c. check whether all atoms lay inside the boundary, if not, add to boundary points, segment triangle into three with the new point 360 361 if (!TesselStruct->InsertStraddlingPoints(mol, LCList)) 361 Log() << Verbose(1) << "Insertion of straddling points failed!" << endl;362 363 Log() << Verbose( 2) << "I created " << TesselStruct->TrianglesOnBoundary.size() << " intermediate triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points." << endl;362 eLog() << Verbose(1) << "Insertion of straddling points failed!" << endl; 363 364 Log() << Verbose(0) << "I created " << TesselStruct->TrianglesOnBoundary.size() << " intermediate triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points." << endl; 364 365 365 366 // 4. Store triangles in tecplot file … … 411 412 // Log() << Verbose(1) << "Correction of concave tesselpoints failed!" << endl; 412 413 413 Log() << Verbose( 2) << "I created " << TesselStruct->TrianglesOnBoundary.size() << " triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points." << endl;414 Log() << Verbose(0) << "I created " << TesselStruct->TrianglesOnBoundary.size() << " triangles with " << TesselStruct->LinesOnBoundary.size() << " lines and " << TesselStruct->PointsOnBoundary.size() << " points." << endl; 414 415 415 416 // 4. Store triangles in tecplot file … … 437 438 if (BoundaryFreeFlag) 438 439 delete[] (BoundaryPoints); 439 440 Log() << Verbose(1) << "End of FindConvexBorder" << endl;441 440 }; 442 441 … … 450 449 bool RemoveAllBoundaryPoints(class Tesselation *&TesselStruct, const molecule * const mol, const char * const filename) 451 450 { 451 Info FunctionInfo(__func__); 452 452 int i=0; 453 453 char number[MAXSTRINGSIZE]; … … 460 460 PointMap::iterator PointRunner; 461 461 while (!TesselStruct->PointsOnBoundary.empty()) { 462 Log() << Verbose( 2) << "Remaining points are: ";462 Log() << Verbose(1) << "Remaining points are: "; 463 463 for (PointMap::iterator PointSprinter = TesselStruct->PointsOnBoundary.begin(); PointSprinter != TesselStruct->PointsOnBoundary.end(); PointSprinter++) 464 Log() << Verbose( 2) << *(PointSprinter->second) << "\t";465 Log() << Verbose( 2) << endl;464 Log() << Verbose(0) << *(PointSprinter->second) << "\t"; 465 Log() << Verbose(0) << endl; 466 466 467 467 PointRunner = TesselStruct->PointsOnBoundary.begin(); … … 503 503 double ConvexizeNonconvexEnvelope(class Tesselation *&TesselStruct, const molecule * const mol, const char * const filename) 504 504 { 505 Info FunctionInfo(__func__); 505 506 double volume = 0; 506 507 class BoundaryPointSet *point = NULL; … … 516 517 int run = 0; 517 518 518 Log() << Verbose(0) << "Begin of ConvexizeNonconvexEnvelope" << endl;519 520 519 // check whether there is something to work on 521 520 if (TesselStruct == NULL) { … … 539 538 for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) { 540 539 line = LineRunner->second; 541 Log() << Verbose( 2) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;540 Log() << Verbose(1) << "INFO: Current line of point " << *point << " is " << *line << "." << endl; 542 541 if (!line->CheckConvexityCriterion()) { 543 542 // remove the point if needed … … 604 603 605 604 // end 606 Log() << Verbose(1) << "Volume is " << volume << "." << endl; 607 Log() << Verbose(0) << "End of ConvexizeNonconvexEnvelope" << endl; 605 Log() << Verbose(0) << "Volume is " << volume << "." << endl; 608 606 return volume; 609 607 }; … … 619 617 double VolumeOfConvexEnvelope(class Tesselation *TesselStruct, class config *configuration) 620 618 { 619 Info FunctionInfo(__func__); 621 620 bool IsAngstroem = configuration->GetIsAngstroem(); 622 621 double volume = 0.; … … 625 624 626 625 // 6a. Every triangle forms a pyramid with the center of gravity as its peak, sum up the volumes 627 Log() << Verbose(1)628 << "Calculating the volume of the pyramids formed out of triangles and center of gravity."629 << endl;630 626 for (TriangleMap::iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++) 631 627 { // go through every triangle, calculate volume of its pyramid with CoG as peak … … 642 638 const double h = x.Norm(); // distance of CoG to triangle 643 639 const double PyramidVolume = (1. / 3.) * G * h; // this formula holds for _all_ pyramids (independent of n-edge base or (not) centered peak) 644 Log() << Verbose( 2) << "Area of triangle is " << setprecision(10) << G << " "640 Log() << Verbose(1) << "Area of triangle is " << setprecision(10) << G << " " 645 641 << (IsAngstroem ? "angstrom" : "atomiclength") << "^2, height is " 646 642 << h << " and the volume is " << PyramidVolume << " " … … 664 660 void StoreTrianglesinFile(const molecule * const mol, const Tesselation *&TesselStruct, const char *filename, const char *extraSuffix) 665 661 { 662 Info FunctionInfo(__func__); 666 663 // 4. Store triangles in tecplot file 667 664 if (filename != NULL) { … … 698 695 void PrepareClustersinWater(config *configuration, molecule *mol, double ClusterVolume, double celldensity) 699 696 { 697 Info FunctionInfo(__func__); 700 698 bool IsAngstroem = true; 701 699 double *GreatestDiameter = NULL; … … 798 796 molecule * FillBoxWithMolecule(MoleculeListClass *List, molecule *filler, config &configuration, double distance[NDIM], double RandomAtomDisplacement, double RandomMolDisplacement, bool DoRandomRotation) 799 797 { 798 Info FunctionInfo(__func__); 800 799 molecule *Filling = new molecule(filler->elemente); 801 800 Vector CurrentPosition; … … 814 813 double phi[NDIM]; 815 814 class Tesselation *TesselStruct[List->ListOfMolecules.size()]; 816 817 Log() << Verbose(0) << "Begin of FillBoxWithMolecule" << endl;818 815 819 816 i=0; … … 877 874 for (int i=0;i<NDIM;i++) 878 875 FillerTranslations.x[i] = RandomMolDisplacement*(rand()/(RAND_MAX/2.) - 1.); 879 Log() << Verbose( 3) << "INFO: Translating this filler by " << FillerTranslations << "." << endl;876 Log() << Verbose(2) << "INFO: Translating this filler by " << FillerTranslations << "." << endl; 880 877 881 878 // go through all atoms … … 938 935 delete(TesselStruct[i]); 939 936 } 940 Log() << Verbose(0) << "End of FillBoxWithMolecule" << endl;941 942 937 return Filling; 943 938 }; … … 954 949 void FindNonConvexBorder(const molecule* const mol, Tesselation *&TesselStruct, const LinkedCell *&LCList, const double RADIUS, const char *filename = NULL) 955 950 { 951 Info FunctionInfo(__func__); 956 952 bool freeLC = false; 957 953 CandidateForTesselation *baseline; … … 959 955 bool OneLoopWithoutSuccessFlag = true; // marks whether we went once through all baselines without finding any without two triangles 960 956 bool TesselationFailFlag = false; 961 962 Log() << Verbose(1) << "Entering search for non convex hull. " << endl; 957 BoundaryTriangleSet *T = NULL; 958 963 959 if (TesselStruct == NULL) { 964 960 Log() << Verbose(1) << "Allocating Tesselation struct ..." << endl; … … 970 966 } 971 967 972 Log() << Verbose(0) << "Begin of FindNonConvexBorder\n";973 974 968 // initialise Linked Cell 975 969 if (LCList == NULL) { … … 984 978 while ((!TesselStruct->OpenLines.empty()) && (OneLoopWithoutSuccessFlag)) { 985 979 // 2a. fill all new OpenLines 986 Log() << Verbose( 0) << "There are " << TesselStruct->OpenLines.size() << " open lines to scan for candidates:" << endl;980 Log() << Verbose(1) << "There are " << TesselStruct->OpenLines.size() << " open lines to scan for candidates:" << endl; 987 981 for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++) 988 Log() << Verbose( 1) << *(Runner->second) << endl;982 Log() << Verbose(2) << *(Runner->second) << endl; 989 983 990 984 for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++) { 991 985 baseline = Runner->second; 992 if (baseline->point == NULL) { 993 Log() << Verbose(0) << "Finding best candidate for open line " << *baseline->BaseLine << endl; 994 TesselationFailFlag = TesselStruct->FindNextSuitableTriangle(*baseline, *(((baseline->BaseLine->triangles.begin()))->second), RADIUS, LCList); //the line is there, so there is a triangle, but only one. 986 if (baseline->pointlist.empty()) { 987 T = (((baseline->BaseLine->triangles.begin()))->second); 988 Log() << Verbose(1) << "Finding best candidate for open line " << *baseline->BaseLine << " of triangle " << *T << endl; 989 TesselationFailFlag = TesselStruct->FindNextSuitableTriangle(*baseline, *T, RADIUS, LCList); //the line is there, so there is a triangle, but only one. 995 990 } 996 991 } … … 998 993 // 2b. search for smallest ShortestAngle among all candidates 999 994 double ShortestAngle = 4.*M_PI; 1000 Log() << Verbose( 0) << "There are " << TesselStruct->OpenLines.size() << " open lines to scan for the best candidates:" << endl;995 Log() << Verbose(1) << "There are " << TesselStruct->OpenLines.size() << " open lines to scan for the best candidates:" << endl; 1001 996 for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++) 1002 Log() << Verbose( 1) << *(Runner->second) << endl;997 Log() << Verbose(2) << *(Runner->second) << endl; 1003 998 1004 999 for (CandidateMap::iterator Runner = TesselStruct->OpenLines.begin(); Runner != TesselStruct->OpenLines.end(); Runner++) { … … 1006 1001 baseline = Runner->second; 1007 1002 ShortestAngle = baseline->ShortestAngle; 1008 //Log() << Verbose( 0) << "New best candidate is " << *baseline->BaseLine << " with point " << *baseline->point << " and angle " << baseline->ShortestAngle << endl;1003 //Log() << Verbose(1) << "New best candidate is " << *baseline->BaseLine << " with point " << *baseline->point << " and angle " << baseline->ShortestAngle << endl; 1009 1004 } 1010 1005 } 1011 if ((ShortestAngle == 4.*M_PI) || (baseline->point == NULL))1006 if ((ShortestAngle == 4.*M_PI) || (baseline->pointlist.empty())) 1012 1007 OneLoopWithoutSuccessFlag = false; 1013 1008 else { … … 1041 1036 // } 1042 1037 1043 // Purges surplus triangles.1044 TesselStruct->RemoveDegeneratedTriangles();1045 1046 // check envelope for consistency1047 CheckListOfBaselines(TesselStruct);1038 // // Purges surplus triangles. 1039 // TesselStruct->RemoveDegeneratedTriangles(); 1040 // 1041 // // check envelope for consistency 1042 // CheckListOfBaselines(TesselStruct); 1048 1043 1049 1044 // write final envelope … … 1053 1048 if (freeLC) 1054 1049 delete(LCList); 1055 Log() << Verbose(0) << "End of FindNonConvexBorder\n";1056 1050 }; 1057 1051 … … 1065 1059 Vector* FindEmbeddingHole(MoleculeListClass *mols, molecule *srcmol) 1066 1060 { 1061 Info FunctionInfo(__func__); 1067 1062 Vector *Center = new Vector; 1068 1063 Center->Zero(); -
src/builder.cpp
r8725ed rf67b6e 1577 1577 molecules->ListOfMolecules.remove(mol); 1578 1578 molecules->DissectMoleculeIntoConnectedSubgraphs(mol,&configuration); 1579 delete(mol); 1579 1580 if (molecules->ListOfMolecules.size() != 0) { 1580 1581 for (MoleculeList::iterator ListRunner = molecules->ListOfMolecules.begin(); ListRunner != molecules->ListOfMolecules.end(); ListRunner++) … … 1807 1808 Log() << Verbose(0) << "Clocks for this operation: " << (end-start) << ", time: " << ((double)(end-start)/CLOCKS_PER_SEC) << "s." << endl; 1808 1809 delete(LCList); 1810 delete(T); 1809 1811 argptr+=2; 1810 1812 } -
src/tesselation.cpp
r8725ed rf67b6e 9 9 10 10 #include "helpers.hpp" 11 #include "info.hpp" 11 12 #include "linkedcell.hpp" 12 13 #include "log.hpp" … … 27 28 Nr(-1) 28 29 { 30 Info FunctionInfo(__func__); 31 Log() << Verbose(1) << "Adding noname." << endl; 29 32 }; 30 33 … … 32 35 * \param *Walker TesselPoint this boundary point represents 33 36 */ 34 BoundaryPointSet::BoundaryPointSet(TesselPoint * Walker) 35 { 36 node = Walker; 37 LinesCount = 0; 38 Nr = Walker->nr; 39 value = 0.; 37 BoundaryPointSet::BoundaryPointSet(TesselPoint * Walker) : 38 LinesCount(0), 39 node(Walker), 40 value(0.), 41 Nr(Walker->nr) 42 { 43 Info FunctionInfo(__func__); 44 Log() << Verbose(1) << "Adding at " << *(node->node) << endl; 40 45 }; 41 46 … … 46 51 BoundaryPointSet::~BoundaryPointSet() 47 52 { 48 //Log() << Verbose(5) << "Erasing point nr. " << Nr << "." << endl; 53 Info FunctionInfo(__func__); 54 //Log() << Verbose(0) << "Erasing point nr. " << Nr << "." << endl; 49 55 if (!lines.empty()) 50 56 eLog() << Verbose(2) << "Memory Leak! I " << *this << " am still connected to some lines." << endl; … … 57 63 void BoundaryPointSet::AddLine(class BoundaryLineSet *line) 58 64 { 59 Log() << Verbose(6) << "Adding " << *this << " to line " << *line << "." 65 Info FunctionInfo(__func__); 66 Log() << Verbose(1) << "Adding " << *this << " to line " << *line << "." 60 67 << endl; 61 68 if (line->endpoints[0] == this) … … 88 95 Nr(-1) 89 96 { 97 Info FunctionInfo(__func__); 90 98 for (int i = 0; i < 2; i++) 91 99 endpoints[i] = NULL; … … 99 107 BoundaryLineSet::BoundaryLineSet(class BoundaryPointSet *Point[2], const int number) 100 108 { 109 Info FunctionInfo(__func__); 101 110 // set number 102 111 Nr = number; … … 109 118 skipped = false; 110 119 // clear triangles list 111 Log() << Verbose( 5) << "New Line with endpoints " << *this << "." << endl;120 Log() << Verbose(0) << "New Line with endpoints " << *this << "." << endl; 112 121 }; 113 122 … … 118 127 BoundaryLineSet::~BoundaryLineSet() 119 128 { 129 Info FunctionInfo(__func__); 120 130 int Numbers[2]; 121 131 … … 136 146 for (LineMap::iterator Runner = erasor.first; Runner != erasor.second; Runner++) 137 147 if ((*Runner).second == this) { 138 //Log() << Verbose( 5) << "Removing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;148 //Log() << Verbose(0) << "Removing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl; 139 149 endpoints[i]->lines.erase(Runner); 140 150 break; … … 142 152 } else { // there's just a single line left 143 153 if (endpoints[i]->lines.erase(Nr)) { 144 //Log() << Verbose( 5) << "Removing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;154 //Log() << Verbose(0) << "Removing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl; 145 155 } 146 156 } 147 157 if (endpoints[i]->lines.empty()) { 148 //Log() << Verbose( 5) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl;158 //Log() << Verbose(0) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl; 149 159 if (endpoints[i] != NULL) { 150 160 delete(endpoints[i]); … … 163 173 void BoundaryLineSet::AddTriangle(class BoundaryTriangleSet *triangle) 164 174 { 165 Log() << Verbose(6) << "Add " << triangle->Nr << " to line " << *this << "." << endl; 175 Info FunctionInfo(__func__); 176 Log() << Verbose(0) << "Add " << triangle->Nr << " to line " << *this << "." << endl; 166 177 triangles.insert(TrianglePair(triangle->Nr, triangle)); 167 178 }; … … 173 184 bool BoundaryLineSet::IsConnectedTo(class BoundaryLineSet *line) 174 185 { 186 Info FunctionInfo(__func__); 175 187 if ((endpoints[0] == line->endpoints[0]) || (endpoints[1] == line->endpoints[0]) || (endpoints[0] == line->endpoints[1]) || (endpoints[1] == line->endpoints[1])) 176 188 return true; … … 187 199 bool BoundaryLineSet::CheckConvexityCriterion() 188 200 { 201 Info FunctionInfo(__func__); 189 202 Vector BaseLineCenter, BaseLineNormal, BaseLine, helper[2], NormalCheck; 190 203 // get the two triangles 191 204 if (triangles.size() != 2) { 192 eLog() << Verbose( 1) << "Baseline " << *this << " is connected to less than two triangles, Tesselation incomplete!" << endl;205 eLog() << Verbose(0) << "Baseline " << *this << " is connected to less than two triangles, Tesselation incomplete!" << endl; 193 206 return true; 194 207 } 195 208 // check normal vectors 196 209 // have a normal vector on the base line pointing outwards 197 //Log() << Verbose( 3) << "INFO: " << *this << " has vectors at " << *(endpoints[0]->node->node) << " and at " << *(endpoints[1]->node->node) << "." << endl;210 //Log() << Verbose(0) << "INFO: " << *this << " has vectors at " << *(endpoints[0]->node->node) << " and at " << *(endpoints[1]->node->node) << "." << endl; 198 211 BaseLineCenter.CopyVector(endpoints[0]->node->node); 199 212 BaseLineCenter.AddVector(endpoints[1]->node->node); … … 201 214 BaseLine.CopyVector(endpoints[0]->node->node); 202 215 BaseLine.SubtractVector(endpoints[1]->node->node); 203 //Log() << Verbose( 3) << "INFO: Baseline is " << BaseLine << " and its center is at " << BaseLineCenter << "." << endl;216 //Log() << Verbose(0) << "INFO: Baseline is " << BaseLine << " and its center is at " << BaseLineCenter << "." << endl; 204 217 205 218 BaseLineNormal.Zero(); … … 209 222 class BoundaryPointSet *node = NULL; 210 223 for(TriangleMap::iterator runner = triangles.begin(); runner != triangles.end(); runner++) { 211 //Log() << Verbose( 3) << "INFO: NormalVector of " << *(runner->second) << " is " << runner->second->NormalVector << "." << endl;224 //Log() << Verbose(0) << "INFO: NormalVector of " << *(runner->second) << " is " << runner->second->NormalVector << "." << endl; 212 225 NormalCheck.AddVector(&runner->second->NormalVector); 213 226 NormalCheck.Scale(sign); … … 216 229 BaseLineNormal.CopyVector(&runner->second->NormalVector); // yes, copy second on top of first 217 230 else { 218 Log() << Verbose(1) << "CRITICAL: Triangle " << *runner->second << " has zero normal vector!" << endl; 219 exit(255); 231 eLog() << Verbose(0) << "Triangle " << *runner->second << " has zero normal vector!" << endl; 220 232 } 221 233 node = runner->second->GetThirdEndpoint(this); 222 234 if (node != NULL) { 223 //Log() << Verbose( 3) << "INFO: Third node for triangle " << *(runner->second) << " is " << *node << " at " << *(node->node->node) << "." << endl;235 //Log() << Verbose(0) << "INFO: Third node for triangle " << *(runner->second) << " is " << *node << " at " << *(node->node->node) << "." << endl; 224 236 helper[i].CopyVector(node->node->node); 225 237 helper[i].SubtractVector(&BaseLineCenter); 226 238 helper[i].MakeNormalVector(&BaseLine); // we want to compare the triangle's heights' angles! 227 //Log() << Verbose( 4) << "INFO: Height vector with respect to baseline is " << helper[i] << "." << endl;239 //Log() << Verbose(0) << "INFO: Height vector with respect to baseline is " << helper[i] << "." << endl; 228 240 i++; 229 241 } else { 230 //eLog() << Verbose(1) << "I cannot find third node in triangle, something's wrong." << endl;242 eLog() << Verbose(1) << "I cannot find third node in triangle, something's wrong." << endl; 231 243 return true; 232 244 } 233 245 } 234 //Log() << Verbose( 3) << "INFO: BaselineNormal is " << BaseLineNormal << "." << endl;246 //Log() << Verbose(0) << "INFO: BaselineNormal is " << BaseLineNormal << "." << endl; 235 247 if (NormalCheck.NormSquared() < MYEPSILON) { 236 Log() << Verbose( 3) << "ACCEPT: Normalvectors of both triangles are the same: convex." << endl;248 Log() << Verbose(0) << "ACCEPT: Normalvectors of both triangles are the same: convex." << endl; 237 249 return true; 238 250 } … … 240 252 double angle = GetAngle(helper[0], helper[1], BaseLineNormal); 241 253 if ((angle - M_PI) > -MYEPSILON) { 242 Log() << Verbose( 3) << "ACCEPT: Angle is greater than pi: convex." << endl;254 Log() << Verbose(0) << "ACCEPT: Angle is greater than pi: convex." << endl; 243 255 return true; 244 256 } else { 245 Log() << Verbose( 3) << "REJECT: Angle is less than pi: concave." << endl;257 Log() << Verbose(0) << "REJECT: Angle is less than pi: concave." << endl; 246 258 return false; 247 259 } … … 254 266 bool BoundaryLineSet::ContainsBoundaryPoint(class BoundaryPointSet *point) 255 267 { 268 Info FunctionInfo(__func__); 256 269 for(int i=0;i<2;i++) 257 270 if (point == endpoints[i]) … … 266 279 class BoundaryPointSet *BoundaryLineSet::GetOtherEndpoint(class BoundaryPointSet *point) 267 280 { 281 Info FunctionInfo(__func__); 268 282 if (endpoints[0] == point) 269 283 return endpoints[1]; … … 291 305 Nr(-1) 292 306 { 307 Info FunctionInfo(__func__); 293 308 for (int i = 0; i < 3; i++) 294 309 { … … 305 320 Nr(number) 306 321 { 322 Info FunctionInfo(__func__); 307 323 // set number 308 324 // set lines 309 Log() << Verbose(5) << "New triangle " << Nr << ":"; 310 for (int i = 0; i < 3; i++) 311 { 312 lines[i] = line[i]; 313 lines[i]->AddTriangle(this); 314 } 325 for (int i = 0; i < 3; i++) { 326 lines[i] = line[i]; 327 lines[i]->AddTriangle(this); 328 } 315 329 // get ascending order of endpoints 316 map<int, class BoundaryPointSet *>OrderMap;330 PointMap OrderMap; 317 331 for (int i = 0; i < 3; i++) 318 332 // for all three lines 319 for (int j = 0; j < 2; j++) 320 { // for both endpoints 321 OrderMap.insert(pair<int, class BoundaryPointSet *> ( 322 line[i]->endpoints[j]->Nr, line[i]->endpoints[j])); 323 // and we don't care whether insertion fails 324 } 333 for (int j = 0; j < 2; j++) { // for both endpoints 334 OrderMap.insert(pair<int, class BoundaryPointSet *> ( 335 line[i]->endpoints[j]->Nr, line[i]->endpoints[j])); 336 // and we don't care whether insertion fails 337 } 325 338 // set endpoints 326 339 int Counter = 0; 327 Log() << Verbose(6) << " with end points "; 328 for (map<int, class BoundaryPointSet *>::iterator runner = OrderMap.begin(); runner 329 != OrderMap.end(); runner++) 330 { 331 endpoints[Counter] = runner->second; 332 Log() << Verbose(0) << " " << *endpoints[Counter]; 333 Counter++; 334 } 335 if (Counter < 3) 336 { 337 eLog() << Verbose(0) << "ERROR! We have a triangle with only two distinct endpoints!" << endl; 338 performCriticalExit(); 339 } 340 Log() << Verbose(0) << "." << endl; 340 Log() << Verbose(0) << "New triangle " << Nr << " with end points: " << endl; 341 for (PointMap::iterator runner = OrderMap.begin(); runner != OrderMap.end(); runner++) { 342 endpoints[Counter] = runner->second; 343 Log() << Verbose(0) << " " << *endpoints[Counter] << endl; 344 Counter++; 345 } 346 if (Counter < 3) { 347 eLog() << Verbose(0) << "We have a triangle with only two distinct endpoints!" << endl; 348 performCriticalExit(); 349 } 341 350 }; 342 351 … … 347 356 BoundaryTriangleSet::~BoundaryTriangleSet() 348 357 { 358 Info FunctionInfo(__func__); 349 359 for (int i = 0; i < 3; i++) { 350 360 if (lines[i] != NULL) { 351 361 if (lines[i]->triangles.erase(Nr)) { 352 //Log() << Verbose( 5) << "Triangle Nr." << Nr << " erased in line " << *lines[i] << "." << endl;362 //Log() << Verbose(0) << "Triangle Nr." << Nr << " erased in line " << *lines[i] << "." << endl; 353 363 } 354 364 if (lines[i]->triangles.empty()) { 355 //Log() << Verbose( 5) << *lines[i] << " is no more attached to any triangle, erasing." << endl;365 //Log() << Verbose(0) << *lines[i] << " is no more attached to any triangle, erasing." << endl; 356 366 delete (lines[i]); 357 367 lines[i] = NULL; … … 359 369 } 360 370 } 361 //Log() << Verbose( 5) << "Erasing triangle Nr." << Nr << " itself." << endl;371 //Log() << Verbose(0) << "Erasing triangle Nr." << Nr << " itself." << endl; 362 372 }; 363 373 … … 368 378 void BoundaryTriangleSet::GetNormalVector(Vector &OtherVector) 369 379 { 380 Info FunctionInfo(__func__); 370 381 // get normal vector 371 382 NormalVector.MakeNormalVector(endpoints[0]->node->node, endpoints[1]->node->node, endpoints[2]->node->node); … … 374 385 if (NormalVector.ScalarProduct(&OtherVector) > 0.) 375 386 NormalVector.Scale(-1.); 387 Log() << Verbose(1) << "Normal Vector is " << NormalVector << "." << endl; 376 388 }; 377 389 … … 390 402 bool BoundaryTriangleSet::GetIntersectionInsideTriangle(Vector *MolCenter, Vector *x, Vector *Intersection) 391 403 { 404 Info FunctionInfo(__func__); 392 405 Vector CrossPoint; 393 406 Vector helper; 394 407 395 408 if (!Intersection->GetIntersectionWithPlane(&NormalVector, endpoints[0]->node->node, MolCenter, x)) { 396 Log() << Verbose(1) << "Alas! Intersection with plane failed - at least numerically - the intersection is not on the plane!" << endl;409 eLog() << Verbose(1) << "Alas! Intersection with plane failed - at least numerically - the intersection is not on the plane!" << endl; 397 410 return false; 398 411 } … … 410 423 } while (CrossPoint.NormSquared() < MYEPSILON); 411 424 if (i==3) { 412 eLog() << Verbose(1) << "Could not find any cross points, something's utterly wrong here!" << endl; 413 exit(255); 425 eLog() << Verbose(0) << "Could not find any cross points, something's utterly wrong here!" << endl; 414 426 } 415 427 CrossPoint.SubtractVector(endpoints[i%3]->node->node); // cross point was returned as absolute vector … … 430 442 bool BoundaryTriangleSet::ContainsBoundaryLine(class BoundaryLineSet *line) 431 443 { 444 Info FunctionInfo(__func__); 432 445 for(int i=0;i<3;i++) 433 446 if (line == lines[i]) … … 442 455 bool BoundaryTriangleSet::ContainsBoundaryPoint(class BoundaryPointSet *point) 443 456 { 457 Info FunctionInfo(__func__); 444 458 for(int i=0;i<3;i++) 445 459 if (point == endpoints[i]) … … 454 468 bool BoundaryTriangleSet::ContainsBoundaryPoint(class TesselPoint *point) 455 469 { 470 Info FunctionInfo(__func__); 456 471 for(int i=0;i<3;i++) 457 472 if (point == endpoints[i]->node) … … 466 481 bool BoundaryTriangleSet::IsPresentTupel(class BoundaryPointSet *Points[3]) 467 482 { 483 Info FunctionInfo(__func__); 468 484 return (((endpoints[0] == Points[0]) 469 485 || (endpoints[0] == Points[1]) … … 487 503 bool BoundaryTriangleSet::IsPresentTupel(class BoundaryTriangleSet *T) 488 504 { 505 Info FunctionInfo(__func__); 489 506 return (((endpoints[0] == T->endpoints[0]) 490 507 || (endpoints[0] == T->endpoints[1]) … … 508 525 class BoundaryPointSet *BoundaryTriangleSet::GetThirdEndpoint(class BoundaryLineSet *line) 509 526 { 527 Info FunctionInfo(__func__); 510 528 // sanity check 511 529 if (!ContainsBoundaryLine(line)) … … 524 542 void BoundaryTriangleSet::GetCenter(Vector *center) 525 543 { 544 Info FunctionInfo(__func__); 526 545 center->Zero(); 527 546 for(int i=0;i<3;i++) … … 536 555 ostream &operator <<(ostream &ost, const BoundaryTriangleSet &a) 537 556 { 538 ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << " at " << *a.endpoints[0]->node->node << "," 539 << a.endpoints[1]->node->Name << " at " << *a.endpoints[1]->node->node << "," << a.endpoints[2]->node->Name << " at " << *a.endpoints[2]->node->node << "]"; 557 ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << "," << a.endpoints[1]->node->Name << "," << a.endpoints[2]->node->Name << "]"; 558 // ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << " at " << *a.endpoints[0]->node->node << "," 559 // << a.endpoints[1]->node->Name << " at " << *a.endpoints[1]->node->node << "," << a.endpoints[2]->node->Name << " at " << *a.endpoints[2]->node->node << "]"; 540 560 return ost; 541 561 }; … … 547 567 TesselPoint::TesselPoint() 548 568 { 569 Info FunctionInfo(__func__); 549 570 node = NULL; 550 571 nr = -1; … … 556 577 TesselPoint::~TesselPoint() 557 578 { 579 Info FunctionInfo(__func__); 558 580 }; 559 581 … … 570 592 ostream & TesselPoint::operator << (ostream &ost) 571 593 { 594 Info FunctionInfo(__func__); 572 595 ost << "[" << (Name) << "|" << this << "]"; 573 596 return ost; … … 581 604 PointCloud::PointCloud() 582 605 { 606 Info FunctionInfo(__func__); 583 607 }; 584 608 … … 587 611 PointCloud::~PointCloud() 588 612 { 613 Info FunctionInfo(__func__); 589 614 }; 590 615 … … 594 619 */ 595 620 CandidateForTesselation::CandidateForTesselation (BoundaryLineSet* line) : 596 point(NULL),597 621 BaseLine(line), 598 622 ShortestAngle(2.*M_PI), 599 623 OtherShortestAngle(2.*M_PI) 600 624 { 625 Info FunctionInfo(__func__); 601 626 }; 602 627 … … 605 630 */ 606 631 CandidateForTesselation::CandidateForTesselation (TesselPoint *candidate, BoundaryLineSet* line, Vector OptCandidateCenter, Vector OtherOptCandidateCenter) : 607 point(candidate),608 632 BaseLine(line), 609 633 ShortestAngle(2.*M_PI), 610 634 OtherShortestAngle(2.*M_PI) 611 635 { 636 Info FunctionInfo(__func__); 612 637 OptCenter.CopyVector(&OptCandidateCenter); 613 638 OtherOptCenter.CopyVector(&OtherOptCandidateCenter); … … 617 642 */ 618 643 CandidateForTesselation::~CandidateForTesselation() { 619 point = NULL;620 644 BaseLine = NULL; 621 645 }; … … 628 652 { 629 653 ost << "[" << a.BaseLine->Nr << "|" << a.BaseLine->endpoints[0]->node->Name << "," << a.BaseLine->endpoints[1]->node->Name << "] with "; 630 if (a.point == NULL)654 if (a.pointlist.empty()) 631 655 ost << "no candidate."; 632 else 633 ost << "candidate " << *(a.point) << " at angle " << (a.ShortestAngle)<< "."; 656 else { 657 ost << "candidate"; 658 if (a.pointlist.size() != 1) 659 ost << "s "; 660 else 661 ost << " "; 662 for (TesselPointList::const_iterator Runner = a.pointlist.begin(); Runner != a.pointlist.end(); Runner++) 663 ost << *(*Runner) << " "; 664 ost << " at angle " << (a.ShortestAngle)<< "."; 665 } 634 666 635 667 return ost; … … 649 681 InternalPointer(PointsOnBoundary.begin()) 650 682 { 683 Info FunctionInfo(__func__); 651 684 } 652 685 ; … … 657 690 Tesselation::~Tesselation() 658 691 { 659 Log() << Verbose(1) << "Free'ing TesselStruct ... " << endl; 692 Info FunctionInfo(__func__); 693 Log() << Verbose(0) << "Free'ing TesselStruct ... " << endl; 660 694 for (TriangleMap::iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) { 661 695 if (runner->second != NULL) { … … 665 699 eLog() << Verbose(1) << "The triangle " << runner->first << " has already been free'd." << endl; 666 700 } 667 Log() << Verbose( 1) << "This envelope was written to file " << TriangleFilesWritten << " times(s)." << endl;701 Log() << Verbose(0) << "This envelope was written to file " << TriangleFilesWritten << " times(s)." << endl; 668 702 } 669 703 ; … … 674 708 Vector * Tesselation::GetCenter(ofstream *out) const 675 709 { 710 Info FunctionInfo(__func__); 676 711 Vector *Center = new Vector(0.,0.,0.); 677 712 int num=0; … … 689 724 TesselPoint * Tesselation::GetPoint() const 690 725 { 726 Info FunctionInfo(__func__); 691 727 return (InternalPointer->second->node); 692 728 }; … … 697 733 TesselPoint * Tesselation::GetTerminalPoint() const 698 734 { 735 Info FunctionInfo(__func__); 699 736 PointMap::const_iterator Runner = PointsOnBoundary.end(); 700 737 Runner--; … … 707 744 void Tesselation::GoToNext() const 708 745 { 746 Info FunctionInfo(__func__); 709 747 if (InternalPointer != PointsOnBoundary.end()) 710 748 InternalPointer++; … … 716 754 void Tesselation::GoToPrevious() const 717 755 { 756 Info FunctionInfo(__func__); 718 757 if (InternalPointer != PointsOnBoundary.begin()) 719 758 InternalPointer--; … … 725 764 void Tesselation::GoToFirst() const 726 765 { 766 Info FunctionInfo(__func__); 727 767 InternalPointer = PointsOnBoundary.begin(); 728 768 }; … … 733 773 void Tesselation::GoToLast() const 734 774 { 775 Info FunctionInfo(__func__); 735 776 InternalPointer = PointsOnBoundary.end(); 736 777 InternalPointer--; … … 742 783 bool Tesselation::IsEmpty() const 743 784 { 785 Info FunctionInfo(__func__); 744 786 return (PointsOnBoundary.empty()); 745 787 }; … … 750 792 bool Tesselation::IsEnd() const 751 793 { 794 Info FunctionInfo(__func__); 752 795 return (InternalPointer == PointsOnBoundary.end()); 753 796 }; … … 762 805 Tesselation::GuessStartingTriangle() 763 806 { 807 Info FunctionInfo(__func__); 764 808 // 4b. create a starting triangle 765 809 // 4b1. create all distances … … 808 852 baseline->second.first->second->node->node, 809 853 baseline->second.second->second->node->node); 810 Log() << Verbose(2) << "Plane vector of candidate triangle is "; 811 PlaneVector.Output(); 812 Log() << Verbose(0) << endl; 854 Log() << Verbose(2) << "Plane vector of candidate triangle is " << PlaneVector << endl; 813 855 // 4. loop over all points 814 856 double sign = 0.; … … 826 868 if (fabs(distance) < 1e-4) // we need to have a small epsilon around 0 which is still ok 827 869 continue; 828 Log() << Verbose(3) << "Projection of " << checker->second->node->Name 829 << " yields distance of " << distance << "." << endl; 870 Log() << Verbose(2) << "Projection of " << checker->second->node->Name << " yields distance of " << distance << "." << endl; 830 871 tmp = distance / fabs(distance); 831 872 // 4b. Any have different sign to than before? (i.e. would lie outside convex hull with this starting triangle) … … 880 921 if (checker == PointsOnBoundary.end()) 881 922 { 882 Log() << Verbose( 0) << "Looks like we have a candidate!" << endl;923 Log() << Verbose(2) << "Looks like we have a candidate!" << endl; 883 924 break; 884 925 } … … 910 951 else 911 952 { 912 Log() << Verbose(1) << "No starting triangle found." << endl; 913 exit(255); 953 eLog() << Verbose(0) << "No starting triangle found." << endl; 914 954 } 915 955 } … … 931 971 void Tesselation::TesselateOnBoundary(const PointCloud * const cloud) 932 972 { 973 Info FunctionInfo(__func__); 933 974 bool flag; 934 975 PointMap::iterator winner; … … 949 990 // get peak point with respect to this base line's only triangle 950 991 BTS = baseline->second->triangles.begin()->second; // there is only one triangle so far 951 Log() << Verbose( 2) << "Current baseline is between " << *(baseline->second) << "." << endl;992 Log() << Verbose(0) << "Current baseline is between " << *(baseline->second) << "." << endl; 952 993 for (int i = 0; i < 3; i++) 953 994 if ((BTS->endpoints[i] != baseline->second->endpoints[0]) && (BTS->endpoints[i] != baseline->second->endpoints[1])) 954 995 peak = BTS->endpoints[i]; 955 Log() << Verbose( 3) << " and has peak " << *peak << "." << endl;996 Log() << Verbose(1) << " and has peak " << *peak << "." << endl; 956 997 957 998 // prepare some auxiliary vectors … … 968 1009 CenterVector.AddVector(BTS->endpoints[i]->node->node); 969 1010 CenterVector.Scale(1. / 3.); 970 Log() << Verbose( 4) << "CenterVector of base triangle is " << CenterVector << endl;1011 Log() << Verbose(2) << "CenterVector of base triangle is " << CenterVector << endl; 971 1012 972 1013 // normal vector of triangle … … 975 1016 BTS->GetNormalVector(NormalVector); 976 1017 NormalVector.CopyVector(&BTS->NormalVector); 977 Log() << Verbose( 4) << "NormalVector of base triangle is " << NormalVector << endl;1018 Log() << Verbose(2) << "NormalVector of base triangle is " << NormalVector << endl; 978 1019 979 1020 // vector in propagation direction (out of triangle) … … 982 1023 TempVector.CopyVector(&CenterVector); 983 1024 TempVector.SubtractVector(baseline->second->endpoints[0]->node->node); // TempVector is vector on triangle plane pointing from one baseline egde towards center! 984 //Log() << Verbose( 2) << "Projection of propagation onto temp: " << PropagationVector.Projection(&TempVector) << "." << endl;1025 //Log() << Verbose(0) << "Projection of propagation onto temp: " << PropagationVector.Projection(&TempVector) << "." << endl; 985 1026 if (PropagationVector.ScalarProduct(&TempVector) > 0) // make sure normal propagation vector points outward from baseline 986 1027 PropagationVector.Scale(-1.); 987 Log() << Verbose( 4) << "PropagationVector of base triangle is " << PropagationVector << endl;1028 Log() << Verbose(2) << "PropagationVector of base triangle is " << PropagationVector << endl; 988 1029 winner = PointsOnBoundary.end(); 989 1030 … … 991 1032 for (PointMap::iterator target = PointsOnBoundary.begin(); target != PointsOnBoundary.end(); target++) { 992 1033 if ((target->second != baseline->second->endpoints[0]) && (target->second != baseline->second->endpoints[1])) { // don't take the same endpoints 993 Log() << Verbose( 3) << "Target point is " << *(target->second) << ":" << endl;1034 Log() << Verbose(1) << "Target point is " << *(target->second) << ":" << endl; 994 1035 995 1036 // first check direction, so that triangles don't intersect … … 998 1039 VirtualNormalVector.ProjectOntoPlane(&NormalVector); 999 1040 TempAngle = VirtualNormalVector.Angle(&PropagationVector); 1000 Log() << Verbose( 4) << "VirtualNormalVector is " << VirtualNormalVector << " and PropagationVector is " << PropagationVector << "." << endl;1041 Log() << Verbose(2) << "VirtualNormalVector is " << VirtualNormalVector << " and PropagationVector is " << PropagationVector << "." << endl; 1001 1042 if (TempAngle > (M_PI/2.)) { // no bends bigger than Pi/2 (90 degrees) 1002 Log() << Verbose( 4) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", bad direction!" << endl;1043 Log() << Verbose(2) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", bad direction!" << endl; 1003 1044 continue; 1004 1045 } else 1005 Log() << Verbose( 4) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", good direction!" << endl;1046 Log() << Verbose(2) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", good direction!" << endl; 1006 1047 1007 1048 // check first and second endpoint (if any connecting line goes to target has at least not more than 1 triangle) … … 1009 1050 LineChecker[1] = baseline->second->endpoints[1]->lines.find(target->first); 1010 1051 if (((LineChecker[0] != baseline->second->endpoints[0]->lines.end()) && (LineChecker[0]->second->triangles.size() == 2))) { 1011 Log() << Verbose( 4) << *(baseline->second->endpoints[0]) << " has line " << *(LineChecker[0]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[0]->second->triangles.size() << " triangles." << endl;1052 Log() << Verbose(2) << *(baseline->second->endpoints[0]) << " has line " << *(LineChecker[0]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[0]->second->triangles.size() << " triangles." << endl; 1012 1053 continue; 1013 1054 } 1014 1055 if (((LineChecker[1] != baseline->second->endpoints[1]->lines.end()) && (LineChecker[1]->second->triangles.size() == 2))) { 1015 Log() << Verbose( 4) << *(baseline->second->endpoints[1]) << " has line " << *(LineChecker[1]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[1]->second->triangles.size() << " triangles." << endl;1056 Log() << Verbose(2) << *(baseline->second->endpoints[1]) << " has line " << *(LineChecker[1]->second) << " to " << *(target->second) << " as endpoint with " << LineChecker[1]->second->triangles.size() << " triangles." << endl; 1016 1057 continue; 1017 1058 } … … 1030 1071 helper.ProjectOntoPlane(&TempVector); 1031 1072 if (fabs(helper.NormSquared()) < MYEPSILON) { 1032 Log() << Verbose( 4) << "Chosen set of vectors is linear dependent." << endl;1073 Log() << Verbose(2) << "Chosen set of vectors is linear dependent." << endl; 1033 1074 continue; 1034 1075 } … … 1047 1088 // calculate angle 1048 1089 TempAngle = NormalVector.Angle(&VirtualNormalVector); 1049 Log() << Verbose( 4) << "NormalVector is " << VirtualNormalVector << " and the angle is " << TempAngle << "." << endl;1090 Log() << Verbose(2) << "NormalVector is " << VirtualNormalVector << " and the angle is " << TempAngle << "." << endl; 1050 1091 if ((SmallestAngle - TempAngle) > MYEPSILON) { // set to new possible winner 1051 1092 SmallestAngle = TempAngle; 1052 1093 winner = target; 1053 Log() << Verbose( 4) << "New winner " << *winner->second->node << " due to smaller angle between normal vectors." << endl;1094 Log() << Verbose(2) << "New winner " << *winner->second->node << " due to smaller angle between normal vectors." << endl; 1054 1095 } else if (fabs(SmallestAngle - TempAngle) < MYEPSILON) { // check the angle to propagation, both possible targets are in one plane! (their normals have same angle) 1055 1096 // hence, check the angles to some normal direction from our base line but in this common plane of both targets... … … 1069 1110 SmallestAngle = TempAngle; 1070 1111 winner = target; 1071 Log() << Verbose( 4) << "New winner " << *winner->second->node << " due to smaller angle " << TempAngle << " to propagation direction." << endl;1112 Log() << Verbose(2) << "New winner " << *winner->second->node << " due to smaller angle " << TempAngle << " to propagation direction." << endl; 1072 1113 } else 1073 Log() << Verbose( 4) << "Keeping old winner " << *winner->second->node << " due to smaller angle to propagation direction." << endl;1114 Log() << Verbose(2) << "Keeping old winner " << *winner->second->node << " due to smaller angle to propagation direction." << endl; 1074 1115 } else 1075 Log() << Verbose( 4) << "Keeping old winner " << *winner->second->node << " due to smaller angle between normal vectors." << endl;1116 Log() << Verbose(2) << "Keeping old winner " << *winner->second->node << " due to smaller angle between normal vectors." << endl; 1076 1117 } 1077 1118 } // end of loop over all boundary points … … 1079 1120 // 5b. The point of the above whose triangle has the greatest angle with the triangle the current line belongs to (it only belongs to one, remember!): New triangle 1080 1121 if (winner != PointsOnBoundary.end()) { 1081 Log() << Verbose( 2) << "Winning target point is " << *(winner->second) << " with angle " << SmallestAngle << "." << endl;1122 Log() << Verbose(0) << "Winning target point is " << *(winner->second) << " with angle " << SmallestAngle << "." << endl; 1082 1123 // create the lins of not yet present 1083 1124 BLS[0] = baseline->second; … … 1109 1150 TrianglesOnBoundaryCount++; 1110 1151 } else { 1111 Log() << Verbose(1) << "I could not determine a winner for this baseline " << *(baseline->second) << "." << endl;1152 eLog() << Verbose(2) << "I could not determine a winner for this baseline " << *(baseline->second) << "." << endl; 1112 1153 } 1113 1154 1114 1155 // 5d. If the set of lines is not yet empty, go to 5. and continue 1115 1156 } else 1116 Log() << Verbose( 2) << "Baseline candidate " << *(baseline->second) << " has a triangle count of " << baseline->second->triangles.size() << "." << endl;1157 Log() << Verbose(0) << "Baseline candidate " << *(baseline->second) << " has a triangle count of " << baseline->second->triangles.size() << "." << endl; 1117 1158 } while (flag); 1118 1159 … … 1129 1170 bool Tesselation::InsertStraddlingPoints(const PointCloud *cloud, const LinkedCell *LC) 1130 1171 { 1172 Info FunctionInfo(__func__); 1131 1173 Vector Intersection, Normal; 1132 1174 TesselPoint *Walker = NULL; … … 1135 1177 bool AddFlag = false; 1136 1178 LinkedCell *BoundaryPoints = NULL; 1137 1138 Log() << Verbose(1) << "Begin of InsertStraddlingPoints" << endl;1139 1179 1140 1180 cloud->GoToFirst(); … … 1147 1187 } 1148 1188 Walker = cloud->GetPoint(); 1149 Log() << Verbose( 2) << "Current point is " << *Walker << "." << endl;1189 Log() << Verbose(0) << "Current point is " << *Walker << "." << endl; 1150 1190 // get the next triangle 1151 1191 triangles = FindClosestTrianglesToPoint(Walker->node, BoundaryPoints); 1152 1192 BTS = triangles->front(); 1153 1193 if ((triangles == NULL) || (BTS->ContainsBoundaryPoint(Walker))) { 1154 Log() << Verbose( 2) << "No triangles found, probably a tesselation point itself." << endl;1194 Log() << Verbose(0) << "No triangles found, probably a tesselation point itself." << endl; 1155 1195 cloud->GoToNext(); 1156 1196 continue; 1157 1197 } else { 1158 1198 } 1159 Log() << Verbose( 2) << "Closest triangle is " << *BTS << "." << endl;1199 Log() << Verbose(0) << "Closest triangle is " << *BTS << "." << endl; 1160 1200 // get the intersection point 1161 1201 if (BTS->GetIntersectionInsideTriangle(Center, Walker->node, &Intersection)) { 1162 Log() << Verbose( 2) << "We have an intersection at " << Intersection << "." << endl;1202 Log() << Verbose(0) << "We have an intersection at " << Intersection << "." << endl; 1163 1203 // we have the intersection, check whether in- or outside of boundary 1164 1204 if ((Center->DistanceSquared(Walker->node) - Center->DistanceSquared(&Intersection)) < -MYEPSILON) { 1165 1205 // inside, next! 1166 Log() << Verbose( 2) << *Walker << " is inside wrt triangle " << *BTS << "." << endl;1206 Log() << Verbose(0) << *Walker << " is inside wrt triangle " << *BTS << "." << endl; 1167 1207 } else { 1168 1208 // outside! 1169 Log() << Verbose( 2) << *Walker << " is outside wrt triangle " << *BTS << "." << endl;1209 Log() << Verbose(0) << *Walker << " is outside wrt triangle " << *BTS << "." << endl; 1170 1210 class BoundaryLineSet *OldLines[3], *NewLines[3]; 1171 1211 class BoundaryPointSet *OldPoints[3], *NewPoint; … … 1177 1217 Normal.CopyVector(&BTS->NormalVector); 1178 1218 // add Walker to boundary points 1179 Log() << Verbose( 2) << "Adding " << *Walker << " to BoundaryPoints." << endl;1219 Log() << Verbose(0) << "Adding " << *Walker << " to BoundaryPoints." << endl; 1180 1220 AddFlag = true; 1181 1221 if (AddBoundaryPoint(Walker,0)) … … 1184 1224 continue; 1185 1225 // remove triangle 1186 Log() << Verbose( 2) << "Erasing triangle " << *BTS << "." << endl;1226 Log() << Verbose(0) << "Erasing triangle " << *BTS << "." << endl; 1187 1227 TrianglesOnBoundary.erase(BTS->Nr); 1188 1228 delete(BTS); … … 1192 1232 BPS[1] = OldPoints[i]; 1193 1233 NewLines[i] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount); 1194 Log() << Verbose( 3) << "Creating new line " << *NewLines[i] << "." << endl;1234 Log() << Verbose(1) << "Creating new line " << *NewLines[i] << "." << endl; 1195 1235 LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, NewLines[i])); // no need for check for unique insertion as BPS[0] is definitely a new one 1196 1236 LinesOnBoundaryCount++; … … 1203 1243 if (NewLines[j]->IsConnectedTo(BLS[0])) { 1204 1244 if (n>2) { 1205 Log() << Verbose(1) << BLS[0] << " connects to all of the new lines?!" << endl;1245 eLog() << Verbose(2) << BLS[0] << " connects to all of the new lines?!" << endl; 1206 1246 return false; 1207 1247 } else … … 1214 1254 BTS->GetNormalVector(Normal); 1215 1255 Normal.Scale(-1.); 1216 Log() << Verbose( 2) << "Created new triangle " << *BTS << "." << endl;1256 Log() << Verbose(0) << "Created new triangle " << *BTS << "." << endl; 1217 1257 TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS)); 1218 1258 TrianglesOnBoundaryCount++; … … 1228 1268 // exit 1229 1269 delete(Center); 1230 Log() << Verbose(1) << "End of InsertStraddlingPoints" << endl;1231 1270 return true; 1232 1271 }; … … 1239 1278 bool Tesselation::AddBoundaryPoint(TesselPoint * Walker, const int n) 1240 1279 { 1280 Info FunctionInfo(__func__); 1241 1281 PointTestPair InsertUnique; 1242 1282 BPS[n] = new class BoundaryPointSet(Walker); … … 1260 1300 void Tesselation::AddTesselationPoint(TesselPoint* Candidate, const int n) 1261 1301 { 1302 Info FunctionInfo(__func__); 1262 1303 PointTestPair InsertUnique; 1263 1304 TPS[n] = new class BoundaryPointSet(Candidate); … … 1267 1308 } else { 1268 1309 delete TPS[n]; 1269 Log() << Verbose( 4) << "Node " << *((InsertUnique.first)->second->node) << " is already present in PointsOnBoundary." << endl;1310 Log() << Verbose(0) << "Node " << *((InsertUnique.first)->second->node) << " is already present in PointsOnBoundary." << endl; 1270 1311 TPS[n] = (InsertUnique.first)->second; 1271 1312 } … … 1280 1321 void Tesselation::SetTesselationPoint(TesselPoint* Candidate, const int n) const 1281 1322 { 1323 Info FunctionInfo(__func__); 1282 1324 PointMap::const_iterator FindPoint = PointsOnBoundary.find(Candidate->nr); 1283 1325 if (FindPoint != PointsOnBoundary.end()) … … 1301 1343 pair<LineMap::iterator,LineMap::iterator> FindPair; 1302 1344 FindPair = a->lines.equal_range(b->node->nr); 1303 Log() << Verbose( 5) << "INFO: There is at least one line between " << *a << " and " << *b << ": " << *(FindLine->second) << "." << endl;1345 Log() << Verbose(1) << "INFO: There is at least one line between " << *a << " and " << *b << ": " << *(FindLine->second) << "." << endl; 1304 1346 1305 1347 for (FindLine = FindPair.first; FindLine != FindPair.second; FindLine++) { … … 1307 1349 if (FindLine->second->triangles.size() < 2) { 1308 1350 insertNewLine = false; 1309 Log() << Verbose( 4) << "Using existing line " << *FindLine->second << endl;1351 Log() << Verbose(0) << "Using existing line " << *FindLine->second << endl; 1310 1352 1311 1353 BPS[0] = FindLine->second->endpoints[0]; … … 1339 1381 void Tesselation::AlwaysAddTesselationTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, const int n) 1340 1382 { 1341 Log() << Verbose(4) << "Adding open line [" << LinesOnBoundaryCount << "|" << *(a->node) << " and " << *(b->node) << "." << endl; 1383 Info FunctionInfo(__func__); 1384 Log() << Verbose(0) << "Adding open line [" << LinesOnBoundaryCount << "|" << *(a->node) << " and " << *(b->node) << "." << endl; 1342 1385 BPS[0] = a; 1343 1386 BPS[1] = b; … … 1357 1400 void Tesselation::AddTesselationTriangle() 1358 1401 { 1402 Info FunctionInfo(__func__); 1359 1403 Log() << Verbose(1) << "Adding triangle to global TrianglesOnBoundary map." << endl; 1360 1404 … … 1375 1419 void Tesselation::AddTesselationTriangle(const int nr) 1376 1420 { 1377 Log() << Verbose(1) << "Adding triangle to global TrianglesOnBoundary map." << endl; 1421 Info FunctionInfo(__func__); 1422 Log() << Verbose(0) << "Adding triangle to global TrianglesOnBoundary map." << endl; 1378 1423 1379 1424 // add triangle to global map … … 1393 1438 void Tesselation::RemoveTesselationTriangle(class BoundaryTriangleSet *triangle) 1394 1439 { 1440 Info FunctionInfo(__func__); 1395 1441 if (triangle == NULL) 1396 1442 return; 1397 1443 for (int i = 0; i < 3; i++) { 1398 1444 if (triangle->lines[i] != NULL) { 1399 Log() << Verbose( 5) << "Removing triangle Nr." << triangle->Nr << " in line " << *triangle->lines[i] << "." << endl;1445 Log() << Verbose(0) << "Removing triangle Nr." << triangle->Nr << " in line " << *triangle->lines[i] << "." << endl; 1400 1446 triangle->lines[i]->triangles.erase(triangle->Nr); 1401 1447 if (triangle->lines[i]->triangles.empty()) { 1402 Log() << Verbose( 5) << *triangle->lines[i] << " is no more attached to any triangle, erasing." << endl;1448 Log() << Verbose(0) << *triangle->lines[i] << " is no more attached to any triangle, erasing." << endl; 1403 1449 RemoveTesselationLine(triangle->lines[i]); 1404 1450 } else { 1405 Log() << Verbose( 5) << *triangle->lines[i] << " is still attached to another triangle: ";1451 Log() << Verbose(0) << *triangle->lines[i] << " is still attached to another triangle: "; 1406 1452 for(TriangleMap::iterator TriangleRunner = triangle->lines[i]->triangles.begin(); TriangleRunner != triangle->lines[i]->triangles.end(); TriangleRunner++) 1407 1453 Log() << Verbose(0) << "[" << (TriangleRunner->second)->Nr << "|" << *((TriangleRunner->second)->endpoints[0]) << ", " << *((TriangleRunner->second)->endpoints[1]) << ", " << *((TriangleRunner->second)->endpoints[2]) << "] \t"; 1408 1454 Log() << Verbose(0) << endl; 1409 1455 // for (int j=0;j<2;j++) { 1410 // Log() << Verbose( 5) << "Lines of endpoint " << *(triangle->lines[i]->endpoints[j]) << ": ";1456 // Log() << Verbose(0) << "Lines of endpoint " << *(triangle->lines[i]->endpoints[j]) << ": "; 1411 1457 // for(LineMap::iterator LineRunner = triangle->lines[i]->endpoints[j]->lines.begin(); LineRunner != triangle->lines[i]->endpoints[j]->lines.end(); LineRunner++) 1412 1458 // Log() << Verbose(0) << "[" << *(LineRunner->second) << "] \t"; … … 1420 1466 1421 1467 if (TrianglesOnBoundary.erase(triangle->Nr)) 1422 Log() << Verbose( 5) << "Removing triangle Nr. " << triangle->Nr << "." << endl;1468 Log() << Verbose(0) << "Removing triangle Nr. " << triangle->Nr << "." << endl; 1423 1469 delete(triangle); 1424 1470 }; … … 1430 1476 void Tesselation::RemoveTesselationLine(class BoundaryLineSet *line) 1431 1477 { 1478 Info FunctionInfo(__func__); 1432 1479 int Numbers[2]; 1433 1480 … … 1450 1497 for (LineMap::iterator Runner = erasor.first; Runner != erasor.second; Runner++) 1451 1498 if ((*Runner).second == line) { 1452 Log() << Verbose( 5) << "Removing Line Nr. " << line->Nr << " in boundary point " << *line->endpoints[i] << "." << endl;1499 Log() << Verbose(0) << "Removing Line Nr. " << line->Nr << " in boundary point " << *line->endpoints[i] << "." << endl; 1453 1500 line->endpoints[i]->lines.erase(Runner); 1454 1501 break; … … 1456 1503 } else { // there's just a single line left 1457 1504 if (line->endpoints[i]->lines.erase(line->Nr)) 1458 Log() << Verbose( 5) << "Removing Line Nr. " << line->Nr << " in boundary point " << *line->endpoints[i] << "." << endl;1505 Log() << Verbose(0) << "Removing Line Nr. " << line->Nr << " in boundary point " << *line->endpoints[i] << "." << endl; 1459 1506 } 1460 1507 if (line->endpoints[i]->lines.empty()) { 1461 Log() << Verbose( 5) << *line->endpoints[i] << " has no more lines it's attached to, erasing." << endl;1508 Log() << Verbose(0) << *line->endpoints[i] << " has no more lines it's attached to, erasing." << endl; 1462 1509 RemoveTesselationPoint(line->endpoints[i]); 1463 1510 } else { 1464 Log() << Verbose( 5) << *line->endpoints[i] << " has still lines it's attached to: ";1511 Log() << Verbose(0) << *line->endpoints[i] << " has still lines it's attached to: "; 1465 1512 for(LineMap::iterator LineRunner = line->endpoints[i]->lines.begin(); LineRunner != line->endpoints[i]->lines.end(); LineRunner++) 1466 1513 Log() << Verbose(0) << "[" << *(LineRunner->second) << "] \t"; … … 1475 1522 1476 1523 if (LinesOnBoundary.erase(line->Nr)) 1477 Log() << Verbose( 5) << "Removing line Nr. " << line->Nr << "." << endl;1524 Log() << Verbose(0) << "Removing line Nr. " << line->Nr << "." << endl; 1478 1525 delete(line); 1479 1526 }; … … 1486 1533 void Tesselation::RemoveTesselationPoint(class BoundaryPointSet *point) 1487 1534 { 1535 Info FunctionInfo(__func__); 1488 1536 if (point == NULL) 1489 1537 return; 1490 1538 if (PointsOnBoundary.erase(point->Nr)) 1491 Log() << Verbose( 5) << "Removing point Nr. " << point->Nr << "." << endl;1539 Log() << Verbose(0) << "Removing point Nr. " << point->Nr << "." << endl; 1492 1540 delete(point); 1493 1541 }; … … 1504 1552 int Tesselation::CheckPresenceOfTriangle(TesselPoint *Candidates[3]) const 1505 1553 { 1554 Info FunctionInfo(__func__); 1506 1555 int adjacentTriangleCount = 0; 1507 1556 class BoundaryPointSet *Points[3]; 1508 1557 1509 Log() << Verbose(2) << "Begin of CheckPresenceOfTriangle" << endl;1510 1558 // builds a triangle point set (Points) of the end points 1511 1559 for (int i = 0; i < 3; i++) { … … 1526 1574 for (; (FindLine != Points[i]->lines.end()) && (FindLine->first == Points[j]->node->nr); FindLine++) { 1527 1575 TriangleMap *triangles = &FindLine->second->triangles; 1528 Log() << Verbose( 3) << "Current line is " << FindLine->first << ": " << *(FindLine->second) << " with triangles " << triangles << "." << endl;1576 Log() << Verbose(1) << "Current line is " << FindLine->first << ": " << *(FindLine->second) << " with triangles " << triangles << "." << endl; 1529 1577 for (TriangleMap::const_iterator FindTriangle = triangles->begin(); FindTriangle != triangles->end(); FindTriangle++) { 1530 1578 if (FindTriangle->second->IsPresentTupel(Points)) { … … 1532 1580 } 1533 1581 } 1534 Log() << Verbose( 3) << "end." << endl;1582 Log() << Verbose(1) << "end." << endl; 1535 1583 } 1536 1584 // Only one of the triangle lines must be considered for the triangle count. 1537 //Log() << Verbose( 2) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;1585 //Log() << Verbose(0) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl; 1538 1586 //return adjacentTriangleCount; 1539 1587 } … … 1542 1590 } 1543 1591 1544 Log() << Verbose(2) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl; 1545 Log() << Verbose(2) << "End of CheckPresenceOfTriangle" << endl; 1592 Log() << Verbose(0) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl; 1546 1593 return adjacentTriangleCount; 1547 1594 }; … … 1557 1604 class BoundaryTriangleSet * Tesselation::GetPresentTriangle(TesselPoint *Candidates[3]) 1558 1605 { 1606 Info FunctionInfo(__func__); 1559 1607 class BoundaryTriangleSet *triangle = NULL; 1560 1608 class BoundaryPointSet *Points[3]; … … 1586 1634 } 1587 1635 // Only one of the triangle lines must be considered for the triangle count. 1588 //Log() << Verbose( 2) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;1636 //Log() << Verbose(0) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl; 1589 1637 //return adjacentTriangleCount; 1590 1638 } … … 1607 1655 void Tesselation::FindStartingTriangle(const double RADIUS, const LinkedCell *LC) 1608 1656 { 1609 Log() << Verbose(1) << "Begin of FindStartingTriangle\n";1657 Info FunctionInfo(__func__); 1610 1658 int i = 0; 1611 1659 TesselPoint* FirstPoint = NULL; … … 1631 1679 for (LC->n[(i+2)%NDIM]=0;LC->n[(i+2)%NDIM]<LC->N[(i+2)%NDIM];LC->n[(i+2)%NDIM]++) { 1632 1680 const LinkedNodes *List = LC->GetCurrentCell(); 1633 //Log() << Verbose( 2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;1681 //Log() << Verbose(1) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl; 1634 1682 if (List != NULL) { 1635 1683 for (LinkedNodes::const_iterator Runner = List->begin();Runner != List->end();Runner++) { 1636 1684 if ((*Runner)->node->x[i] > maxCoordinate[i]) { 1637 Log() << Verbose( 2) << "New maximal for axis " << i << " node is " << *(*Runner) << " at " << *(*Runner)->node << "." << endl;1685 Log() << Verbose(1) << "New maximal for axis " << i << " node is " << *(*Runner) << " at " << *(*Runner)->node << "." << endl; 1638 1686 maxCoordinate[i] = (*Runner)->node->x[i]; 1639 1687 MaxPoint[i] = (*Runner); … … 1646 1694 } 1647 1695 1648 Log() << Verbose( 2) << "Found maximum coordinates: ";1696 Log() << Verbose(1) << "Found maximum coordinates: "; 1649 1697 for (int i=0;i<NDIM;i++) 1650 1698 Log() << Verbose(0) << i << ": " << *MaxPoint[i] << "\t"; … … 1652 1700 1653 1701 BTS = NULL; 1654 CandidateList *OptCandidates = new CandidateList();1655 1702 for (int k=0;k<NDIM;k++) { 1656 1703 Oben.Zero(); 1657 1704 Oben.x[k] = 1.; 1658 1705 FirstPoint = MaxPoint[k]; 1659 Log() << Verbose( 1) << "Coordinates of start node at " << *FirstPoint->node << "." << endl;1706 Log() << Verbose(0) << "Coordinates of start node at " << *FirstPoint->node << "." << endl; 1660 1707 1661 1708 double ShortestAngle; … … 1688 1735 1689 1736 // adding point 1 and point 2 and add the line between them 1690 Log() << Verbose( 1) << "Coordinates of start node at " << *FirstPoint->node << "." << endl;1737 Log() << Verbose(0) << "Coordinates of start node at " << *FirstPoint->node << "." << endl; 1691 1738 AddTesselationPoint(FirstPoint, 0); 1692 Log() << Verbose( 1) << "Found second point is at " << *SecondPoint->node << ".\n";1739 Log() << Verbose(0) << "Found second point is at " << *SecondPoint->node << ".\n"; 1693 1740 AddTesselationPoint(SecondPoint, 1); 1694 1741 AddTesselationLine(TPS[0], TPS[1], 0); 1695 1742 1696 //Log() << Verbose( 2) << "INFO: OldSphereCenter is at " << helper << ".\n";1697 FindThirdPointForTesselation(Oben, SearchDirection, helper, BLS[0], NULL, *&OptCandidates, &ShortestAngle, RADIUS, LC);1698 Log() << Verbose(1) << "List of third Points is ";1699 for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) {1700 Log() << Verbose(0) << " " << *(*it)->point;1701 }1702 Log() << Verbose(0) << endl;1703 1704 for ( CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) {1743 //Log() << Verbose(1) << "INFO: OldSphereCenter is at " << helper << ".\n"; 1744 CandidateForTesselation OptCandidates(BLS[0]); 1745 FindThirdPointForTesselation(Oben, SearchDirection, helper, OptCandidates, NULL, RADIUS, LC); 1746 Log() << Verbose(0) << "List of third Points is:" << endl; 1747 for (TesselPointList::iterator it = OptCandidates.pointlist.begin(); it != OptCandidates.pointlist.end(); it++) { 1748 Log() << Verbose(0) << " " << *(*it) << endl; 1749 } 1750 1751 for (TesselPointList::iterator it = OptCandidates.pointlist.begin(); it != OptCandidates.pointlist.end(); it++) { 1705 1752 // add third triangle point 1706 AddTesselationPoint((*it) ->point, 2);1753 AddTesselationPoint((*it), 2); 1707 1754 // add the second and third line 1708 1755 AddTesselationLine(TPS[1], TPS[2], 1); … … 1712 1759 AddTesselationTriangle(); 1713 1760 // ... and calculate its normal vector (with correct orientation) 1714 (*it)->OptCenter.Scale(-1.); 1715 Log() << Verbose(2) << "Anti-Oben is currently " << (*it)->OptCenter << "." << endl; 1716 BTS->GetNormalVector((*it)->OptCenter); // vector to compare with should point inwards 1761 OptCandidates.OptCenter.Scale(-1.); 1762 Log() << Verbose(1) << "Anti-Oben is currently " << OptCandidates.OptCenter << "." << endl; 1763 BTS->GetNormalVector(OptCandidates.OptCenter); // vector to compare with should point inwards 1764 OptCandidates.OptCenter.Scale(-1.); 1717 1765 Log() << Verbose(0) << "==> Found starting triangle consists of " << *FirstPoint << ", " << *SecondPoint << " and " 1718 << *(*it)->point<< " with normal vector " << BTS->NormalVector << ".\n";1766 << (*it) << " with normal vector " << BTS->NormalVector << ".\n"; 1719 1767 1720 1768 // // if we do not reach the end with the next step of iteration, we need to setup a new first line … … 1727 1775 // AddTesselationLine(TPS[0], TPS[1], 0); 1728 1776 // } 1729 Log() << Verbose( 2) << "Projection is " << BTS->NormalVector.ScalarProduct(&Oben) << "." << endl;1777 Log() << Verbose(1) << "Projection is " << BTS->NormalVector.ScalarProduct(&Oben) << "." << endl; 1730 1778 } 1731 1779 if (BTS != NULL) // we have created one starting triangle … … 1733 1781 else { 1734 1782 // remove all candidates from the list and then the list itself 1735 class CandidateForTesselation *remover = NULL; 1736 for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) { 1737 remover = *it; 1738 delete(remover); 1739 } 1740 OptCandidates->clear(); 1741 } 1742 } 1743 1744 // remove all candidates from the list and then the list itself 1745 class CandidateForTesselation *remover = NULL; 1746 for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) { 1747 remover = *it; 1748 delete(remover); 1749 } 1750 delete(OptCandidates); 1751 Log() << Verbose(1) << "End of FindStartingTriangle\n"; 1783 } 1784 } 1752 1785 }; 1753 1786 1754 1787 /** Checks for a given baseline and a third point candidate whether baselines of the found triangle don't have even better candidates. 1755 1788 * This is supposed to prevent early closing of the tesselation. 1756 * \param *BaseRay baseline, i.e. not \a *OptCandidate1789 * \param CandidateLine CandidateForTesselation with baseline and shortestangle , i.e. not \a *OptCandidate 1757 1790 * \param *ThirdNode third point in triangle, not in BoundaryLineSet::endpoints 1758 * \param ShortestAngle path length on this circle band for the current \a *ThirdNode1759 1791 * \param RADIUS radius of sphere 1760 1792 * \param *LC LinkedCell structure 1761 1793 * \return true - there is a better candidate (smaller angle than \a ShortestAngle), false - no better TesselPoint candidate found 1762 1794 */ 1763 bool Tesselation::HasOtherBaselineBetterCandidate(const BoundaryLineSet * const BaseRay, const TesselPoint * const ThirdNode, double ShortestAngle, double RADIUS, const LinkedCell * const LC) const 1764 { 1765 bool result = false; 1766 Vector CircleCenter; 1767 Vector CirclePlaneNormal; 1768 Vector OldSphereCenter; 1769 Vector SearchDirection; 1770 Vector helper; 1771 TesselPoint *OtherOptCandidate = NULL; 1772 double OtherShortestAngle = 2.*M_PI; // This will indicate the quadrant. 1773 double radius, CircleRadius; 1774 BoundaryLineSet *Line = NULL; 1775 BoundaryTriangleSet *T = NULL; 1776 1777 Log() << Verbose(1) << "Begin of HasOtherBaselineBetterCandidate" << endl; 1778 1779 // check both other lines 1780 PointMap::const_iterator FindPoint = PointsOnBoundary.find(ThirdNode->nr); 1781 if (FindPoint != PointsOnBoundary.end()) { 1782 for (int i=0;i<2;i++) { 1783 LineMap::const_iterator FindLine = (FindPoint->second)->lines.find(BaseRay->endpoints[0]->node->nr); 1784 if (FindLine != (FindPoint->second)->lines.end()) { 1785 Line = FindLine->second; 1786 Log() << Verbose(1) << "Found line " << *Line << "." << endl; 1787 if (Line->triangles.size() == 1) { 1788 T = Line->triangles.begin()->second; 1789 // construct center of circle 1790 CircleCenter.CopyVector(Line->endpoints[0]->node->node); 1791 CircleCenter.AddVector(Line->endpoints[1]->node->node); 1792 CircleCenter.Scale(0.5); 1793 1794 // construct normal vector of circle 1795 CirclePlaneNormal.CopyVector(Line->endpoints[0]->node->node); 1796 CirclePlaneNormal.SubtractVector(Line->endpoints[1]->node->node); 1797 1798 // calculate squared radius of circle 1799 radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal); 1800 if (radius/4. < RADIUS*RADIUS) { 1801 CircleRadius = RADIUS*RADIUS - radius/4.; 1802 CirclePlaneNormal.Normalize(); 1803 //Log() << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl; 1804 1805 // construct old center 1806 GetCenterofCircumcircle(&OldSphereCenter, *T->endpoints[0]->node->node, *T->endpoints[1]->node->node, *T->endpoints[2]->node->node); 1807 helper.CopyVector(&T->NormalVector); // normal vector ensures that this is correct center of the two possible ones 1808 radius = Line->endpoints[0]->node->node->DistanceSquared(&OldSphereCenter); 1809 helper.Scale(sqrt(RADIUS*RADIUS - radius)); 1810 OldSphereCenter.AddVector(&helper); 1811 OldSphereCenter.SubtractVector(&CircleCenter); 1812 //Log() << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl; 1813 1814 // construct SearchDirection 1815 SearchDirection.MakeNormalVector(&T->NormalVector, &CirclePlaneNormal); 1816 helper.CopyVector(Line->endpoints[0]->node->node); 1817 helper.SubtractVector(ThirdNode->node); 1818 if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards! 1819 SearchDirection.Scale(-1.); 1820 SearchDirection.ProjectOntoPlane(&OldSphereCenter); 1821 SearchDirection.Normalize(); 1822 Log() << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl; 1823 if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) { 1824 // rotated the wrong way! 1825 eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl; 1826 } 1827 1828 // add third point 1829 CandidateList *OptCandidates = new CandidateList(); 1830 FindThirdPointForTesselation(T->NormalVector, SearchDirection, OldSphereCenter, Line, ThirdNode, OptCandidates, &OtherShortestAngle, RADIUS, LC); 1831 for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) { 1832 if (((*it)->point == BaseRay->endpoints[0]->node) || ((*it)->point == BaseRay->endpoints[1]->node)) // skip if it's the same triangle than suggested 1833 continue; 1834 Log() << Verbose(1) << " Third point candidate is " << *(*it)->point 1835 << " with circumsphere's center at " << (*it)->OptCenter << "." << endl; 1836 Log() << Verbose(1) << " Baseline is " << *BaseRay << endl; 1837 1838 // check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2) 1839 TesselPoint *PointCandidates[3]; 1840 PointCandidates[0] = (*it)->point; 1841 PointCandidates[1] = BaseRay->endpoints[0]->node; 1842 PointCandidates[2] = BaseRay->endpoints[1]->node; 1843 bool check=false; 1844 int existentTrianglesCount = CheckPresenceOfTriangle(PointCandidates); 1845 // If there is no triangle, add it regularly. 1846 if (existentTrianglesCount == 0) { 1847 SetTesselationPoint((*it)->point, 0); 1848 SetTesselationPoint(BaseRay->endpoints[0]->node, 1); 1849 SetTesselationPoint(BaseRay->endpoints[1]->node, 2); 1850 1851 if (CheckLineCriteriaForDegeneratedTriangle((const BoundaryPointSet ** const )TPS)) { 1852 OtherOptCandidate = (*it)->point; 1853 check = true; 1854 } 1855 } else if ((existentTrianglesCount >= 1) && (existentTrianglesCount <= 3)) { // If there is a planar region within the structure, we need this triangle a second time. 1856 SetTesselationPoint((*it)->point, 0); 1857 SetTesselationPoint(BaseRay->endpoints[0]->node, 1); 1858 SetTesselationPoint(BaseRay->endpoints[1]->node, 2); 1859 1860 // We demand that at most one new degenerate line is created and that this line also already exists (which has to be the case due to existentTrianglesCount == 1) 1861 // i.e. at least one of the three lines must be present with TriangleCount <= 1 1862 if (CheckLineCriteriaForDegeneratedTriangle((const BoundaryPointSet ** const)TPS)) { 1863 OtherOptCandidate = (*it)->point; 1864 check = true; 1865 } 1866 } 1867 1868 if (check) { 1869 if (ShortestAngle > OtherShortestAngle) { 1870 Log() << Verbose(1) << "There is a better candidate than " << *ThirdNode << " with " << ShortestAngle << " from baseline " << *Line << ": " << *OtherOptCandidate << " with " << OtherShortestAngle << "." << endl; 1871 result = true; 1872 break; 1873 } 1874 } 1875 } 1876 delete(OptCandidates); 1877 if (result) 1878 break; 1879 } else { 1880 Log() << Verbose(1) << "Circumcircle for base line " << *Line << " and base triangle " << T << " is too big!" << endl; 1881 } 1882 } else { 1883 eLog() << Verbose(2) << "Baseline is connected to two triangles already?" << endl; 1884 } 1885 } else { 1886 Log() << Verbose(2) << "No present baseline between " << BaseRay->endpoints[0] << " and candidate " << *ThirdNode << "." << endl; 1887 } 1888 } 1889 } else { 1890 eLog() << Verbose(1) << "Could not find the TesselPoint " << *ThirdNode << "." << endl; 1891 } 1892 1893 Log() << Verbose(1) << "End of HasOtherBaselineBetterCandidate" << endl; 1894 1895 return result; 1896 }; 1795 //bool Tesselation::HasOtherBaselineBetterCandidate(CandidateForTesselation &CandidateLine, const TesselPoint * const ThirdNode, double RADIUS, const LinkedCell * const LC) const 1796 //{ 1797 // Info FunctionInfo(__func__); 1798 // bool result = false; 1799 // Vector CircleCenter; 1800 // Vector CirclePlaneNormal; 1801 // Vector OldSphereCenter; 1802 // Vector SearchDirection; 1803 // Vector helper; 1804 // TesselPoint *OtherOptCandidate = NULL; 1805 // double OtherShortestAngle = 2.*M_PI; // This will indicate the quadrant. 1806 // double radius, CircleRadius; 1807 // BoundaryLineSet *Line = NULL; 1808 // BoundaryTriangleSet *T = NULL; 1809 // 1810 // // check both other lines 1811 // PointMap::const_iterator FindPoint = PointsOnBoundary.find(ThirdNode->nr); 1812 // if (FindPoint != PointsOnBoundary.end()) { 1813 // for (int i=0;i<2;i++) { 1814 // LineMap::const_iterator FindLine = (FindPoint->second)->lines.find(BaseRay->endpoints[0]->node->nr); 1815 // if (FindLine != (FindPoint->second)->lines.end()) { 1816 // Line = FindLine->second; 1817 // Log() << Verbose(0) << "Found line " << *Line << "." << endl; 1818 // if (Line->triangles.size() == 1) { 1819 // T = Line->triangles.begin()->second; 1820 // // construct center of circle 1821 // CircleCenter.CopyVector(Line->endpoints[0]->node->node); 1822 // CircleCenter.AddVector(Line->endpoints[1]->node->node); 1823 // CircleCenter.Scale(0.5); 1824 // 1825 // // construct normal vector of circle 1826 // CirclePlaneNormal.CopyVector(Line->endpoints[0]->node->node); 1827 // CirclePlaneNormal.SubtractVector(Line->endpoints[1]->node->node); 1828 // 1829 // // calculate squared radius of circle 1830 // radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal); 1831 // if (radius/4. < RADIUS*RADIUS) { 1832 // CircleRadius = RADIUS*RADIUS - radius/4.; 1833 // CirclePlaneNormal.Normalize(); 1834 // //Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl; 1835 // 1836 // // construct old center 1837 // GetCenterofCircumcircle(&OldSphereCenter, *T->endpoints[0]->node->node, *T->endpoints[1]->node->node, *T->endpoints[2]->node->node); 1838 // helper.CopyVector(&T->NormalVector); // normal vector ensures that this is correct center of the two possible ones 1839 // radius = Line->endpoints[0]->node->node->DistanceSquared(&OldSphereCenter); 1840 // helper.Scale(sqrt(RADIUS*RADIUS - radius)); 1841 // OldSphereCenter.AddVector(&helper); 1842 // OldSphereCenter.SubtractVector(&CircleCenter); 1843 // //Log() << Verbose(1) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl; 1844 // 1845 // // construct SearchDirection 1846 // SearchDirection.MakeNormalVector(&T->NormalVector, &CirclePlaneNormal); 1847 // helper.CopyVector(Line->endpoints[0]->node->node); 1848 // helper.SubtractVector(ThirdNode->node); 1849 // if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards! 1850 // SearchDirection.Scale(-1.); 1851 // SearchDirection.ProjectOntoPlane(&OldSphereCenter); 1852 // SearchDirection.Normalize(); 1853 // Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl; 1854 // if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) { 1855 // // rotated the wrong way! 1856 // eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl; 1857 // } 1858 // 1859 // // add third point 1860 // FindThirdPointForTesselation(T->NormalVector, SearchDirection, OldSphereCenter, OptCandidates, ThirdNode, RADIUS, LC); 1861 // for (TesselPointList::iterator it = OptCandidates.pointlist.begin(); it != OptCandidates.pointlist.end(); ++it) { 1862 // if (((*it) == BaseRay->endpoints[0]->node) || ((*it) == BaseRay->endpoints[1]->node)) // skip if it's the same triangle than suggested 1863 // continue; 1864 // Log() << Verbose(0) << " Third point candidate is " << (*it) 1865 // << " with circumsphere's center at " << (*it)->OptCenter << "." << endl; 1866 // Log() << Verbose(0) << " Baseline is " << *BaseRay << endl; 1867 // 1868 // // check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2) 1869 // TesselPoint *PointCandidates[3]; 1870 // PointCandidates[0] = (*it); 1871 // PointCandidates[1] = BaseRay->endpoints[0]->node; 1872 // PointCandidates[2] = BaseRay->endpoints[1]->node; 1873 // bool check=false; 1874 // int existentTrianglesCount = CheckPresenceOfTriangle(PointCandidates); 1875 // // If there is no triangle, add it regularly. 1876 // if (existentTrianglesCount == 0) { 1877 // SetTesselationPoint((*it), 0); 1878 // SetTesselationPoint(BaseRay->endpoints[0]->node, 1); 1879 // SetTesselationPoint(BaseRay->endpoints[1]->node, 2); 1880 // 1881 // if (CheckLineCriteriaForDegeneratedTriangle((const BoundaryPointSet ** const )TPS)) { 1882 // OtherOptCandidate = (*it); 1883 // check = true; 1884 // } 1885 // } else if ((existentTrianglesCount >= 1) && (existentTrianglesCount <= 3)) { // If there is a planar region within the structure, we need this triangle a second time. 1886 // SetTesselationPoint((*it), 0); 1887 // SetTesselationPoint(BaseRay->endpoints[0]->node, 1); 1888 // SetTesselationPoint(BaseRay->endpoints[1]->node, 2); 1889 // 1890 // // We demand that at most one new degenerate line is created and that this line also already exists (which has to be the case due to existentTrianglesCount == 1) 1891 // // i.e. at least one of the three lines must be present with TriangleCount <= 1 1892 // if (CheckLineCriteriaForDegeneratedTriangle((const BoundaryPointSet ** const)TPS)) { 1893 // OtherOptCandidate = (*it); 1894 // check = true; 1895 // } 1896 // } 1897 // 1898 // if (check) { 1899 // if (ShortestAngle > OtherShortestAngle) { 1900 // Log() << Verbose(0) << "There is a better candidate than " << *ThirdNode << " with " << ShortestAngle << " from baseline " << *Line << ": " << *OtherOptCandidate << " with " << OtherShortestAngle << "." << endl; 1901 // result = true; 1902 // break; 1903 // } 1904 // } 1905 // } 1906 // delete(OptCandidates); 1907 // if (result) 1908 // break; 1909 // } else { 1910 // Log() << Verbose(0) << "Circumcircle for base line " << *Line << " and base triangle " << T << " is too big!" << endl; 1911 // } 1912 // } else { 1913 // eLog() << Verbose(2) << "Baseline is connected to two triangles already?" << endl; 1914 // } 1915 // } else { 1916 // Log() << Verbose(1) << "No present baseline between " << BaseRay->endpoints[0] << " and candidate " << *ThirdNode << "." << endl; 1917 // } 1918 // } 1919 // } else { 1920 // eLog() << Verbose(1) << "Could not find the TesselPoint " << *ThirdNode << "." << endl; 1921 // } 1922 // 1923 // return result; 1924 //}; 1897 1925 1898 1926 /** This function finds a triangle to a line, adjacent to an existing one. … … 1906 1934 bool Tesselation::FindNextSuitableTriangle(CandidateForTesselation &CandidateLine, BoundaryTriangleSet &T, const double& RADIUS, const LinkedCell *LC) 1907 1935 { 1908 Log() << Verbose(0) << "Begin of FindNextSuitableTriangle\n";1936 Info FunctionInfo(__func__); 1909 1937 bool result = true; 1910 CandidateList *OptCandidates = new CandidateList();1911 1938 1912 1939 Vector CircleCenter; … … 1917 1944 TesselPoint *ThirdNode = NULL; 1918 1945 LineMap::iterator testline; 1919 double ShortestAngle = 2.*M_PI; // This will indicate the quadrant.1920 1946 double radius, CircleRadius; 1921 1947 1922 Log() << Verbose( 1) << "Current baseline is " << *CandidateLine.BaseLine << " of triangle " << T << "." << endl;1948 Log() << Verbose(0) << "Current baseline is " << *CandidateLine.BaseLine << " of triangle " << T << "." << endl; 1923 1949 for (int i=0;i<3;i++) 1924 1950 if ((T.endpoints[i]->node != CandidateLine.BaseLine->endpoints[0]->node) && (T.endpoints[i]->node != CandidateLine.BaseLine->endpoints[1]->node)) … … 1939 1965 CircleRadius = RADIUS*RADIUS - radius/4.; 1940 1966 CirclePlaneNormal.Normalize(); 1941 //Log() << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;1967 Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl; 1942 1968 1943 1969 // construct old center … … 1948 1974 OldSphereCenter.AddVector(&helper); 1949 1975 OldSphereCenter.SubtractVector(&CircleCenter); 1950 //Log() << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;1976 Log() << Verbose(1) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl; 1951 1977 1952 1978 // construct SearchDirection … … 1958 1984 SearchDirection.ProjectOntoPlane(&OldSphereCenter); 1959 1985 SearchDirection.Normalize(); 1960 Log() << Verbose( 2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;1986 Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl; 1961 1987 if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) { 1962 1988 // rotated the wrong way! … … 1965 1991 1966 1992 // add third point 1967 FindThirdPointForTesselation(T.NormalVector, SearchDirection, OldSphereCenter, CandidateLine .BaseLine, ThirdNode, OptCandidates, &ShortestAngle, RADIUS, LC);1993 FindThirdPointForTesselation(T.NormalVector, SearchDirection, OldSphereCenter, CandidateLine, ThirdNode, RADIUS, LC); 1968 1994 1969 1995 } else { 1970 Log() << Verbose( 1) << "Circumcircle for base line " << *CandidateLine.BaseLine << " and base triangle " << T << " is too big!" << endl;1971 } 1972 1973 if ( OptCandidates->begin() == OptCandidates->end()) {1996 Log() << Verbose(0) << "Circumcircle for base line " << *CandidateLine.BaseLine << " and base triangle " << T << " is too big!" << endl; 1997 } 1998 1999 if (CandidateLine.pointlist.empty()) { 1974 2000 eLog() << Verbose(2) << "Could not find a suitable candidate." << endl; 1975 2001 return false; 1976 2002 } 1977 Log() << Verbose(1) << "Third Points are "; 1978 for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) { 1979 Log() << Verbose(1) << " " << *(*it)->point << endl; 1980 } 1981 1982 BoundaryLineSet *BaseRay = CandidateLine.BaseLine; 1983 for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) { 1984 Log() << Verbose(1) << " Third point candidate is " << *(*it)->point 1985 << " with circumsphere's center at " << (*it)->OptCenter << "." << endl; 1986 Log() << Verbose(1) << " Baseline is " << *BaseRay << endl; 1987 1988 // check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2) 1989 TesselPoint *PointCandidates[3]; 1990 PointCandidates[0] = (*it)->point; 1991 PointCandidates[1] = BaseRay->endpoints[0]->node; 1992 PointCandidates[2] = BaseRay->endpoints[1]->node; 1993 int existentTrianglesCount = CheckPresenceOfTriangle(PointCandidates); 1994 1995 BTS = NULL; 1996 // check for present edges and whether we reach better candidates from them 1997 //if (HasOtherBaselineBetterCandidate(BaseRay, (*it)->point, ShortestAngle, RADIUS, LC) ) { 1998 if (0) { 1999 result = false; 2000 break; 2001 } else { 2002 // If there is no triangle, add it regularly. 2003 if (existentTrianglesCount == 0) { 2004 AddTesselationPoint((*it)->point, 0); 2005 AddTesselationPoint(BaseRay->endpoints[0]->node, 1); 2006 AddTesselationPoint(BaseRay->endpoints[1]->node, 2); 2007 2008 if (CheckLineCriteriaForDegeneratedTriangle((const BoundaryPointSet ** const )TPS)) { 2009 CandidateLine.point = (*it)->point; 2010 CandidateLine.OptCenter.CopyVector(&((*it)->OptCenter)); 2011 CandidateLine.OtherOptCenter.CopyVector(&((*it)->OtherOptCenter)); 2012 CandidateLine.ShortestAngle = ShortestAngle; 2013 } else { 2014 // eLog() << Verbose(2) << "This triangle consisting of "; 2015 // Log() << Verbose(0) << *(*it)->point << ", "; 2016 // Log() << Verbose(0) << *BaseRay->endpoints[0]->node << " and "; 2017 // Log() << Verbose(0) << *BaseRay->endpoints[1]->node << " "; 2018 // Log() << Verbose(0) << "exists and is not added, as it 0x80000000006fc150(does not seem helpful!" << endl; 2019 result = false; 2020 } 2021 } else if ((existentTrianglesCount >= 1) && (existentTrianglesCount <= 3)) { // If there is a planar region within the structure, we need this triangle a second time. 2022 AddTesselationPoint((*it)->point, 0); 2023 AddTesselationPoint(BaseRay->endpoints[0]->node, 1); 2024 AddTesselationPoint(BaseRay->endpoints[1]->node, 2); 2025 2026 // We demand that at most one new degenerate line is created and that this line also already exists (which has to be the case due to existentTrianglesCount == 1) 2027 // i.e. at least one of the three lines must be present with TriangleCount <= 1 2028 if (CheckLineCriteriaForDegeneratedTriangle((const BoundaryPointSet ** const)TPS) || CandidateLine.BaseLine->skipped) { 2029 CandidateLine.point = (*it)->point; 2030 CandidateLine.OptCenter.CopyVector(&(*it)->OptCenter); 2031 CandidateLine.OtherOptCenter.CopyVector(&(*it)->OtherOptCenter); 2032 CandidateLine.ShortestAngle = ShortestAngle+2.*M_PI; 2033 2034 } else { 2035 // eLog() << Verbose(2) << "This triangle consisting of " << *(*it)->point << ", " << *BaseRay->endpoints[0]->node << " and " << *BaseRay->endpoints[1]->node << " " << "exists and is not added, as it does not seem helpful!" << endl; 2036 result = false; 2037 } 2038 } else { 2039 // Log() << Verbose(1) << "This triangle consisting of "; 2040 // Log() << Verbose(0) << *(*it)->point << ", "; 2041 // Log() << Verbose(0) << *BaseRay->endpoints[0]->node << " and "; 2042 // Log() << Verbose(0) << *BaseRay->endpoints[1]->node << " "; 2043 // Log() << Verbose(0) << "is invalid!" << endl; 2044 result = false; 2045 } 2046 } 2047 2048 // set baseline to new ray from ref point (here endpoints[0]->node) to current candidate (here (*it)->point)) 2049 BaseRay = BLS[0]; 2050 if ((BTS != NULL) && (BTS->NormalVector.NormSquared() < MYEPSILON)) { 2051 eLog() << Verbose(1) << "Triangle " << *BTS << " has zero normal vector!" << endl; 2052 exit(255); 2053 } 2054 2055 } 2056 2057 // remove all candidates from the list and then the list itself 2058 class CandidateForTesselation *remover = NULL; 2059 for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) { 2060 remover = *it; 2061 delete(remover); 2062 } 2063 delete(OptCandidates); 2064 Log() << Verbose(0) << "End of FindNextSuitableTriangle\n"; 2003 Log() << Verbose(0) << "Third Points are: " << endl; 2004 for (TesselPointList::iterator it = CandidateLine.pointlist.begin(); it != CandidateLine.pointlist.end(); ++it) { 2005 Log() << Verbose(0) << " " << *(*it) << endl; 2006 } 2007 2008 return true; 2009 2010 // BoundaryLineSet *BaseRay = CandidateLine.BaseLine; 2011 // for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) { 2012 // Log() << Verbose(0) << "Third point candidate is " << *(*it)->point 2013 // << " with circumsphere's center at " << (*it)->OptCenter << "." << endl; 2014 // Log() << Verbose(0) << "Baseline is " << *BaseRay << endl; 2015 // 2016 // // check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2) 2017 // TesselPoint *PointCandidates[3]; 2018 // PointCandidates[0] = (*it)->point; 2019 // PointCandidates[1] = BaseRay->endpoints[0]->node; 2020 // PointCandidates[2] = BaseRay->endpoints[1]->node; 2021 // int existentTrianglesCount = CheckPresenceOfTriangle(PointCandidates); 2022 // 2023 // BTS = NULL; 2024 // // check for present edges and whether we reach better candidates from them 2025 // //if (HasOtherBaselineBetterCandidate(BaseRay, (*it)->point, ShortestAngle, RADIUS, LC) ) { 2026 // if (0) { 2027 // result = false; 2028 // break; 2029 // } else { 2030 // // If there is no triangle, add it regularly. 2031 // if (existentTrianglesCount == 0) { 2032 // AddTesselationPoint((*it)->point, 0); 2033 // AddTesselationPoint(BaseRay->endpoints[0]->node, 1); 2034 // AddTesselationPoint(BaseRay->endpoints[1]->node, 2); 2035 // 2036 // if (CheckLineCriteriaForDegeneratedTriangle((const BoundaryPointSet ** const )TPS)) { 2037 // CandidateLine.point = (*it)->point; 2038 // CandidateLine.OptCenter.CopyVector(&((*it)->OptCenter)); 2039 // CandidateLine.OtherOptCenter.CopyVector(&((*it)->OtherOptCenter)); 2040 // CandidateLine.ShortestAngle = ShortestAngle; 2041 // } else { 2042 //// eLog() << Verbose(1) << "This triangle consisting of "; 2043 //// Log() << Verbose(0) << *(*it)->point << ", "; 2044 //// Log() << Verbose(0) << *BaseRay->endpoints[0]->node << " and "; 2045 //// Log() << Verbose(0) << *BaseRay->endpoints[1]->node << " "; 2046 //// Log() << Verbose(0) << "exists and is not added, as it 0x80000000006fc150(does not seem helpful!" << endl; 2047 // result = false; 2048 // } 2049 // } else if ((existentTrianglesCount >= 1) && (existentTrianglesCount <= 3)) { // If there is a planar region within the structure, we need this triangle a second time. 2050 // AddTesselationPoint((*it)->point, 0); 2051 // AddTesselationPoint(BaseRay->endpoints[0]->node, 1); 2052 // AddTesselationPoint(BaseRay->endpoints[1]->node, 2); 2053 // 2054 // // We demand that at most one new degenerate line is created and that this line also already exists (which has to be the case due to existentTrianglesCount == 1) 2055 // // i.e. at least one of the three lines must be present with TriangleCount <= 1 2056 // if (CheckLineCriteriaForDegeneratedTriangle((const BoundaryPointSet ** const)TPS) || CandidateLine.BaseLine->skipped) { 2057 // CandidateLine.point = (*it)->point; 2058 // CandidateLine.OptCenter.CopyVector(&(*it)->OptCenter); 2059 // CandidateLine.OtherOptCenter.CopyVector(&(*it)->OtherOptCenter); 2060 // CandidateLine.ShortestAngle = ShortestAngle+2.*M_PI; 2061 // 2062 // } else { 2063 //// eLog() << Verbose(1) << "This triangle consisting of " << *(*it)->point << ", " << *BaseRay->endpoints[0]->node << " and " << *BaseRay->endpoints[1]->node << " " << "exists and is not added, as it does not seem helpful!" << endl; 2064 // result = false; 2065 // } 2066 // } else { 2067 //// Log() << Verbose(1) << "This triangle consisting of "; 2068 //// Log() << Verbose(0) << *(*it)->point << ", "; 2069 //// Log() << Verbose(0) << *BaseRay->endpoints[0]->node << " and "; 2070 //// Log() << Verbose(0) << *BaseRay->endpoints[1]->node << " "; 2071 //// Log() << Verbose(0) << "is invalid!" << endl; 2072 // result = false; 2073 // } 2074 // } 2075 // 2076 // // set baseline to new ray from ref point (here endpoints[0]->node) to current candidate (here (*it)->point)) 2077 // BaseRay = BLS[0]; 2078 // if ((BTS != NULL) && (BTS->NormalVector.NormSquared() < MYEPSILON)) { 2079 // eLog() << Verbose(1) << "Triangle " << *BTS << " has zero normal vector!" << endl; 2080 // exit(255); 2081 // } 2082 // 2083 // } 2084 // 2085 // // remove all candidates from the list and then the list itself 2086 // class CandidateForTesselation *remover = NULL; 2087 // for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) { 2088 // remover = *it; 2089 // delete(remover); 2090 // } 2091 // delete(OptCandidates); 2065 2092 return result; 2066 2093 }; 2067 2094 2068 2095 /** Adds the present line and candidate point from \a &CandidateLine to the Tesselation. 2069 * \param &CandidateLine triangle to add 2070 */ 2071 void Tesselation::AddCandidateTriangle(CandidateForTesselation &CandidateLine) 2072 { 2096 * \param CandidateLine triangle to add 2097 * \NOTE we need the copy operator here as the original CandidateForTesselation is removed in AddTesselationLine() 2098 */ 2099 void Tesselation::AddCandidateTriangle(CandidateForTesselation CandidateLine) 2100 { 2101 Info FunctionInfo(__func__); 2073 2102 Vector Center; 2074 Log() << Verbose(2) << "BaseLine is :" << *(CandidateLine.BaseLine) << " with candidate " << *(CandidateLine.point) << "." << endl; 2075 // add the points 2076 AddTesselationPoint(CandidateLine.point, 0); 2077 AddTesselationPoint(CandidateLine.BaseLine->endpoints[0]->node, 1); 2078 AddTesselationPoint(CandidateLine.BaseLine->endpoints[1]->node, 2); 2079 2080 Center.CopyVector(&CandidateLine.OptCenter); 2081 // add the lines 2082 AddTesselationLine(TPS[0], TPS[1], 0); 2083 AddTesselationLine(TPS[0], TPS[2], 1); 2084 AddTesselationLine(TPS[1], TPS[2], 2); 2085 2086 // add the triangles 2087 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount); 2088 AddTesselationTriangle(); 2089 Center.Scale(-1.); 2090 BTS->GetNormalVector(Center); 2091 2092 Log() << Verbose(0) << "--> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " for this triangle ... " << endl; 2103 BoundaryLineSet *BaseLine = CandidateLine.BaseLine; 2104 2105 // go through all candidates (in degenerate n-nodes case we may have to add multiple triangles) 2106 for (TesselPointList::iterator Runner = CandidateLine.pointlist.begin(); Runner != CandidateLine.pointlist.end(); Runner++) { 2107 Log() << Verbose(2) << "BaseLine is :" << *(BaseLine) << " with current candidate " << *(*Runner) << "." << endl; 2108 // add the points 2109 AddTesselationPoint((*Runner), 0); 2110 AddTesselationPoint(BaseLine->endpoints[0]->node, 1); 2111 AddTesselationPoint(BaseLine->endpoints[1]->node, 2); 2112 2113 Center.CopyVector(&CandidateLine.OptCenter); 2114 // add the lines 2115 AddTesselationLine(TPS[0], TPS[1], 0); 2116 AddTesselationLine(TPS[0], TPS[2], 1); 2117 AddTesselationLine(TPS[1], TPS[2], 2); 2118 2119 // add the triangles 2120 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount); 2121 AddTesselationTriangle(); 2122 Center.Scale(-1.); 2123 BTS->GetNormalVector(Center); 2124 2125 Log() << Verbose(0) << "--> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector << "." << endl; 2126 BaseLine = BLS[0]; // go to next baseline, in order to add in star formation from CandidateLine.BaseLine->endpoints[0]->node 2127 } 2093 2128 }; 2094 2129 … … 2102 2137 class BoundaryPointSet *Tesselation::IsConvexRectangle(class BoundaryLineSet *Base) 2103 2138 { 2139 Info FunctionInfo(__func__); 2104 2140 class BoundaryPointSet *Spot = NULL; 2105 2141 class BoundaryLineSet *OtherBase; … … 2113 2149 OtherBase = new class BoundaryLineSet(BPS,-1); 2114 2150 2115 Log() << Verbose( 3) << "INFO: Current base line is " << *Base << "." << endl;2116 Log() << Verbose( 3) << "INFO: Other base line is " << *OtherBase << "." << endl;2151 Log() << Verbose(1) << "INFO: Current base line is " << *Base << "." << endl; 2152 Log() << Verbose(1) << "INFO: Other base line is " << *OtherBase << "." << endl; 2117 2153 2118 2154 // get the closest point on each line to the other line … … 2134 2170 delete(ClosestPoint); 2135 2171 if ((distance[0] * distance[1]) > 0) { // have same sign? 2136 Log() << Verbose( 3) << "REJECT: Both SKPs have same sign: " << distance[0] << " and " << distance[1] << ". " << *Base << "' rectangle is concave." << endl;2172 Log() << Verbose(1) << "REJECT: Both SKPs have same sign: " << distance[0] << " and " << distance[1] << ". " << *Base << "' rectangle is concave." << endl; 2137 2173 if (distance[0] < distance[1]) { 2138 2174 Spot = Base->endpoints[0]; … … 2142 2178 return Spot; 2143 2179 } else { // different sign, i.e. we are in between 2144 Log() << Verbose( 3) << "ACCEPT: Rectangle of triangles of base line " << *Base << " is convex." << endl;2180 Log() << Verbose(0) << "ACCEPT: Rectangle of triangles of base line " << *Base << " is convex." << endl; 2145 2181 return NULL; 2146 2182 } … … 2150 2186 void Tesselation::PrintAllBoundaryPoints(ofstream *out) const 2151 2187 { 2188 Info FunctionInfo(__func__); 2152 2189 // print all lines 2153 Log() << Verbose( 1) << "Printing all boundary points for debugging:" << endl;2190 Log() << Verbose(0) << "Printing all boundary points for debugging:" << endl; 2154 2191 for (PointMap::const_iterator PointRunner = PointsOnBoundary.begin();PointRunner != PointsOnBoundary.end(); PointRunner++) 2155 Log() << Verbose( 2) << *(PointRunner->second) << endl;2192 Log() << Verbose(0) << *(PointRunner->second) << endl; 2156 2193 }; 2157 2194 2158 2195 void Tesselation::PrintAllBoundaryLines(ofstream *out) const 2159 2196 { 2197 Info FunctionInfo(__func__); 2160 2198 // print all lines 2161 Log() << Verbose( 1) << "Printing all boundary lines for debugging:" << endl;2199 Log() << Verbose(0) << "Printing all boundary lines for debugging:" << endl; 2162 2200 for (LineMap::const_iterator LineRunner = LinesOnBoundary.begin(); LineRunner != LinesOnBoundary.end(); LineRunner++) 2163 Log() << Verbose( 2) << *(LineRunner->second) << endl;2201 Log() << Verbose(0) << *(LineRunner->second) << endl; 2164 2202 }; 2165 2203 2166 2204 void Tesselation::PrintAllBoundaryTriangles(ofstream *out) const 2167 2205 { 2206 Info FunctionInfo(__func__); 2168 2207 // print all triangles 2169 Log() << Verbose( 1) << "Printing all boundary triangles for debugging:" << endl;2208 Log() << Verbose(0) << "Printing all boundary triangles for debugging:" << endl; 2170 2209 for (TriangleMap::const_iterator TriangleRunner = TrianglesOnBoundary.begin(); TriangleRunner != TrianglesOnBoundary.end(); TriangleRunner++) 2171 Log() << Verbose( 2) << *(TriangleRunner->second) << endl;2210 Log() << Verbose(0) << *(TriangleRunner->second) << endl; 2172 2211 }; 2173 2212 … … 2179 2218 double Tesselation::PickFarthestofTwoBaselines(class BoundaryLineSet *Base) 2180 2219 { 2220 Info FunctionInfo(__func__); 2181 2221 class BoundaryLineSet *OtherBase; 2182 2222 Vector *ClosestPoint[2]; … … 2190 2230 OtherBase = new class BoundaryLineSet(BPS,-1); 2191 2231 2192 Log() << Verbose( 3) << "INFO: Current base line is " << *Base << "." << endl;2193 Log() << Verbose( 3) << "INFO: Other base line is " << *OtherBase << "." << endl;2232 Log() << Verbose(0) << "INFO: Current base line is " << *Base << "." << endl; 2233 Log() << Verbose(0) << "INFO: Other base line is " << *OtherBase << "." << endl; 2194 2234 2195 2235 // get the closest point on each line to the other line … … 2211 2251 2212 2252 if (Distance.NormSquared() < MYEPSILON) { // check for intersection 2213 Log() << Verbose( 3) << "REJECT: Both lines have an intersection: Nothing to do." << endl;2253 Log() << Verbose(0) << "REJECT: Both lines have an intersection: Nothing to do." << endl; 2214 2254 return false; 2215 2255 } else { // check for sign against BaseLineNormal … … 2221 2261 } 2222 2262 for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) { 2223 Log() << Verbose( 4) << "INFO: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << "." << endl;2263 Log() << Verbose(1) << "INFO: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << "." << endl; 2224 2264 BaseLineNormal.AddVector(&(runner->second->NormalVector)); 2225 2265 } … … 2227 2267 2228 2268 if (Distance.ScalarProduct(&BaseLineNormal) > MYEPSILON) { // Distance points outwards, hence OtherBase higher than Base -> flip 2229 Log() << Verbose( 2) << "ACCEPT: Other base line would be higher: Flipping baseline." << endl;2269 Log() << Verbose(0) << "ACCEPT: Other base line would be higher: Flipping baseline." << endl; 2230 2270 // calculate volume summand as a general tetraeder 2231 2271 return volume; 2232 2272 } else { // Base higher than OtherBase -> do nothing 2233 Log() << Verbose( 2) << "REJECT: Base line is higher: Nothing to do." << endl;2273 Log() << Verbose(0) << "REJECT: Base line is higher: Nothing to do." << endl; 2234 2274 return 0.; 2235 2275 } … … 2246 2286 class BoundaryLineSet * Tesselation::FlipBaseline(class BoundaryLineSet *Base) 2247 2287 { 2288 Info FunctionInfo(__func__); 2248 2289 class BoundaryLineSet *OldLines[4], *NewLine; 2249 2290 class BoundaryPointSet *OldPoints[2]; … … 2252 2293 int i,m; 2253 2294 2254 Log() << Verbose(1) << "Begin of FlipBaseline" << endl;2255 2256 2295 // calculate NormalVector for later use 2257 2296 BaseLineNormal.Zero(); … … 2261 2300 } 2262 2301 for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) { 2263 Log() << Verbose( 4) << "INFO: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << "." << endl;2302 Log() << Verbose(1) << "INFO: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << "." << endl; 2264 2303 BaseLineNormal.AddVector(&(runner->second->NormalVector)); 2265 2304 } … … 2274 2313 i=0; 2275 2314 m=0; 2276 Log() << Verbose( 3) << "The four old lines are: ";2315 Log() << Verbose(0) << "The four old lines are: "; 2277 2316 for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) 2278 2317 for (int j=0;j<3;j++) // all of their endpoints and baselines … … 2282 2321 } 2283 2322 Log() << Verbose(0) << endl; 2284 Log() << Verbose( 3) << "The two old points are: ";2323 Log() << Verbose(0) << "The two old points are: "; 2285 2324 for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) 2286 2325 for (int j=0;j<3;j++) // all of their endpoints and baselines … … 2308 2347 2309 2348 // remove triangles and baseline removes itself 2310 Log() << Verbose( 3) << "INFO: Deleting baseline " << *Base << " from global list." << endl;2349 Log() << Verbose(0) << "INFO: Deleting baseline " << *Base << " from global list." << endl; 2311 2350 OldBaseLineNr = Base->Nr; 2312 2351 m=0; 2313 2352 for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) { 2314 Log() << Verbose( 3) << "INFO: Deleting triangle " << *(runner->second) << "." << endl;2353 Log() << Verbose(0) << "INFO: Deleting triangle " << *(runner->second) << "." << endl; 2315 2354 OldTriangleNrs[m++] = runner->second->Nr; 2316 2355 RemoveTesselationTriangle(runner->second); … … 2322 2361 NewLine = new class BoundaryLineSet(BPS, OldBaseLineNr); 2323 2362 LinesOnBoundary.insert(LinePair(OldBaseLineNr, NewLine)); // no need for check for unique insertion as NewLine is definitely a new one 2324 Log() << Verbose( 3) << "INFO: Created new baseline " << *NewLine << "." << endl;2363 Log() << Verbose(0) << "INFO: Created new baseline " << *NewLine << "." << endl; 2325 2364 2326 2365 // construct new triangles with flipped baseline … … 2337 2376 BTS->GetNormalVector(BaseLineNormal); 2338 2377 AddTesselationTriangle(OldTriangleNrs[0]); 2339 Log() << Verbose( 3) << "INFO: Created new triangle " << *BTS << "." << endl;2378 Log() << Verbose(0) << "INFO: Created new triangle " << *BTS << "." << endl; 2340 2379 2341 2380 BLS[0] = (i==2 ? OldLines[3] : OldLines[2]); … … 2345 2384 BTS->GetNormalVector(BaseLineNormal); 2346 2385 AddTesselationTriangle(OldTriangleNrs[1]); 2347 Log() << Verbose( 3) << "INFO: Created new triangle " << *BTS << "." << endl;2386 Log() << Verbose(0) << "INFO: Created new triangle " << *BTS << "." << endl; 2348 2387 } else { 2349 Log() << Verbose(1) << "The four old lines do not connect, something's utterly wrong here!" << endl;2388 eLog() << Verbose(0) << "The four old lines do not connect, something's utterly wrong here!" << endl; 2350 2389 return NULL; 2351 2390 } 2352 2391 2353 Log() << Verbose(1) << "End of FlipBaseline" << endl;2354 2392 return NewLine; 2355 2393 }; … … 2366 2404 void Tesselation::FindSecondPointForTesselation(TesselPoint* a, Vector Oben, TesselPoint*& OptCandidate, double Storage[3], double RADIUS, const LinkedCell *LC) 2367 2405 { 2368 Log() << Verbose(2) << "Begin of FindSecondPointForTesselation" << endl;2406 Info FunctionInfo(__func__); 2369 2407 Vector AngleCheck; 2370 2408 class TesselPoint* Candidate = NULL; … … 2387 2425 Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1; 2388 2426 } 2389 Log() << Verbose( 3) << "LC Intervals from [" << N[0] << "<->" << LC->N[0] << ", " << N[1] << "<->" << LC->N[1] << ", " << N[2] << "<->" << LC->N[2] << "] :"2427 Log() << Verbose(0) << "LC Intervals from [" << N[0] << "<->" << LC->N[0] << ", " << N[1] << "<->" << LC->N[1] << ", " << N[2] << "<->" << LC->N[2] << "] :" 2390 2428 << " [" << Nlower[0] << "," << Nupper[0] << "], " << " [" << Nlower[1] << "," << Nupper[1] << "], " << " [" << Nlower[2] << "," << Nupper[2] << "], " << endl; 2391 2429 … … 2394 2432 for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) { 2395 2433 const LinkedNodes *List = LC->GetCurrentCell(); 2396 //Log() << Verbose( 2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;2434 //Log() << Verbose(1) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl; 2397 2435 if (List != NULL) { 2398 2436 for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) { … … 2425 2463 angle = AngleCheck.Angle(&Oben); 2426 2464 if (angle < Storage[0]) { 2427 //Log() << Verbose( 3) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);2428 Log() << Verbose( 3) << "Current candidate is " << *Candidate << ": Is a better candidate with distance " << norm << " and angle " << angle << " to oben " << Oben << ".\n";2465 //Log() << Verbose(1) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]); 2466 Log() << Verbose(1) << "Current candidate is " << *Candidate << ": Is a better candidate with distance " << norm << " and angle " << angle << " to oben " << Oben << ".\n"; 2429 2467 OptCandidate = Candidate; 2430 2468 Storage[0] = angle; 2431 //Log() << Verbose( 3) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[2]);2469 //Log() << Verbose(1) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[2]); 2432 2470 } else { 2433 //Log() << Verbose( 3) << "Current candidate is " << *Candidate << ": Looses with angle " << angle << " to a better candidate " << *OptCandidate << endl;2471 //Log() << Verbose(1) << "Current candidate is " << *Candidate << ": Looses with angle " << angle << " to a better candidate " << *OptCandidate << endl; 2434 2472 } 2435 2473 } else { 2436 //Log() << Verbose( 3) << "Current candidate is " << *Candidate << ": Refused due to Radius " << norm << endl;2474 //Log() << Verbose(1) << "Current candidate is " << *Candidate << ": Refused due to Radius " << norm << endl; 2437 2475 } 2438 2476 } else { 2439 //Log() << Verbose( 3) << "Current candidate is " << *Candidate << ": Candidate is equal to first endpoint." << *a << "." << endl;2477 //Log() << Verbose(1) << "Current candidate is " << *Candidate << ": Candidate is equal to first endpoint." << *a << "." << endl; 2440 2478 } 2441 2479 } 2442 2480 } else { 2443 Log() << Verbose( 3) << "Linked cell list is empty." << endl;2481 Log() << Verbose(0) << "Linked cell list is empty." << endl; 2444 2482 } 2445 2483 } 2446 Log() << Verbose(2) << "End of FindSecondPointForTesselation" << endl;2447 2484 }; 2448 2485 … … 2473 2510 * @param SearchDirection general direction where to search for the next point, relative to center of BaseLine 2474 2511 * @param OldSphereCenter center of sphere for base triangle, relative to center of BaseLine, giving null angle for the parameter circle 2475 * @param BaseLine BoundaryLineSet with the current base line2512 * @param CandidateLine CandidateForTesselation with the current base line and list of candidates and ShortestAngle 2476 2513 * @param ThirdNode third point to avoid in search 2477 * @param candidates list of equally good candidates to return2478 * @param ShortestAngle the current path length on this circle band for the current OptCandidate2479 2514 * @param RADIUS radius of sphere 2480 2515 * @param *LC LinkedCell structure with neighbouring points 2481 2516 */ 2482 void Tesselation::FindThirdPointForTesselation(Vector &NormalVector, Vector &SearchDirection, Vector &OldSphereCenter, class BoundaryLineSet *BaseLine, const class TesselPoint * const ThirdNode, CandidateList* &candidates, double *ShortestAngle, const double RADIUS, const LinkedCell *LC) const 2483 { 2517 void Tesselation::FindThirdPointForTesselation(Vector &NormalVector, Vector &SearchDirection, Vector &OldSphereCenter, CandidateForTesselation &CandidateLine, const class TesselPoint * const ThirdNode, const double RADIUS, const LinkedCell *LC) const 2518 { 2519 Info FunctionInfo(__func__); 2484 2520 Vector CircleCenter; // center of the circle, i.e. of the band of sphere's centers 2485 2521 Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in … … 2494 2530 int N[NDIM], Nlower[NDIM], Nupper[NDIM]; 2495 2531 TesselPoint *Candidate = NULL; 2496 CandidateForTesselation *optCandidate = NULL; 2497 2498 Log() << Verbose(1) << "Begin of FindThirdPointForTesselation" << endl; 2499 2500 Log() << Verbose(2) << "INFO: NormalVector of BaseTriangle is " << NormalVector << "." << endl; 2532 2533 Log() << Verbose(1) << "INFO: NormalVector of BaseTriangle is " << NormalVector << "." << endl; 2501 2534 2502 2535 // construct center of circle 2503 CircleCenter.CopyVector( BaseLine->endpoints[0]->node->node);2504 CircleCenter.AddVector( BaseLine->endpoints[1]->node->node);2536 CircleCenter.CopyVector(CandidateLine.BaseLine->endpoints[0]->node->node); 2537 CircleCenter.AddVector(CandidateLine.BaseLine->endpoints[1]->node->node); 2505 2538 CircleCenter.Scale(0.5); 2506 2539 2507 2540 // construct normal vector of circle 2508 CirclePlaneNormal.CopyVector( BaseLine->endpoints[0]->node->node);2509 CirclePlaneNormal.SubtractVector( BaseLine->endpoints[1]->node->node);2541 CirclePlaneNormal.CopyVector(CandidateLine.BaseLine->endpoints[0]->node->node); 2542 CirclePlaneNormal.SubtractVector(CandidateLine.BaseLine->endpoints[1]->node->node); 2510 2543 2511 2544 // calculate squared radius TesselPoint *ThirdNode,f circle … … 2514 2547 CircleRadius = RADIUS*RADIUS - radius/4.; 2515 2548 CirclePlaneNormal.Normalize(); 2516 //Log() << Verbose( 2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;2549 //Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl; 2517 2550 2518 2551 // test whether old center is on the band's plane … … 2523 2556 radius = OldSphereCenter.ScalarProduct(&OldSphereCenter); 2524 2557 if (fabs(radius - CircleRadius) < HULLEPSILON) { 2525 //Log() << Verbose( 2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl;2558 //Log() << Verbose(1) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl; 2526 2559 2527 2560 // check SearchDirection 2528 //Log() << Verbose( 2) << "INFO: SearchDirection is " << SearchDirection << "." << endl;2561 //Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl; 2529 2562 if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) { // rotated the wrong way! 2530 2563 eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are not orthogonal!" << endl; … … 2535 2568 for(int i=0;i<NDIM;i++) // store indices of this cell 2536 2569 N[i] = LC->n[i]; 2537 //Log() << Verbose( 2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;2570 //Log() << Verbose(1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl; 2538 2571 } else { 2539 2572 eLog() << Verbose(1) << "Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl; … … 2541 2574 } 2542 2575 // then go through the current and all neighbouring cells and check the contained points for possible candidates 2543 //Log() << Verbose( 2) << "LC Intervals:";2576 //Log() << Verbose(1) << "LC Intervals:"; 2544 2577 for (int i=0;i<NDIM;i++) { 2545 2578 Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0; … … 2552 2585 for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) { 2553 2586 const LinkedNodes *List = LC->GetCurrentCell(); 2554 //Log() << Verbose( 2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;2587 //Log() << Verbose(1) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl; 2555 2588 if (List != NULL) { 2556 2589 for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) { … … 2558 2591 2559 2592 // check for three unique points 2560 //Log() << Verbose(2) << "INFO: Current Candidate is " << *Candidate << " at " << Candidate->node<< "." << endl;2561 if ((Candidate != BaseLine->endpoints[0]->node) && (Candidate !=BaseLine->endpoints[1]->node) ){2593 Log() << Verbose(2) << "INFO: Current Candidate is " << *Candidate << " at " << *(Candidate->node) << "." << endl; 2594 if ((Candidate != CandidateLine.BaseLine->endpoints[0]->node) && (Candidate != CandidateLine.BaseLine->endpoints[1]->node) ){ 2562 2595 2563 2596 // construct both new centers 2564 GetCenterofCircumcircle(&NewSphereCenter, * BaseLine->endpoints[0]->node->node, *BaseLine->endpoints[1]->node->node, *Candidate->node);2597 GetCenterofCircumcircle(&NewSphereCenter, *CandidateLine.BaseLine->endpoints[0]->node->node, *CandidateLine.BaseLine->endpoints[1]->node->node, *Candidate->node); 2565 2598 OtherNewSphereCenter.CopyVector(&NewSphereCenter); 2566 2599 2567 if ((NewNormalVector.MakeNormalVector( BaseLine->endpoints[0]->node->node,BaseLine->endpoints[1]->node->node, Candidate->node))2600 if ((NewNormalVector.MakeNormalVector(CandidateLine.BaseLine->endpoints[0]->node->node, CandidateLine.BaseLine->endpoints[1]->node->node, Candidate->node)) 2568 2601 && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON) 2569 2602 ) { 2570 2603 helper.CopyVector(&NewNormalVector); 2571 //Log() << Verbose(2) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl;2572 radius = BaseLine->endpoints[0]->node->node->DistanceSquared(&NewSphereCenter);2604 Log() << Verbose(1) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl; 2605 radius = CandidateLine.BaseLine->endpoints[0]->node->node->DistanceSquared(&NewSphereCenter); 2573 2606 if (radius < RADIUS*RADIUS) { 2574 2607 helper.Scale(sqrt(RADIUS*RADIUS - radius)); 2575 //Log() << Verbose(2) << "INFO: Distance of NewCircleCenter to NewSphereCenter is " << helper.Norm() << " with sphere radius " << RADIUS << "." << endl;2608 Log() << Verbose(2) << "INFO: Distance of NewCircleCenter to NewSphereCenter is " << helper.Norm() << " with sphere radius " << RADIUS << "." << endl; 2576 2609 NewSphereCenter.AddVector(&helper); 2577 2610 NewSphereCenter.SubtractVector(&CircleCenter); 2578 //Log() << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl;2611 Log() << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl; 2579 2612 2580 2613 // OtherNewSphereCenter is created by the same vector just in the other direction … … 2582 2615 OtherNewSphereCenter.AddVector(&helper); 2583 2616 OtherNewSphereCenter.SubtractVector(&CircleCenter); 2584 //Log() << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl;2617 Log() << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl; 2585 2618 2586 2619 alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, NormalVector, SearchDirection); … … 2589 2622 // if there is a better candidate, drop the current list and add the new candidate 2590 2623 // otherwise ignore the new candidate and keep the list 2591 if (*ShortestAngle > (alpha - HULLEPSILON)) { 2592 optCandidate = new CandidateForTesselation(Candidate, BaseLine, OptCandidateCenter, OtherOptCandidateCenter); 2624 if (CandidateLine.ShortestAngle > (alpha - HULLEPSILON)) { 2593 2625 if (fabs(alpha - Otheralpha) > MYEPSILON) { 2594 optCandidate->OptCenter.CopyVector(&NewSphereCenter);2595 optCandidate->OtherOptCenter.CopyVector(&OtherNewSphereCenter);2626 CandidateLine.OptCenter.CopyVector(&NewSphereCenter); 2627 CandidateLine.OtherOptCenter.CopyVector(&OtherNewSphereCenter); 2596 2628 } else { 2597 optCandidate->OptCenter.CopyVector(&OtherNewSphereCenter);2598 optCandidate->OtherOptCenter.CopyVector(&NewSphereCenter);2629 CandidateLine.OptCenter.CopyVector(&OtherNewSphereCenter); 2630 CandidateLine.OtherOptCenter.CopyVector(&NewSphereCenter); 2599 2631 } 2600 2632 // if there is an equal candidate, add it to the list without clearing the list 2601 if (( *ShortestAngle - HULLEPSILON) < alpha) {2602 candidates->push_back(optCandidate);2603 Log() << Verbose( 2) << "ACCEPT: We have found an equally good candidate: " << *(optCandidate->point) << " with "2604 << alpha << " and circumsphere's center at " << optCandidate->OptCenter << "." << endl;2633 if ((CandidateLine.ShortestAngle - HULLEPSILON) < alpha) { 2634 CandidateLine.pointlist.push_back(Candidate); 2635 Log() << Verbose(0) << "ACCEPT: We have found an equally good candidate: " << *(Candidate) << " with " 2636 << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << "." << endl; 2605 2637 } else { 2606 2638 // remove all candidates from the list and then the list itself 2607 class CandidateForTesselation *remover = NULL; 2608 for (CandidateList::iterator it = candidates->begin(); it != candidates->end(); ++it) { 2609 remover = *it; 2610 delete(remover); 2611 } 2612 candidates->clear(); 2613 candidates->push_back(optCandidate); 2614 Log() << Verbose(2) << "ACCEPT: We have found a better candidate: " << *(optCandidate->point) << " with " 2615 << alpha << " and circumsphere's center at " << optCandidate->OptCenter << "." << endl; 2639 CandidateLine.pointlist.clear(); 2640 CandidateLine.pointlist.push_back(Candidate); 2641 Log() << Verbose(0) << "ACCEPT: We have found a better candidate: " << *(Candidate) << " with " 2642 << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << "." << endl; 2616 2643 } 2617 *ShortestAngle = alpha;2618 //Log() << Verbose(2) << "INFO: There are " << candidates->size() << " candidates in the list now." << endl;2644 CandidateLine.ShortestAngle = alpha; 2645 Log() << Verbose(0) << "INFO: There are " << CandidateLine.pointlist.size() << " candidates in the list now." << endl; 2619 2646 } else { 2620 if (( optCandidate != NULL) && (optCandidate->point != NULL)) {2621 //Log() << Verbose(2) << "REJECT: Old candidate " << *(optCandidate->point) << " with " << *ShortestAngle << " is better than new one " << *Candidate << " with " << alpha << " ." << endl;2647 if ((Candidate != NULL) && (CandidateLine.pointlist.begin() != CandidateLine.pointlist.end())) { 2648 Log() << Verbose(1) << "REJECT: Old candidate " << *(Candidate) << " with " << CandidateLine.ShortestAngle << " is better than new one " << *Candidate << " with " << alpha << " ." << endl; 2622 2649 } else { 2623 //Log() << Verbose(2) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl;2650 Log() << Verbose(1) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl; 2624 2651 } 2625 2652 } 2626 2653 2627 2654 } else { 2628 //Log() << Verbose(2) << "REJECT: NewSphereCenter " << NewSphereCenter << " for " << *Candidate << " is too far away: " << radius << "." << endl;2655 Log() << Verbose(1) << "REJECT: NewSphereCenter " << NewSphereCenter << " for " << *Candidate << " is too far away: " << radius << "." << endl; 2629 2656 } 2630 2657 } else { 2631 //Log() << Verbose(2) << "REJECT: Three points from " << *BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl;2658 Log() << Verbose(1) << "REJECT: Three points from " << *CandidateLine.BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl; 2632 2659 } 2633 2660 } else { 2634 2661 if (ThirdNode != NULL) { 2635 //Log() << Verbose(2) << "REJECT: Base triangle " << *BaseLine << " and " << *ThirdNode << " contains Candidate " << *Candidate << "." << endl;2662 Log() << Verbose(1) << "REJECT: Base triangle " << *CandidateLine.BaseLine << " and " << *ThirdNode << " contains Candidate " << *Candidate << "." << endl; 2636 2663 } else { 2637 //Log() << Verbose(2) << "REJECT: Base triangle " << *BaseLine << " contains Candidate " << *Candidate << "." << endl;2664 Log() << Verbose(1) << "REJECT: Base triangle " << *CandidateLine.BaseLine << " contains Candidate " << *Candidate << "." << endl; 2638 2665 } 2639 2666 } … … 2646 2673 } else { 2647 2674 if (ThirdNode != NULL) 2648 Log() << Verbose( 2) << "Circumcircle for base line " << *BaseLine << " and third node " << *ThirdNode << " is too big!" << endl;2675 Log() << Verbose(1) << "Circumcircle for base line " << *CandidateLine.BaseLine << " and third node " << *ThirdNode << " is too big!" << endl; 2649 2676 else 2650 Log() << Verbose(2) << "Circumcircle for base line " << *BaseLine << " is too big!" << endl; 2651 } 2652 2653 //Log() << Verbose(2) << "INFO: Sorting candidate list ..." << endl; 2654 if (candidates->size() > 1) { 2655 candidates->unique(); 2656 candidates->sort(SortCandidates); 2657 } 2658 2659 Log() << Verbose(1) << "End of FindThirdPointForTesselation" << endl; 2677 Log() << Verbose(1) << "Circumcircle for base line " << *CandidateLine.BaseLine << " is too big!" << endl; 2678 } 2679 2680 Log() << Verbose(1) << "INFO: Sorting candidate list ..." << endl; 2681 if (CandidateLine.pointlist.size() > 1) { 2682 CandidateLine.pointlist.unique(); 2683 CandidateLine.pointlist.sort(); //SortCandidates); 2684 } 2660 2685 }; 2661 2686 … … 2667 2692 class BoundaryPointSet *Tesselation::GetCommonEndpoint(const BoundaryLineSet * line1, const BoundaryLineSet * line2) const 2668 2693 { 2694 Info FunctionInfo(__func__); 2669 2695 const BoundaryLineSet * lines[2] = { line1, line2 }; 2670 2696 class BoundaryPointSet *node = NULL; … … 2680 2706 { // if insertion fails, we have common endpoint 2681 2707 node = OrderTest.first->second; 2682 Log() << Verbose( 5) << "Common endpoint of lines " << *line12708 Log() << Verbose(1) << "Common endpoint of lines " << *line1 2683 2709 << " and " << *line2 << " is: " << *node << "." << endl; 2684 2710 j = 2; … … 2697 2723 list<BoundaryTriangleSet*> * Tesselation::FindClosestTrianglesToPoint(const Vector *x, const LinkedCell* LC) const 2698 2724 { 2725 Info FunctionInfo(__func__); 2699 2726 TesselPoint *trianglePoints[3]; 2700 2727 TesselPoint *SecondPoint = NULL; … … 2702 2729 2703 2730 if (LinesOnBoundary.empty()) { 2704 Log() << Verbose(0) << "Error: There is no tesselation structure to compare the point with, please create one first.";2731 eLog() << Verbose(1) << "Error: There is no tesselation structure to compare the point with, please create one first."; 2705 2732 return NULL; 2706 2733 } … … 2710 2737 // check whether closest point is "too close" :), then it's inside 2711 2738 if (trianglePoints[0] == NULL) { 2712 Log() << Verbose( 2) << "Is the only point, no one else is closeby." << endl;2739 Log() << Verbose(0) << "Is the only point, no one else is closeby." << endl; 2713 2740 return NULL; 2714 2741 } 2715 2742 if (trianglePoints[0]->node->DistanceSquared(x) < MYEPSILON) { 2716 Log() << Verbose( 3) << "Point is right on a tesselation point, no nearest triangle." << endl;2743 Log() << Verbose(1) << "Point is right on a tesselation point, no nearest triangle." << endl; 2717 2744 PointMap::const_iterator PointRunner = PointsOnBoundary.find(trianglePoints[0]->nr); 2718 2745 triangles = new list<BoundaryTriangleSet*>; … … 2746 2773 eLog() << Verbose(1) << "IsInnerPoint encounters serious error, point " << i << " not found." << endl; 2747 2774 } 2748 //Log() << Verbose( 2) << "List of triangle points:" << endl;2749 //Log() << Verbose( 3) << *trianglePoints[i] << endl;2775 //Log() << Verbose(1) << "List of triangle points:" << endl; 2776 //Log() << Verbose(2) << *trianglePoints[i] << endl; 2750 2777 } 2751 2778 2752 2779 triangles = FindTriangles(trianglePoints); 2753 Log() << Verbose( 2) << "List of possible triangles:" << endl;2780 Log() << Verbose(1) << "List of possible triangles:" << endl; 2754 2781 for(list<BoundaryTriangleSet*>::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++) 2755 Log() << Verbose( 3) << **Runner << endl;2782 Log() << Verbose(2) << **Runner << endl; 2756 2783 2757 2784 delete(connectedClosestPoints); 2758 2785 } else { 2759 2786 triangles = NULL; 2760 Log() << Verbose(1) << "There is no circle of connected points!" << endl;2787 eLog() << Verbose(2) << "There is no circle of connected points!" << endl; 2761 2788 } 2762 2789 } … … 2778 2805 class BoundaryTriangleSet * Tesselation::FindClosestTriangleToPoint(const Vector *x, const LinkedCell* LC) const 2779 2806 { 2807 Info FunctionInfo(__func__); 2780 2808 class BoundaryTriangleSet *result = NULL; 2781 2809 list<BoundaryTriangleSet*> *triangles = FindClosestTrianglesToPoint(x, LC); … … 2787 2815 if (triangles->size() == 1) { // there is no degenerate case 2788 2816 result = triangles->front(); 2789 Log() << Verbose( 2) << "Normal Vector of this triangle is " << result->NormalVector << "." << endl;2817 Log() << Verbose(1) << "Normal Vector of this triangle is " << result->NormalVector << "." << endl; 2790 2818 } else { 2791 2819 result = triangles->front(); 2792 2820 result->GetCenter(&Center); 2793 2821 Center.SubtractVector(x); 2794 Log() << Verbose( 2) << "Normal Vector of this front side is " << result->NormalVector << "." << endl;2822 Log() << Verbose(1) << "Normal Vector of this front side is " << result->NormalVector << "." << endl; 2795 2823 if (Center.ScalarProduct(&result->NormalVector) < 0) { 2796 2824 result = triangles->back(); 2797 Log() << Verbose( 2) << "Normal Vector of this back side is " << result->NormalVector << "." << endl;2825 Log() << Verbose(1) << "Normal Vector of this back side is " << result->NormalVector << "." << endl; 2798 2826 if (Center.ScalarProduct(&result->NormalVector) < 0) { 2799 2827 eLog() << Verbose(1) << "Front and back side yield NormalVector in wrong direction!" << endl; … … 2814 2842 bool Tesselation::IsInnerPoint(const Vector &Point, const LinkedCell* const LC) const 2815 2843 { 2844 Info FunctionInfo(__func__); 2816 2845 class BoundaryTriangleSet *result = FindClosestTriangleToPoint(&Point, LC); 2817 2846 Vector Center; … … 2823 2852 2824 2853 result->GetCenter(&Center); 2825 Log() << Verbose( 3) << "INFO: Central point of the triangle is " << Center << "." << endl;2854 Log() << Verbose(2) << "INFO: Central point of the triangle is " << Center << "." << endl; 2826 2855 Center.SubtractVector(&Point); 2827 Log() << Verbose( 3) << "INFO: Vector from center to point to test is " << Center << "." << endl;2856 Log() << Verbose(2) << "INFO: Vector from center to point to test is " << Center << "." << endl; 2828 2857 if (Center.ScalarProduct(&result->NormalVector) > -MYEPSILON) { 2829 2858 Log() << Verbose(1) << Point << " is an inner point." << endl; … … 2844 2873 bool Tesselation::IsInnerPoint(const TesselPoint * const Point, const LinkedCell* const LC) const 2845 2874 { 2875 Info FunctionInfo(__func__); 2846 2876 return IsInnerPoint(*(Point->node), LC); 2847 2877 } … … 2855 2885 set<TesselPoint*> * Tesselation::GetAllConnectedPoints(const TesselPoint* const Point) const 2856 2886 { 2887 Info FunctionInfo(__func__); 2857 2888 set<TesselPoint*> *connectedPoints = new set<TesselPoint*>; 2858 2889 class BoundaryPointSet *ReferencePoint = NULL; 2859 2890 TesselPoint* current; 2860 2891 bool takePoint = false; 2861 2862 Log() << Verbose(3) << "Begin of GetAllConnectedPoints" << endl;2863 2892 2864 2893 // find the respective boundary point … … 2867 2896 ReferencePoint = PointRunner->second; 2868 2897 } else { 2869 Log() << Verbose(2) << "GetAllConnectedPoints() could not find the BoundaryPoint belonging to " << *Point << "." << endl;2898 eLog() << Verbose(2) << "GetAllConnectedPoints() could not find the BoundaryPoint belonging to " << *Point << "." << endl; 2870 2899 ReferencePoint = NULL; 2871 2900 } … … 2891 2920 2892 2921 if (takePoint) { 2893 Log() << Verbose( 5) << "INFO: Endpoint " << *current << " of line " << *(findLines->second) << " is enlisted." << endl;2922 Log() << Verbose(1) << "INFO: Endpoint " << *current << " of line " << *(findLines->second) << " is enlisted." << endl; 2894 2923 connectedPoints->insert(current); 2895 2924 } … … 2903 2932 } 2904 2933 2905 Log() << Verbose(3) << "End of GetAllConnectedPoints" << endl;2906 2934 return connectedPoints; 2907 2935 }; … … 2921 2949 list<TesselPoint*> * Tesselation::GetCircleOfConnectedPoints(const TesselPoint* const Point, const Vector * const Reference) const 2922 2950 { 2951 Info FunctionInfo(__func__); 2923 2952 map<double, TesselPoint*> anglesOfPoints; 2924 2953 set<TesselPoint*> *connectedPoints = GetAllConnectedPoints(Point); … … 2931 2960 2932 2961 if (connectedPoints == NULL) { 2933 Log() << Verbose(2) << "Could not find any connected points!" << endl;2962 eLog() << Verbose(2) << "Could not find any connected points!" << endl; 2934 2963 delete(connectedCircle); 2935 2964 return NULL; 2936 2965 } 2937 Log() << Verbose(2) << "Begin of GetCircleOfConnectedPoints" << endl;2938 2966 2939 2967 // calculate central point … … 2943 2971 // << "; scale factor " << 1.0/connectedPoints.size(); 2944 2972 center.Scale(1.0/connectedPoints->size()); 2945 Log() << Verbose( 4) << "INFO: Calculated center of all circle points is " << center << "." << endl;2973 Log() << Verbose(1) << "INFO: Calculated center of all circle points is " << center << "." << endl; 2946 2974 2947 2975 // projection plane of the circle is at the closes Point and normal is pointing away from center of all circle points … … 2949 2977 PlaneNormal.SubtractVector(¢er); 2950 2978 PlaneNormal.Normalize(); 2951 Log() << Verbose( 4) << "INFO: Calculated plane normal of circle is " << PlaneNormal << "." << endl;2979 Log() << Verbose(1) << "INFO: Calculated plane normal of circle is " << PlaneNormal << "." << endl; 2952 2980 2953 2981 // construct one orthogonal vector … … 2958 2986 } 2959 2987 if ((Reference == NULL) || (AngleZero.NormSquared() < MYEPSILON )) { 2960 Log() << Verbose( 4) << "Using alternatively " << *(*connectedPoints->begin())->node << " as angle 0 referencer." << endl;2988 Log() << Verbose(1) << "Using alternatively " << *(*connectedPoints->begin())->node << " as angle 0 referencer." << endl; 2961 2989 AngleZero.CopyVector((*connectedPoints->begin())->node); 2962 2990 AngleZero.SubtractVector(Point->node); … … 2967 2995 } 2968 2996 } 2969 Log() << Verbose( 4) << "INFO: Reference vector on this plane representing angle 0 is " << AngleZero << "." << endl;2997 Log() << Verbose(1) << "INFO: Reference vector on this plane representing angle 0 is " << AngleZero << "." << endl; 2970 2998 if (AngleZero.NormSquared() > MYEPSILON) 2971 2999 OrthogonalVector.MakeNormalVector(&PlaneNormal, &AngleZero); 2972 3000 else 2973 3001 OrthogonalVector.MakeNormalVector(&PlaneNormal); 2974 Log() << Verbose( 4) << "INFO: OrthogonalVector on plane is " << OrthogonalVector << "." << endl;3002 Log() << Verbose(1) << "INFO: OrthogonalVector on plane is " << OrthogonalVector << "." << endl; 2975 3003 2976 3004 // go through all connected points and calculate angle … … 2980 3008 helper.ProjectOntoPlane(&PlaneNormal); 2981 3009 double angle = GetAngle(helper, AngleZero, OrthogonalVector); 2982 Log() << Verbose( 3) << "INFO: Calculated angle is " << angle << " for point " << **listRunner << "." << endl;3010 Log() << Verbose(0) << "INFO: Calculated angle is " << angle << " for point " << **listRunner << "." << endl; 2983 3011 anglesOfPoints.insert(pair<double, TesselPoint*>(angle, (*listRunner))); 2984 3012 } … … 2989 3017 2990 3018 delete(connectedPoints); 2991 2992 Log() << Verbose(2) << "End of GetCircleOfConnectedPoints" << endl;2993 3019 2994 3020 return connectedCircle; … … 3003 3029 list<list<TesselPoint*> *> * Tesselation::GetPathsOfConnectedPoints(const TesselPoint* const Point) const 3004 3030 { 3031 Info FunctionInfo(__func__); 3005 3032 map<double, TesselPoint*> anglesOfPoints; 3006 3033 list<list<TesselPoint*> *> *ListOfPaths = new list<list<TesselPoint*> *>; … … 3047 3074 StartLine = CurrentLine; 3048 3075 CurrentPoint = CurrentLine->GetOtherEndpoint(ReferencePoint); 3049 Log() << Verbose( 3)<< "INFO: Beginning path retrieval at " << *CurrentPoint << " of line " << *CurrentLine << "." << endl;3076 Log() << Verbose(1)<< "INFO: Beginning path retrieval at " << *CurrentPoint << " of line " << *CurrentLine << "." << endl; 3050 3077 do { 3051 3078 // push current one 3052 Log() << Verbose( 3) << "INFO: Putting " << *CurrentPoint << " at end of path." << endl;3079 Log() << Verbose(1) << "INFO: Putting " << *CurrentPoint << " at end of path." << endl; 3053 3080 connectedPath->push_back(CurrentPoint->node); 3054 3081 3055 3082 // find next triangle 3056 3083 for (TriangleMap::iterator Runner = CurrentLine->triangles.begin(); Runner != CurrentLine->triangles.end(); Runner++) { 3057 Log() << Verbose( 3) << "INFO: Inspecting triangle " << *Runner->second << "." << endl;3084 Log() << Verbose(1) << "INFO: Inspecting triangle " << *Runner->second << "." << endl; 3058 3085 if ((Runner->second != triangle)) { // look for first triangle not equal to old one 3059 3086 triangle = Runner->second; … … 3062 3089 if (!TriangleRunner->second) { 3063 3090 TriangleRunner->second = true; 3064 Log() << Verbose( 3) << "INFO: Connecting triangle is " << *triangle << "." << endl;3091 Log() << Verbose(1) << "INFO: Connecting triangle is " << *triangle << "." << endl; 3065 3092 break; 3066 3093 } else { 3067 Log() << Verbose( 3) << "INFO: Skipping " << *triangle << ", as we have already visited it." << endl;3094 Log() << Verbose(1) << "INFO: Skipping " << *triangle << ", as we have already visited it." << endl; 3068 3095 triangle = NULL; 3069 3096 } … … 3080 3107 if ((triangle->lines[i] != CurrentLine) && (triangle->lines[i]->ContainsBoundaryPoint(ReferencePoint))) { // not the current line and still containing Point 3081 3108 CurrentLine = triangle->lines[i]; 3082 Log() << Verbose( 3) << "INFO: Connecting line is " << *CurrentLine << "." << endl;3109 Log() << Verbose(1) << "INFO: Connecting line is " << *CurrentLine << "." << endl; 3083 3110 break; 3084 3111 } … … 3094 3121 } while (CurrentLine != StartLine); 3095 3122 // last point is missing, as it's on start line 3096 Log() << Verbose( 3) << "INFO: Putting " << *CurrentPoint << " at end of path." << endl;3123 Log() << Verbose(1) << "INFO: Putting " << *CurrentPoint << " at end of path." << endl; 3097 3124 if (StartLine->GetOtherEndpoint(ReferencePoint)->node != connectedPath->back()) 3098 3125 connectedPath->push_back(StartLine->GetOtherEndpoint(ReferencePoint)->node); … … 3100 3127 ListOfPaths->push_back(connectedPath); 3101 3128 } else { 3102 Log() << Verbose( 3) << "INFO: Skipping " << *runner->second << ", as we have already visited it." << endl;3129 Log() << Verbose(1) << "INFO: Skipping " << *runner->second << ", as we have already visited it." << endl; 3103 3130 } 3104 3131 } … … 3118 3145 list<list<TesselPoint*> *> * Tesselation::GetClosedPathsOfConnectedPoints(const TesselPoint* const Point) const 3119 3146 { 3147 Info FunctionInfo(__func__); 3120 3148 list<list<TesselPoint*> *> *ListofPaths = GetPathsOfConnectedPoints(Point); 3121 3149 list<list<TesselPoint*> *> *ListofClosedPaths = new list<list<TesselPoint*> *>; … … 3131 3159 connectedPath = *ListRunner; 3132 3160 3133 Log() << Verbose( 2) << "INFO: Current path is " << connectedPath << "." << endl;3161 Log() << Verbose(1) << "INFO: Current path is " << connectedPath << "." << endl; 3134 3162 3135 3163 // go through list, look for reappearance of starting Point and count … … 3141 3169 if ((*CircleRunner == *CircleStart) && (CircleRunner != CircleStart)) { // is not the very first point 3142 3170 // we have a closed circle from Marker to new Marker 3143 Log() << Verbose( 3) << count+1 << ". closed path consists of: ";3171 Log() << Verbose(1) << count+1 << ". closed path consists of: "; 3144 3172 newPath = new list<TesselPoint*>; 3145 3173 list<TesselPoint*>::iterator CircleSprinter = Marker; … … 3157 3185 } 3158 3186 } 3159 Log() << Verbose( 3) << "INFO: " << count << " closed additional path(s) have been created." << endl;3187 Log() << Verbose(1) << "INFO: " << count << " closed additional path(s) have been created." << endl; 3160 3188 3161 3189 // delete list of paths … … 3179 3207 set<BoundaryTriangleSet*> *Tesselation::GetAllTriangles(const BoundaryPointSet * const Point) const 3180 3208 { 3209 Info FunctionInfo(__func__); 3181 3210 set<BoundaryTriangleSet*> *connectedTriangles = new set<BoundaryTriangleSet*>; 3182 3211 … … 3217 3246 return 0.; 3218 3247 } else 3219 Log() << Verbose( 2) << "Removing point " << *point << " from tesselated boundary ..." << endl;3248 Log() << Verbose(0) << "Removing point " << *point << " from tesselated boundary ..." << endl; 3220 3249 3221 3250 // copy old location for the volume … … 3247 3276 NormalVector.Zero(); 3248 3277 for (map<class BoundaryTriangleSet *, int>::iterator Runner = Candidates.begin(); Runner != Candidates.end(); Runner++) { 3249 Log() << Verbose( 3) << "INFO: Removing triangle " << *(Runner->first) << "." << endl;3278 Log() << Verbose(1) << "INFO: Removing triangle " << *(Runner->first) << "." << endl; 3250 3279 NormalVector.SubtractVector(&Runner->first->NormalVector); // has to point inward 3251 3280 RemoveTesselationTriangle(Runner->first); … … 3277 3306 smallestangle = 0.; 3278 3307 for (MiddleNode = connectedPath->begin(); MiddleNode != connectedPath->end(); MiddleNode++) { 3279 Log() << Verbose( 3) << "INFO: MiddleNode is " << **MiddleNode << "." << endl;3308 Log() << Verbose(1) << "INFO: MiddleNode is " << **MiddleNode << "." << endl; 3280 3309 // construct vectors to next and previous neighbour 3281 3310 StartNode = MiddleNode; … … 3305 3334 MiddleNode = EndNode; 3306 3335 if (MiddleNode == connectedPath->end()) { 3307 Log() << Verbose(1) << "CRITICAL: Could not find a smallest angle!" << endl;3308 exit(255);3336 eLog() << Verbose(0) << "CRITICAL: Could not find a smallest angle!" << endl; 3337 performCriticalExit(); 3309 3338 } 3310 3339 StartNode = MiddleNode; … … 3315 3344 if (EndNode == connectedPath->end()) 3316 3345 EndNode = connectedPath->begin(); 3317 Log() << Verbose( 4) << "INFO: StartNode is " << **StartNode << "." << endl;3318 Log() << Verbose( 4) << "INFO: MiddleNode is " << **MiddleNode << "." << endl;3319 Log() << Verbose( 4) << "INFO: EndNode is " << **EndNode << "." << endl;3320 Log() << Verbose( 3) << "INFO: Attempting to create triangle " << (*StartNode)->Name << ", " << (*MiddleNode)->Name << " and " << (*EndNode)->Name << "." << endl;3346 Log() << Verbose(2) << "INFO: StartNode is " << **StartNode << "." << endl; 3347 Log() << Verbose(2) << "INFO: MiddleNode is " << **MiddleNode << "." << endl; 3348 Log() << Verbose(2) << "INFO: EndNode is " << **EndNode << "." << endl; 3349 Log() << Verbose(1) << "INFO: Attempting to create triangle " << (*StartNode)->Name << ", " << (*MiddleNode)->Name << " and " << (*EndNode)->Name << "." << endl; 3321 3350 TriangleCandidates[0] = *StartNode; 3322 3351 TriangleCandidates[1] = *MiddleNode; … … 3324 3353 triangle = GetPresentTriangle(TriangleCandidates); 3325 3354 if (triangle != NULL) { 3326 eLog() << Verbose( 2) << "New triangle already present, skipping!" << endl;3355 eLog() << Verbose(0) << "New triangle already present, skipping!" << endl; 3327 3356 StartNode++; 3328 3357 MiddleNode++; … … 3336 3365 continue; 3337 3366 } 3338 Log() << Verbose( 5) << "Adding new triangle points."<< endl;3367 Log() << Verbose(3) << "Adding new triangle points."<< endl; 3339 3368 AddTesselationPoint(*StartNode, 0); 3340 3369 AddTesselationPoint(*MiddleNode, 1); 3341 3370 AddTesselationPoint(*EndNode, 2); 3342 Log() << Verbose( 5) << "Adding new triangle lines."<< endl;3371 Log() << Verbose(3) << "Adding new triangle lines."<< endl; 3343 3372 AddTesselationLine(TPS[0], TPS[1], 0); 3344 3373 AddTesselationLine(TPS[0], TPS[2], 1); … … 3355 3384 // prepare nodes for next triangle 3356 3385 StartNode = EndNode; 3357 Log() << Verbose( 4) << "Removing " << **MiddleNode << " from closed path, remaining points: " << connectedPath->size() << "." << endl;3386 Log() << Verbose(2) << "Removing " << **MiddleNode << " from closed path, remaining points: " << connectedPath->size() << "." << endl; 3358 3387 connectedPath->remove(*MiddleNode); // remove the middle node (it is surrounded by triangles) 3359 3388 if (connectedPath->size() == 2) { // we are done … … 3362 3391 break; 3363 3392 } else if (connectedPath->size() < 2) { // something's gone wrong! 3364 Log() << Verbose(1) << "CRITICAL: There are only two endpoints left!" << endl;3365 exit(255);3393 eLog() << Verbose(0) << "CRITICAL: There are only two endpoints left!" << endl; 3394 performCriticalExit(); 3366 3395 } else { 3367 3396 MiddleNode = StartNode; … … 3391 3420 if (maxgain != 0) { 3392 3421 volume += maxgain; 3393 Log() << Verbose( 3) << "Flipping baseline with highest volume" << **Candidate << "." << endl;3422 Log() << Verbose(1) << "Flipping baseline with highest volume" << **Candidate << "." << endl; 3394 3423 OtherBase = FlipBaseline(*Candidate); 3395 3424 NewLines.erase(Candidate); … … 3402 3431 delete(connectedPath); 3403 3432 } 3404 Log() << Verbose( 1) << count << " triangles were created." << endl;3433 Log() << Verbose(0) << count << " triangles were created." << endl; 3405 3434 } else { 3406 3435 while (!ListOfClosedPaths->empty()) { … … 3410 3439 delete(connectedPath); 3411 3440 } 3412 Log() << Verbose( 1) << "No need to create any triangles." << endl;3441 Log() << Verbose(0) << "No need to create any triangles." << endl; 3413 3442 } 3414 3443 delete(ListOfClosedPaths); 3415 3444 3416 Log() << Verbose( 1) << "Removed volume is " << volume << "." << endl;3445 Log() << Verbose(0) << "Removed volume is " << volume << "." << endl; 3417 3446 3418 3447 return volume; … … 3431 3460 list<BoundaryTriangleSet*> *Tesselation::FindTriangles(const TesselPoint* const Points[3]) const 3432 3461 { 3462 Info FunctionInfo(__func__); 3433 3463 list<BoundaryTriangleSet*> *result = new list<BoundaryTriangleSet*>; 3434 3464 LineMap::const_iterator FindLine; … … 3479 3509 map<int, int> * Tesselation::FindAllDegeneratedLines() 3480 3510 { 3511 Info FunctionInfo(__func__); 3481 3512 map<int, class BoundaryLineSet *> AllLines; 3482 3513 map<int, int> * DegeneratedLines = new map<int, int>; … … 3484 3515 // sanity check 3485 3516 if (LinesOnBoundary.empty()) { 3486 Log() << Verbose(1) << "Warning:FindAllDegeneratedTriangles() was called without any tesselation structure.";3517 eLog() << Verbose(2) << "FindAllDegeneratedTriangles() was called without any tesselation structure."; 3487 3518 return DegeneratedLines; 3488 3519 } … … 3500 3531 AllLines.clear(); 3501 3532 3502 Log() << Verbose( 1) << "FindAllDegeneratedLines() found " << DegeneratedLines->size() << " lines." << endl;3533 Log() << Verbose(0) << "FindAllDegeneratedLines() found " << DegeneratedLines->size() << " lines." << endl; 3503 3534 map<int,int>::iterator it; 3504 3535 for (it = DegeneratedLines->begin(); it != DegeneratedLines->end(); it++) 3505 Log() << Verbose( 2) << (*it).first << " => " << (*it).second << endl;3536 Log() << Verbose(0) << (*it).first << " => " << (*it).second << endl; 3506 3537 3507 3538 return DegeneratedLines; … … 3516 3547 map<int, int> * Tesselation::FindAllDegeneratedTriangles() 3517 3548 { 3549 Info FunctionInfo(__func__); 3518 3550 map<int, int> * DegeneratedLines = FindAllDegeneratedLines(); 3519 3551 map<int, int> * DegeneratedTriangles = new map<int, int>; … … 3543 3575 delete(DegeneratedLines); 3544 3576 3545 Log() << Verbose( 1) << "FindAllDegeneratedTriangles() found " << DegeneratedTriangles->size() << " triangles:" << endl;3577 Log() << Verbose(0) << "FindAllDegeneratedTriangles() found " << DegeneratedTriangles->size() << " triangles:" << endl; 3546 3578 map<int,int>::iterator it; 3547 3579 for (it = DegeneratedTriangles->begin(); it != DegeneratedTriangles->end(); it++) 3548 Log() << Verbose( 2) << (*it).first << " => " << (*it).second << endl;3580 Log() << Verbose(0) << (*it).first << " => " << (*it).second << endl; 3549 3581 3550 3582 return DegeneratedTriangles; … … 3557 3589 void Tesselation::RemoveDegeneratedTriangles() 3558 3590 { 3591 Info FunctionInfo(__func__); 3559 3592 map<int, int> * DegeneratedTriangles = FindAllDegeneratedTriangles(); 3560 3593 TriangleMap::iterator finder; 3561 3594 BoundaryTriangleSet *triangle = NULL, *partnerTriangle = NULL; 3562 3595 int count = 0; 3563 3564 Log() << Verbose(1) << "Begin of RemoveDegeneratedTriangles" << endl;3565 3596 3566 3597 for (map<int, int>::iterator TriangleKeyRunner = DegeneratedTriangles->begin(); … … 3621 3652 // erase the pair 3622 3653 count += (int) DegeneratedTriangles->erase(triangle->Nr); 3623 Log() << Verbose( 1) << "RemoveDegeneratedTriangles() removes triangle " << *triangle << "." << endl;3654 Log() << Verbose(0) << "RemoveDegeneratedTriangles() removes triangle " << *triangle << "." << endl; 3624 3655 RemoveTesselationTriangle(triangle); 3625 3656 count += (int) DegeneratedTriangles->erase(partnerTriangle->Nr); 3626 Log() << Verbose( 1) << "RemoveDegeneratedTriangles() removes triangle " << *partnerTriangle << "." << endl;3657 Log() << Verbose(0) << "RemoveDegeneratedTriangles() removes triangle " << *partnerTriangle << "." << endl; 3627 3658 RemoveTesselationTriangle(partnerTriangle); 3628 3659 } else { 3629 Log() << Verbose( 1) << "RemoveDegeneratedTriangles() does not remove triangle " << *triangle3660 Log() << Verbose(0) << "RemoveDegeneratedTriangles() does not remove triangle " << *triangle 3630 3661 << " and its partner " << *partnerTriangle << " because it is essential for at" 3631 3662 << " least one of the endpoints to be kept in the tesselation structure." << endl; … … 3636 3667 LastTriangle = NULL; 3637 3668 3638 Log() << Verbose(1) << "RemoveDegeneratedTriangles() removed " << count << " triangles:" << endl; 3639 Log() << Verbose(1) << "End of RemoveDegeneratedTriangles" << endl; 3669 Log() << Verbose(0) << "RemoveDegeneratedTriangles() removed " << count << " triangles:" << endl; 3640 3670 } 3641 3671 … … 3650 3680 void Tesselation::AddBoundaryPointByDegeneratedTriangle(class TesselPoint *point, LinkedCell *LC) 3651 3681 { 3652 Log() << Verbose(2) << "Begin of AddBoundaryPointByDegeneratedTriangle" << endl; 3653 3682 Info FunctionInfo(__func__); 3654 3683 // find nearest boundary point 3655 3684 class TesselPoint *BackupPoint = NULL; … … 3667 3696 return; 3668 3697 } 3669 Log() << Verbose( 2) << "Nearest point on boundary is " << NearestPoint->Name << "." << endl;3698 Log() << Verbose(0) << "Nearest point on boundary is " << NearestPoint->Name << "." << endl; 3670 3699 3671 3700 // go through its lines and find the best one to split … … 3700 3729 3701 3730 // create new triangle to connect point (connects automatically with the missing spot of the chosen line) 3702 Log() << Verbose( 5) << "Adding new triangle points."<< endl;3731 Log() << Verbose(2) << "Adding new triangle points."<< endl; 3703 3732 AddTesselationPoint((BestLine->endpoints[0]->node), 0); 3704 3733 AddTesselationPoint((BestLine->endpoints[1]->node), 1); 3705 3734 AddTesselationPoint(point, 2); 3706 Log() << Verbose( 5) << "Adding new triangle lines."<< endl;3735 Log() << Verbose(2) << "Adding new triangle lines."<< endl; 3707 3736 AddTesselationLine(TPS[0], TPS[1], 0); 3708 3737 AddTesselationLine(TPS[0], TPS[2], 1); … … 3711 3740 BTS->GetNormalVector(TempTriangle->NormalVector); 3712 3741 BTS->NormalVector.Scale(-1.); 3713 Log() << Verbose( 3) << "INFO: NormalVector of new triangle is " << BTS->NormalVector << "." << endl;3742 Log() << Verbose(1) << "INFO: NormalVector of new triangle is " << BTS->NormalVector << "." << endl; 3714 3743 AddTesselationTriangle(); 3715 3744 3716 3745 // create other side of this triangle and close both new sides of the first created triangle 3717 Log() << Verbose( 5) << "Adding new triangle points."<< endl;3746 Log() << Verbose(2) << "Adding new triangle points."<< endl; 3718 3747 AddTesselationPoint((BestLine->endpoints[0]->node), 0); 3719 3748 AddTesselationPoint((BestLine->endpoints[1]->node), 1); 3720 3749 AddTesselationPoint(point, 2); 3721 Log() << Verbose( 5) << "Adding new triangle lines."<< endl;3750 Log() << Verbose(2) << "Adding new triangle lines."<< endl; 3722 3751 AddTesselationLine(TPS[0], TPS[1], 0); 3723 3752 AddTesselationLine(TPS[0], TPS[2], 1); … … 3725 3754 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount); 3726 3755 BTS->GetNormalVector(TempTriangle->NormalVector); 3727 Log() << Verbose( 3) << "INFO: NormalVector of other new triangle is " << BTS->NormalVector << "." << endl;3756 Log() << Verbose(1) << "INFO: NormalVector of other new triangle is " << BTS->NormalVector << "." << endl; 3728 3757 AddTesselationTriangle(); 3729 3758 … … 3732 3761 if ((BTS->lines[i]->ContainsBoundaryPoint(BestLine->endpoints[0])) && (BTS->lines[i]->ContainsBoundaryPoint(BestLine->endpoints[1]))) { 3733 3762 if (BestLine == BTS->lines[i]){ 3734 Log() << Verbose(1) << "CRITICAL:BestLine is same as found line, something's wrong here!" << endl;3735 exit(255);3763 eLog() << Verbose(0) << "BestLine is same as found line, something's wrong here!" << endl; 3764 performCriticalExit(); 3736 3765 } 3737 3766 BTS->lines[i]->triangles.insert( pair<int, class BoundaryTriangleSet *> (TempTriangle->Nr, TempTriangle) ); … … 3740 3769 } 3741 3770 } 3742 3743 // exit3744 Log() << Verbose(2) << "End of AddBoundaryPointByDegeneratedTriangle" << endl;3745 3771 }; 3746 3772 … … 3752 3778 void Tesselation::Output(const char *filename, const PointCloud * const cloud) 3753 3779 { 3780 Info FunctionInfo(__func__); 3754 3781 ofstream *tempstream = NULL; 3755 3782 string NameofTempFile; … … 3764 3791 NameofTempFile.erase(npos, 1); 3765 3792 NameofTempFile.append(TecplotSuffix); 3766 Log() << Verbose( 1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";3793 Log() << Verbose(0) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n"; 3767 3794 tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc); 3768 3795 WriteTecplotFile(tempstream, this, cloud, TriangleFilesWritten); … … 3778 3805 NameofTempFile.erase(npos, 1); 3779 3806 NameofTempFile.append(Raster3DSuffix); 3780 Log() << Verbose( 1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";3807 Log() << Verbose(0) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n"; 3781 3808 tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc); 3782 3809 WriteRaster3dFile(tempstream, this, cloud); -
src/tesselation.hpp
r8725ed rf67b6e 65 65 #define DistanceMultiMap multimap <double, pair < PointMap::iterator, PointMap::iterator> > 66 66 #define DistanceMultiMapPair pair <double, pair < PointMap::iterator, PointMap::iterator> > 67 68 #define TesselPointList list <TesselPoint *> 67 69 68 70 /********************************************** declarations *******************************/ … … 192 194 ~CandidateForTesselation(); 193 195 194 TesselPoint *point;196 TesselPointList pointlist; 195 197 BoundaryLineSet *BaseLine; 196 198 Vector OptCenter; … … 218 220 void AddTesselationTriangle(); 219 221 void AddTesselationTriangle(const int nr); 220 void AddCandidateTriangle(CandidateForTesselation &CandidateLine);222 void AddCandidateTriangle(CandidateForTesselation CandidateLine); 221 223 void RemoveTesselationTriangle(class BoundaryTriangleSet *triangle); 222 224 void RemoveTesselationLine(class BoundaryLineSet *line); … … 227 229 void FindStartingTriangle(const double RADIUS, const LinkedCell *LC); 228 230 void FindSecondPointForTesselation(class TesselPoint* a, Vector Oben, class TesselPoint*& OptCandidate, double Storage[3], double RADIUS, const LinkedCell *LC); 229 void FindThirdPointForTesselation(Vector &NormalVector, Vector &SearchDirection, Vector &OldSphereCenter, class BoundaryLineSet *BaseLine, const class TesselPoint * const ThirdNode, CandidateList* &candidates, double *ShortestAngle, const double RADIUS, const LinkedCell *LC) const;231 void FindThirdPointForTesselation(Vector &NormalVector, Vector &SearchDirection, Vector &OldSphereCenter, CandidateForTesselation &CandidateLine, const class TesselPoint * const ThirdNode, const double RADIUS, const LinkedCell *LC) const; 230 232 bool FindNextSuitableTriangle(CandidateForTesselation &CandidateLine, BoundaryTriangleSet &T, const double& RADIUS, const LinkedCell *LC); 231 233 int CheckPresenceOfTriangle(class TesselPoint *Candidates[3]) const; … … 296 298 mutable PointMap::const_iterator InternalPointer; 297 299 298 bool HasOtherBaselineBetterCandidate(const BoundaryLineSet * const BaseRay, const TesselPoint * const OptCandidate, double ShortestAngle, double RADIUS, const LinkedCell * const LC) const;300 //bool HasOtherBaselineBetterCandidate(const BoundaryLineSet * const BaseRay, const TesselPoint * const OptCandidate, double ShortestAngle, double RADIUS, const LinkedCell * const LC) const; 299 301 }; 300 302 -
src/tesselationhelpers.cpp
r8725ed rf67b6e 8 8 #include <fstream> 9 9 10 #include "info.hpp" 10 11 #include "linkedcell.hpp" 11 12 #include "log.hpp" … … 15 16 #include "verbose.hpp" 16 17 17 double DetGet(gsl_matrix * const A, const int inPlace) { 18 double DetGet(gsl_matrix * const A, const int inPlace) 19 { 20 Info FunctionInfo(__func__); 18 21 /* 19 22 inPlace = 1 => A is replaced with the LU decomposed copy. … … 45 48 void GetSphere(Vector * const center, const Vector &a, const Vector &b, const Vector &c, const double RADIUS) 46 49 { 50 Info FunctionInfo(__func__); 47 51 gsl_matrix *A = gsl_matrix_calloc(3,3); 48 52 double m11, m12, m13, m14; … … 111 115 const double HalfplaneIndicator, const double AlternativeIndicator, const double alpha, const double beta, const double gamma, const double RADIUS, const double Umkreisradius) 112 116 { 117 Info FunctionInfo(__func__); 113 118 Vector TempNormal, helper; 114 119 double Restradius; 115 120 Vector OtherCenter; 116 Log() << Verbose(3) << "Begin of GetCenterOfSphere.\n";117 121 Center->Zero(); 118 122 helper.CopyVector(&a); … … 128 132 Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma))); 129 133 NewUmkreismittelpunkt->CopyVector(Center); 130 Log() << Verbose( 4) << "Center of new circumference is " << *NewUmkreismittelpunkt << ".\n";134 Log() << Verbose(1) << "Center of new circumference is " << *NewUmkreismittelpunkt << ".\n"; 131 135 // Here we calculated center of circumscribing circle, using barycentric coordinates 132 Log() << Verbose( 4) << "Center of circumference is " << *Center << " in direction " << *Direction << ".\n";136 Log() << Verbose(1) << "Center of circumference is " << *Center << " in direction " << *Direction << ".\n"; 133 137 134 138 TempNormal.CopyVector(&a); … … 154 158 TempNormal.Normalize(); 155 159 Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius); 156 Log() << Verbose( 4) << "Height of center of circumference to center of sphere is " << Restradius << ".\n";160 Log() << Verbose(1) << "Height of center of circumference to center of sphere is " << Restradius << ".\n"; 157 161 TempNormal.Scale(Restradius); 158 Log() << Verbose( 4) << "Shift vector to sphere of circumference is " << TempNormal << ".\n";162 Log() << Verbose(1) << "Shift vector to sphere of circumference is " << TempNormal << ".\n"; 159 163 160 164 Center->AddVector(&TempNormal); 161 Log() << Verbose( 0) << "Center of sphere of circumference is " << *Center << ".\n";165 Log() << Verbose(1) << "Center of sphere of circumference is " << *Center << ".\n"; 162 166 GetSphere(&OtherCenter, a, b, c, RADIUS); 163 Log() << Verbose(0) << "OtherCenter of sphere of circumference is " << OtherCenter << ".\n"; 164 Log() << Verbose(3) << "End of GetCenterOfSphere.\n"; 167 Log() << Verbose(1) << "OtherCenter of sphere of circumference is " << OtherCenter << ".\n"; 165 168 }; 166 169 … … 174 177 void GetCenterofCircumcircle(Vector * const Center, const Vector &a, const Vector &b, const Vector &c) 175 178 { 179 Info FunctionInfo(__func__); 176 180 Vector helper; 177 181 double alpha, beta, gamma; … … 186 190 beta = M_PI - SideC.Angle(&SideA); 187 191 gamma = M_PI - SideA.Angle(&SideB); 188 //Log() << Verbose( 3) << "INFO: alpha = " << alpha/M_PI*180. << ", beta = " << beta/M_PI*180. << ", gamma = " << gamma/M_PI*180. << "." << endl;192 //Log() << Verbose(1) << "INFO: alpha = " << alpha/M_PI*180. << ", beta = " << beta/M_PI*180. << ", gamma = " << gamma/M_PI*180. << "." << endl; 189 193 if (fabs(M_PI - alpha - beta - gamma) > HULLEPSILON) { 190 194 eLog() << Verbose(1) << "GetCenterofCircumcircle: Sum of angles " << (alpha+beta+gamma)/M_PI*180. << " > 180 degrees by " << fabs(M_PI - alpha - beta - gamma)/M_PI*180. << "!" << endl; … … 219 223 double GetPathLengthonCircumCircle(const Vector &CircleCenter, const Vector &CirclePlaneNormal, const double CircleRadius, const Vector &NewSphereCenter, const Vector &OldSphereCenter, const Vector &NormalVector, const Vector &SearchDirection) 220 224 { 225 Info FunctionInfo(__func__); 221 226 Vector helper; 222 227 double radius, alpha; … … 236 241 if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON) // acos is not unique on [0, 2.*M_PI), hence extra check to decide between two half intervals 237 242 alpha = 2.*M_PI - alpha; 238 //Log() << Verbose( 2) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << OldSphereCenter << " and resulting angle is " << alpha << "." << endl;243 //Log() << Verbose(1) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << OldSphereCenter << " and resulting angle is " << alpha << "." << endl; 239 244 radius = helper.Distance(&OldSphereCenter); 240 245 helper.ProjectOntoPlane(&NormalVector); 241 246 // check whether new center is somewhat away or at least right over the current baseline to prevent intersecting triangles 242 247 if ((radius > HULLEPSILON) || (helper.Norm() < HULLEPSILON)) { 243 //Log() << Verbose( 2) << "INFO: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << "." << endl;248 //Log() << Verbose(1) << "INFO: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << "." << endl; 244 249 return alpha; 245 250 } else { … … 264 269 double MinIntersectDistance(const gsl_vector * x, void *params) 265 270 { 271 Info FunctionInfo(__func__); 266 272 double retval = 0; 267 273 struct Intersection *I = (struct Intersection *)params; … … 284 290 285 291 retval = HeightA.ScalarProduct(&HeightA) + HeightB.ScalarProduct(&HeightB); 286 //Log() << Verbose( 2) << "MinIntersectDistance called, result: " << retval << endl;292 //Log() << Verbose(1) << "MinIntersectDistance called, result: " << retval << endl; 287 293 288 294 return retval; … … 304 310 bool existsIntersection(const Vector &point1, const Vector &point2, const Vector &point3, const Vector &point4) 305 311 { 312 Info FunctionInfo(__func__); 306 313 bool result; 307 314 … … 351 358 352 359 if (status == GSL_SUCCESS) { 353 Log() << Verbose( 2) << "converged to minimum" << endl;360 Log() << Verbose(1) << "converged to minimum" << endl; 354 361 } 355 362 } while (status == GSL_CONTINUE && iter < 100); … … 376 383 t2 = HeightB.ScalarProduct(&SideB)/SideB.ScalarProduct(&SideB); 377 384 378 Log() << Verbose( 2) << "Intersection " << intersection << " is at "385 Log() << Verbose(1) << "Intersection " << intersection << " is at " 379 386 << t1 << " for (" << point1 << "," << point2 << ") and at " 380 387 << t2 << " for (" << point3 << "," << point4 << "): "; 381 388 382 389 if (((t1 >= 0) && (t1 <= 1)) && ((t2 >= 0) && (t2 <= 1))) { 383 Log() << Verbose( 0) << "true intersection." << endl;390 Log() << Verbose(1) << "true intersection." << endl; 384 391 result = true; 385 392 } else { 386 Log() << Verbose( 0) << "intersection out of region of interest." << endl;393 Log() << Verbose(1) << "intersection out of region of interest." << endl; 387 394 result = false; 388 395 } … … 407 414 double GetAngle(const Vector &point, const Vector &reference, const Vector &OrthogonalVector) 408 415 { 416 Info FunctionInfo(__func__); 409 417 if (reference.IsZero()) 410 418 return M_PI; … … 418 426 } 419 427 420 Log() << Verbose( 4) << "INFO: " << point << " has angle " << phi << " with respect to reference " << reference << "." << endl;428 Log() << Verbose(1) << "INFO: " << point << " has angle " << phi << " with respect to reference " << reference << "." << endl; 421 429 422 430 return phi; … … 433 441 double CalculateVolumeofGeneralTetraeder(const Vector &a, const Vector &b, const Vector &c, const Vector &d) 434 442 { 443 Info FunctionInfo(__func__); 435 444 Vector Point, TetraederVector[3]; 436 445 double volume; … … 456 465 bool CheckLineCriteriaForDegeneratedTriangle(const BoundaryPointSet * const nodes[3]) 457 466 { 467 Info FunctionInfo(__func__); 458 468 bool result = false; 459 469 int counter = 0; … … 482 492 } 483 493 if ((!result) && (counter > 1)) { 484 Log() << Verbose( 2) << "INFO: Degenerate triangle is ok, at least two, here " << counter << ", existing lines are used." << endl;494 Log() << Verbose(1) << "INFO: Degenerate triangle is ok, at least two, here " << counter << ", existing lines are used." << endl; 485 495 result = true; 486 496 } … … 489 499 490 500 491 /** Sort function for the candidate list. 492 */ 493 bool SortCandidates(const CandidateForTesselation* candidate1, const CandidateForTesselation* candidate2) 494 { 495 Vector BaseLineVector, OrthogonalVector, helper; 496 if (candidate1->BaseLine != candidate2->BaseLine) { // sanity check 497 eLog() << Verbose(1) << "sortCandidates was called for two different baselines: " << candidate1->BaseLine << " and " << candidate2->BaseLine << "." << endl; 498 //return false; 499 exit(1); 500 } 501 // create baseline vector 502 BaseLineVector.CopyVector(candidate1->BaseLine->endpoints[1]->node->node); 503 BaseLineVector.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node); 504 BaseLineVector.Normalize(); 505 506 // create normal in-plane vector to cope with acos() non-uniqueness on [0,2pi] (note that is pointing in the "right" direction already, hence ">0" test!) 507 helper.CopyVector(candidate1->BaseLine->endpoints[0]->node->node); 508 helper.SubtractVector(candidate1->point->node); 509 OrthogonalVector.CopyVector(&helper); 510 helper.VectorProduct(&BaseLineVector); 511 OrthogonalVector.SubtractVector(&helper); 512 OrthogonalVector.Normalize(); 513 514 // calculate both angles and correct with in-plane vector 515 helper.CopyVector(candidate1->point->node); 516 helper.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node); 517 double phi = BaseLineVector.Angle(&helper); 518 if (OrthogonalVector.ScalarProduct(&helper) > 0) { 519 phi = 2.*M_PI - phi; 520 } 521 helper.CopyVector(candidate2->point->node); 522 helper.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node); 523 double psi = BaseLineVector.Angle(&helper); 524 if (OrthogonalVector.ScalarProduct(&helper) > 0) { 525 psi = 2.*M_PI - psi; 526 } 527 528 Log() << Verbose(2) << *candidate1->point << " has angle " << phi << endl; 529 Log() << Verbose(2) << *candidate2->point << " has angle " << psi << endl; 530 531 // return comparison 532 return phi < psi; 533 }; 501 ///** Sort function for the candidate list. 502 // */ 503 //bool SortCandidates(const CandidateForTesselation* candidate1, const CandidateForTesselation* candidate2) 504 //{ 505 // Info FunctionInfo(__func__); 506 // Vector BaseLineVector, OrthogonalVector, helper; 507 // if (candidate1->BaseLine != candidate2->BaseLine) { // sanity check 508 // eLog() << Verbose(1) << "sortCandidates was called for two different baselines: " << candidate1->BaseLine << " and " << candidate2->BaseLine << "." << endl; 509 // //return false; 510 // exit(1); 511 // } 512 // // create baseline vector 513 // BaseLineVector.CopyVector(candidate1->BaseLine->endpoints[1]->node->node); 514 // BaseLineVector.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node); 515 // BaseLineVector.Normalize(); 516 // 517 // // create normal in-plane vector to cope with acos() non-uniqueness on [0,2pi] (note that is pointing in the "right" direction already, hence ">0" test!) 518 // helper.CopyVector(candidate1->BaseLine->endpoints[0]->node->node); 519 // helper.SubtractVector(candidate1->point->node); 520 // OrthogonalVector.CopyVector(&helper); 521 // helper.VectorProduct(&BaseLineVector); 522 // OrthogonalVector.SubtractVector(&helper); 523 // OrthogonalVector.Normalize(); 524 // 525 // // calculate both angles and correct with in-plane vector 526 // helper.CopyVector(candidate1->point->node); 527 // helper.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node); 528 // double phi = BaseLineVector.Angle(&helper); 529 // if (OrthogonalVector.ScalarProduct(&helper) > 0) { 530 // phi = 2.*M_PI - phi; 531 // } 532 // helper.CopyVector(candidate2->point->node); 533 // helper.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node); 534 // double psi = BaseLineVector.Angle(&helper); 535 // if (OrthogonalVector.ScalarProduct(&helper) > 0) { 536 // psi = 2.*M_PI - psi; 537 // } 538 // 539 // Log() << Verbose(1) << *candidate1->point << " has angle " << phi << endl; 540 // Log() << Verbose(1) << *candidate2->point << " has angle " << psi << endl; 541 // 542 // // return comparison 543 // return phi < psi; 544 //}; 534 545 535 546 /** … … 543 554 TesselPoint* FindSecondClosestPoint(const Vector* Point, const LinkedCell* const LC) 544 555 { 556 Info FunctionInfo(__func__); 545 557 TesselPoint* closestPoint = NULL; 546 558 TesselPoint* secondClosestPoint = NULL; … … 553 565 for(int i=0;i<NDIM;i++) // store indices of this cell 554 566 N[i] = LC->n[i]; 555 Log() << Verbose( 2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;567 Log() << Verbose(1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl; 556 568 557 569 LC->GetNeighbourBounds(Nlower, Nupper); 558 //Log() << Verbose( 0) << endl;570 //Log() << Verbose(1) << endl; 559 571 for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++) 560 572 for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++) 561 573 for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) { 562 574 const LinkedNodes *List = LC->GetCurrentCell(); 563 //Log() << Verbose( 3) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;575 //Log() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl; 564 576 if (List != NULL) { 565 577 for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) { … … 597 609 TesselPoint* FindClosestPoint(const Vector* Point, TesselPoint *&SecondPoint, const LinkedCell* const LC) 598 610 { 611 Info FunctionInfo(__func__); 599 612 TesselPoint* closestPoint = NULL; 600 613 SecondPoint = NULL; … … 607 620 for(int i=0;i<NDIM;i++) // store indices of this cell 608 621 N[i] = LC->n[i]; 609 Log() << Verbose( 3) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;622 Log() << Verbose(1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl; 610 623 611 624 LC->GetNeighbourBounds(Nlower, Nupper); 612 //Log() << Verbose( 0) << endl;625 //Log() << Verbose(1) << endl; 613 626 for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++) 614 627 for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++) 615 628 for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) { 616 629 const LinkedNodes *List = LC->GetCurrentCell(); 617 //Log() << Verbose( 3) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;630 //Log() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl; 618 631 if (List != NULL) { 619 632 for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) { … … 626 639 distance = currentNorm; 627 640 closestPoint = (*Runner); 628 //Log() << Verbose( 2) << "INFO: New Nearest Neighbour is " << *closestPoint << "." << endl;641 //Log() << Verbose(1) << "INFO: New Nearest Neighbour is " << *closestPoint << "." << endl; 629 642 } else if (currentNorm < secondDistance) { 630 643 secondDistance = currentNorm; 631 644 SecondPoint = (*Runner); 632 //Log() << Verbose( 2) << "INFO: New Second Nearest Neighbour is " << *SecondPoint << "." << endl;645 //Log() << Verbose(1) << "INFO: New Second Nearest Neighbour is " << *SecondPoint << "." << endl; 633 646 } 634 647 } … … 640 653 // output 641 654 if (closestPoint != NULL) { 642 Log() << Verbose( 2) << "Closest point is " << *closestPoint;655 Log() << Verbose(1) << "Closest point is " << *closestPoint; 643 656 if (SecondPoint != NULL) 644 657 Log() << Verbose(0) << " and second closest is " << *SecondPoint; … … 656 669 Vector * GetClosestPointBetweenLine(const BoundaryLineSet * const Base, const BoundaryLineSet * const OtherBase) 657 670 { 671 Info FunctionInfo(__func__); 658 672 // construct the plane of the two baselines (i.e. take both their directional vectors) 659 673 Vector Normal; … … 666 680 Normal.VectorProduct(&OtherBaseline); 667 681 Normal.Normalize(); 668 Log() << Verbose( 4) << "First direction is " << Baseline << ", second direction is " << OtherBaseline << ", normal of intersection plane is " << Normal << "." << endl;682 Log() << Verbose(1) << "First direction is " << Baseline << ", second direction is " << OtherBaseline << ", normal of intersection plane is " << Normal << "." << endl; 669 683 670 684 // project one offset point of OtherBase onto this plane (and add plane offset vector) … … 683 697 Normal.CopyVector(Intersection); 684 698 Normal.SubtractVector(Base->endpoints[0]->node->node); 685 Log() << Verbose( 3) << "Found closest point on " << *Base << " at " << *Intersection << ", factor in line is " << fabs(Normal.ScalarProduct(&Baseline)/Baseline.NormSquared()) << "." << endl;699 Log() << Verbose(1) << "Found closest point on " << *Base << " at " << *Intersection << ", factor in line is " << fabs(Normal.ScalarProduct(&Baseline)/Baseline.NormSquared()) << "." << endl; 686 700 687 701 return Intersection; … … 696 710 double DistanceToTrianglePlane(const Vector *x, const BoundaryTriangleSet * const triangle) 697 711 { 712 Info FunctionInfo(__func__); 698 713 double distance = 0.; 699 714 if (x == NULL) { … … 712 727 void WriteVrmlFile(ofstream * const vrmlfile, const Tesselation * const Tess, const PointCloud * const cloud) 713 728 { 729 Info FunctionInfo(__func__); 714 730 TesselPoint *Walker = NULL; 715 731 int i; … … 755 771 void IncludeSphereinRaster3D(ofstream * const rasterfile, const Tesselation * const Tess, const PointCloud * const cloud) 756 772 { 773 Info FunctionInfo(__func__); 757 774 Vector helper; 758 775 … … 783 800 void WriteRaster3dFile(ofstream * const rasterfile, const Tesselation * const Tess, const PointCloud * const cloud) 784 801 { 802 Info FunctionInfo(__func__); 785 803 TesselPoint *Walker = NULL; 786 804 int i; … … 828 846 void WriteTecplotFile(ofstream * const tecplot, const Tesselation * const TesselStruct, const PointCloud * const cloud, const int N) 829 847 { 848 Info FunctionInfo(__func__); 830 849 if ((tecplot != NULL) && (TesselStruct != NULL)) { 831 850 // write header … … 848 867 849 868 // print atom coordinates 850 Log() << Verbose(2) << "The following triangles were created:";851 869 int Counter = 1; 852 870 TesselPoint *Walker = NULL; … … 858 876 *tecplot << endl; 859 877 // print connectivity 878 Log() << Verbose(1) << "The following triangles were created:" << endl; 860 879 for (TriangleMap::const_iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++) { 861 Log() << Verbose( 0) << " " << runner->second->endpoints[0]->node->Name << "<->" << runner->second->endpoints[1]->node->Name << "<->" << runner->second->endpoints[2]->node->Name;880 Log() << Verbose(1) << " " << runner->second->endpoints[0]->node->Name << "<->" << runner->second->endpoints[1]->node->Name << "<->" << runner->second->endpoints[2]->node->Name << endl; 862 881 *tecplot << LookupList[runner->second->endpoints[0]->node->nr] << " " << LookupList[runner->second->endpoints[1]->node->nr] << " " << LookupList[runner->second->endpoints[2]->node->nr] << endl; 863 882 } 864 883 delete[] (LookupList); 865 Log() << Verbose(0) << endl;866 884 } 867 885 }; … … 874 892 void CalculateConcavityPerBoundaryPoint(const Tesselation * const TesselStruct) 875 893 { 894 Info FunctionInfo(__func__); 876 895 class BoundaryPointSet *point = NULL; 877 896 class BoundaryLineSet *line = NULL; 878 897 879 //Log() << Verbose(2) << "Begin of CalculateConcavityPerBoundaryPoint" << endl;880 898 // calculate remaining concavity 881 899 for (PointMap::const_iterator PointRunner = TesselStruct->PointsOnBoundary.begin(); PointRunner != TesselStruct->PointsOnBoundary.end(); PointRunner++) { … … 885 903 for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) { 886 904 line = LineRunner->second; 887 //Log() << Verbose( 2) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;905 //Log() << Verbose(1) << "INFO: Current line of point " << *point << " is " << *line << "." << endl; 888 906 if (!line->CheckConvexityCriterion()) 889 907 point->value += 1; 890 908 } 891 909 } 892 //Log() << Verbose(2) << "End of CalculateConcavityPerBoundaryPoint" << endl;893 910 }; 894 911 … … 901 918 bool CheckListOfBaselines(const Tesselation * const TesselStruct) 902 919 { 920 Info FunctionInfo(__func__); 903 921 LineMap::const_iterator testline; 904 922 bool result = false; … … 908 926 for (testline = TesselStruct->LinesOnBoundary.begin(); testline != TesselStruct->LinesOnBoundary.end(); testline++) { 909 927 if (testline->second->triangles.size() != 2) { 910 Log() << Verbose( 1) << *testline->second << "\t" << testline->second->triangles.size() << endl;928 Log() << Verbose(2) << *testline->second << "\t" << testline->second->triangles.size() << endl; 911 929 counter++; 912 930 }
Note:
See TracChangeset
for help on using the changeset viewer.